{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Parallel Tempering with MPI" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Parallel tempering is a technique for improving the efficiency of MCMC, especially in the case of multi-modal target distributions. The core idea is to build a sequence of distributions starting with an easy-to-sample reference (e.g. the prior) and ending with the difficult-to-sample target (e.g. the posterior). Given the target $\\pi$ and the reference $\\pi_0$, this sequence is constructed using a \"temperature\" parameter $0 < \\beta < 1$:\n", "$$\\pi_\\beta = \\pi_0^{1-\\beta} \\cdot \\pi^\\beta$$\n", "Here, for $\\beta = 0$ (also referred to as \"cold\" chain) the distribution is equal to the target, for $\\beta = 1$ (also referred to as \"hot\" chain) it is equal to the reference. By running multiple chains with an ascending sequence of temperatures $\\beta_0 = 0, \\beta_1, \\dots, \\beta_n = 1$ and allowing theses chains to pass states based on Metropolis-Hastings, the simpler properties of hotter distributions improve the exploration of the state space and can result in faster convergence of the cold chain to the target distribution.\n", "\n", "*hopsy* implements parallel tempering. This notebook illustrates it by sampling a mixture of Gaussians that has two distinct modes. Depending on the starting point, vanilla MCMC approaches have trouble capturing the multi-modality. This is because once the chain has found a mode, Metropolis-Hastings proposal are very unlikely to propose high-density points in the other mode. With parallel tempering, the hotter chains are less effected by this and can better sample broadly. By passing these states on to the colder chains, other modes can be explored." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import hopsy\n", "import multiprocessing\n", "import ipyparallel as ipp\n", "\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The sampling code with all imports must be defined within a file or function like here. The code is later called by an MPI process, enabling communication of the tempered chains in *HOPS*. Notice, that for the definition of the Markov chain, a synchronized `RandomNumberGenerator` (must have the same seed for each process!) and an exchange probability must be given.\n", "\n", "Important: For technical reasons it is not possible to use multiple chains from each MPI process. A different way of generating replica chains for computing statistics is shown below. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def run_tempered_chain(n_samples: int, chain_idx: int):\n", " import hopsy\n", " import numpy as np\n", " from mpi4py import MPI\n", "\n", " class GaussianMixture:\n", " def __init__(self, mu1, mu2):\n", " epsilon = 0.05\n", " cov = epsilon * np.eye(2, 2)\n", " self.model1 = hopsy.Gaussian(mean=mu1, covariance=cov)\n", " self.model2 = hopsy.Gaussian(mean=mu2, covariance=cov)\n", "\n", " def compute_negative_log_likelihood(self, x):\n", " return -np.log(\n", " np.exp(-self.model1.compute_negative_log_likelihood(x))\n", " + np.exp(-self.model2.compute_negative_log_likelihood(x))\n", " )\n", " \n", " A = np.array([[1, 0], [0, 1], [-1, 0], [0, -1]])\n", " b = np.array([1, 1, 1, 1])\n", "\n", " model = GaussianMixture(np.ones(2).reshape(2, 1), -np.ones(2).reshape(2, 1))\n", " problem = hopsy.Problem(A, b, model)\n", "\n", " comm = MPI.COMM_WORLD\n", " rank = comm.Get_rank()\n", "\n", " syncRng = hopsy.RandomNumberGenerator(42)\n", " mc = hopsy.MarkovChain(\n", " proposal=hopsy.GaussianCoordinateHitAndRunProposal,\n", " problem=problem,\n", " parallel_tempering_sync_rng=syncRng,\n", " exchange_attempt_probability=0.15,\n", " starting_point=0.9 * np.ones(2))\n", " mc.proposal.stepsize = 0.25\n", "\n", " rng = hopsy.RandomNumberGenerator(rank + chain_idx + 11)\n", "\n", " acc_rate, samples = hopsy.sample(markov_chains=mc, rngs=rng, n_samples=n_samples)\n", " return acc_rate, samples, rank" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function implements running MPI processes from within notebooks using `ipyparallel`. An different strategy would be to place your code in a Python file and invoking the processes using `mpiexec` from the command line." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def run_tempered_chains(n_samples: int, n_temps: int, chain_idx: int):\n", " with ipp.Cluster(engines='mpi', n=n_temps) as rc:\n", " view = rc.broadcast_view()\n", " result = view.apply_sync(run_tempered_chain, n_samples, chain_idx)\n", " result = sorted(result, key=lambda x: x[2])\n", "\n", " return [r[1] for r in result]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because this approach is hard to debug, it is generally a good idea to run a single instance of tempered chains for a short time before running large numbers of samples." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting 5 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "eab8159d166941b4876eec905042b5c7", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/5 [00:00<?, ?engine/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Stopping engine(s): 1713946637\n", "engine set stopped 1713946637: {'exit_code': 1, 'pid': 48750, 'identifier': 'ipengine-1713946635-uenb-1713946637-48604'}\n", "Stopping controller\n", "Controller stopped: {'exit_code': 0, 'pid': 48691, 'identifier': 'ipcontroller-1713946635-uenb-48604'}\n" ] } ], "source": [ "_ = run_tempered_chains(10, 5, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using multiprocessing, multiple replicas of tempered chains can be run in parallel. This is simply done by using the `multiprocessing` module and calling the spawning function from before in different subprocesses." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting 5 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>\n", "Starting 5 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>\n", "Starting 5 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "afd31492153c41619b7bab68e7415bd5", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/5 [00:00<?, ?engine/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Starting 5 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ca10fd4dcdd941c1982a6b94ede97669", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/5 [00:00<?, ?engine/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "055259ac1f16495b8fd3d558e3ac1e6b", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/5 [00:00<?, ?engine/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "74255025e58d40aeb6244e48d247ecbf", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/5 [00:00<?, ?engine/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Stopping engine(s): 1713946647\n", "Stopping engine(s): 1713946647\n", "Stopping engine(s): 1713946647\n", "Stopping engine(s): 1713946647\n", "Stopping controller\n", "Stopping controller\n", "Stopping controller\n", "Stopping controller\n" ] } ], "source": [ "n_samples = 550000\n", "n_replicas = 4\n", "n_temps = 5\n", "\n", "with multiprocessing.Pool(n_replicas) as pool:\n", " result = pool.starmap(run_tempered_chains, [(n_samples, n_temps, chain_idx) for chain_idx in range(n_replicas)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assemble cold-chain samples and compute statistics:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([[1.0005234 , 1.00049931]]), array([[6946.92752069, 6910.3586694 ]]))" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "samples = np.array([result[i][0][0] for i in range(n_replicas)])\n", "hopsy.rhat(samples), hopsy.ess(samples)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cold chains seem to have converged. The following plot illustrates posteriors of all tempered distributions. The multi-modality is very well captured by all replicas." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAB8UAAAHqCAYAAACdjp8kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABxSUlEQVR4nO3deZhU5Zk47KcbpaFVUERoMC7gviAuhI67iUwQjYljFmV+o+gYl3whiYPGyIxLxCRoNARNHE2MipmMozEumYwGTYyMkwQxEpvEBeMCEZdGxbArCpzvjwwdq7ugu7qr6pyquu/rqkvr1Kmq53RR9ZzzPu9SlyRJEgAAAAAAAABQherTDgAAAAAAAAAASkVRHAAAAAAAAICqpSgOAAAAAAAAQNVSFAcAAAAAAACgaimKAwAAAAAAAFC1FMUBAAAAAAAAqFqK4gAAAAAAAABULUVxAAAAAAAAAKqWojgAAAAAAAAAVUtRHCiahQsXRl1dXcyYMaNt21e/+tWoq6tLLygAqHHyMwBki9wMANkjP0P1UxQHqs7KlSvj0ksvjWOOOSYGDBjQ4WQGACi/3/3udzFx4sTYZ599Yosttogdd9wxPvOZz8Sf/vSntEMDgJr01FNPxac//ekYPnx4NDY2xsCBA+OII46In/3sZ2mHBgD8n69//etRV1cX++67b9qhQMVTFAdK6qKLLoq33367rO/55ptvxpQpU+KZZ56JkSNHlvW9AaASpJGfr7zyyrjrrrvi6KOPjmuuuSbOOuuseOSRR+LAAw+MJ598sqyxAEDWpJGb//znP8eKFStiwoQJcc0118TFF18cEREf//jH4/vf/35ZYwGALEojP7/fyy+/HN/4xjdiiy22SC0GqCabpR0AUH6rVq0qWyLdbLPNYrPNyvtTM2TIkHjttdeiqakpHn/88fjgBz9Y1vcHgO6o9vw8adKkuO2226J3795t20466aQYMWJEXHHFFfGjH/2orPEAQGeqPTcfe+yxceyxx+ZsmzhxYhx00EExbdq0OOuss8oaDwB0RbXn5/c7//zz40Mf+lCsW7cu3nzzzdTigGphpDhUuQ3rnjz99NPxD//wD7HNNtvEYYcd1vb4j370ozjooIOib9++MWDAgDj55JNj0aJFOa9x1FFHxb777htz586NQw45JPr27RvDhg2LG264ocvv396PfvSjGD16dDQ2NsY222wTRxxxRDz44INtj//0pz+N4447LoYOHRoNDQ2xyy67xOWXXx7r1q3r9D0bGhqiqamp0/0AIC21mJ8POeSQnIJ4RMRuu+0W++yzTzzzzDOdPh8ASqkWc3M+vXr1ih122CGWLl3arecDQDHVcn5+5JFH4ic/+UlMnz69y88BNk1RHGrEpz/96Vi9enV84xvfiDPPPDMi/roeyamnnhq77bZbTJs2Lc4999x46KGH4ogjjuhwAfyXv/wljj322DjooIPim9/8ZnzgAx+Iz33uc3HzzTcXHMtll10Wp5xySmy++eYxZcqUuOyyy2KHHXaIX/3qV237zJgxI7bccsuYNGlSXHPNNXHQQQfFJZdcEhdeeGGP/g4AkCW1np+TJInFixfHwIEDu/V8ACi2WszNq1atijfffDNeeOGF+Pa3vx0///nP4+ijjy44XgAolVrLz+vWrYsvfOEL8dnPfjZGjBhRcIzARiRAVbv00kuTiEjGjx+fs33hwoVJr169kq9//es52//4xz8mm222Wc72I488MomI5Fvf+lbbtjVr1iT7779/MmjQoOTdd99NkiRJFixYkEREcsstt3R4/w2ee+65pL6+Pvn7v//7ZN26dTnvvX79+rb/X716dYdjOfvss5PGxsbknXfe6fLx/+53v+sQEwCkrdbz8wb//u//nkREctNNNxX8XAAoplrOzWeffXYSEUlEJPX19cmnPvWp5K233urScwGglGo1P3/3u99N+vfvn7z++uttx7DPPvt0+jxg04wUhxpxzjnn5Ny/++67Y/369fGZz3wm3nzzzbZbU1NT7LbbbvHwww/n7L/ZZpvF2Wef3Xa/d+/ecfbZZ8frr78ec+fO7XIc9957b6xfvz4uueSSqK/P/Ql6/1Q0ffv2bfv/FStWxJtvvhmHH354rF69OubPn9/l9wOALKvl/Dx//vz4/Oc/HwcffHBMmDChoOcCQKnUYm4+99xz4xe/+EXceuutMW7cuFi3bl28++67XY4VAEqtlvLzkiVL4pJLLomLL744tttuuy7HBnRus7QDAMpj2LBhOfefe+65SJIkdtttt7z7b7755jn3hw4dGltssUXOtt133z0iIhYuXBgf+tCHuhTHCy+8EPX19bH33ntvcr+nnnoqLrroovjVr34Vy5cvz3ls2bJlXXovAMi6Ws3Pra2tcdxxx0X//v3jJz/5SfTq1avLzwWAUqrF3LznnnvGnnvuGRERp556anz0ox+N448/PubMmZN3HVUAKLdays8XXXRRDBgwIL7whS90KSag6xTFoUa8v3daRMT69eujrq4ufv7zn+dtiN5yyy3LFVoHS5cujSOPPDL69esXU6ZMiV122SX69OkTv//97+MrX/lKrF+/PrXYAKCYajE/L1u2LMaNGxdLly6N//3f/42hQ4eWOHIA6LpazM3tfepTn4qzzz47/vSnP8Uee+xR5KgBoHC1kp+fe+65+P73vx/Tp0+PV199tW37O++8E++9914sXLgw+vXrFwMGDCjHoUDVURSHGrXLLrtEkiQxbNiwtl5xm/Lqq6/GqlWrcnrU/elPf4qIiJ133rmg912/fn08/fTTsf/+++fdZ9asWbFkyZK4++6744gjjmjbvmDBgi6/DwBUomrPz++8804cf/zx8ac//Sl++ctfdtq7HgDSVu25OZ+33347IszSBkB2VWt+fuWVV2L9+vXxxS9+Mb74xS92eHzYsGHxpS99KaZPn97lmIG/saY41KgTTzwxevXqFZdddlkkSZLzWJIksWTJkpxta9euje9973tt999999343ve+F9ttt10cdNBBXX7fE044Ierr62PKlCkdesVtiGND7773x/Xuu+/Gv/3bv3X5fQCgElVzfl63bl2cdNJJMXv27Ljzzjvj4IMP7nJ8AJCWas7Nr7/+eodt7733Xvzwhz+Mvn376rwGQGZVa37ed99945577ulw22effWLHHXeMe+65J84444wuxwvkMlIcatQuu+wSX/va12Ly5MmxcOHCOOGEE2KrrbaKBQsWxD333BNnnXVWnH/++W37Dx06NK688spYuHBh7L777nHHHXdES0tLfP/73++wRsum7LrrrvGv//qvcfnll8fhhx8eJ554YjQ0NMTvfve7GDp0aEydOjUOOeSQ2GabbWLChAnxxS9+Merq6uLf//3fO5zgbMp3v/vdWLp0ads0Mz/72c/i5ZdfjoiIL3zhC9G/f/8uvxYAlEs15+fzzjsv/uu//iuOP/74eOutt+JHP/pRzuP/+I//2OV4AaBcqjk3n3322bF8+fI44ogjYvvtt4/W1tb4j//4j5g/f35861vfSnXqWQDYlGrNzwMHDowTTjihw/YNI8PzPQZ0naI41LALL7wwdt999/j2t78dl112WURE7LDDDvHRj340Pv7xj+fsu80228Stt94aX/jCF+LGG2+MwYMHx3e/+90488wzC37fKVOmxLBhw+I73/lO/Ou//ms0NjbGfvvtF6ecckpERGy77bbx3//933HeeefFRRddFNtss0384z/+Yxx99NExduzYLr3H1VdfHX/+85/b7t99991x9913R8RfG90VxQHIqmrNzy0tLRHx145qP/vZzzo8rigOQFZVa24+6aST4qabborrr78+lixZEltttVUcdNBBceWVV3Y4LgDImmrNz0Dp1CWFDL0EatJRRx0Vb775Zjz55JNphwIA/B/5GQCyRW4GgOyRn4ENrCkOAAAAAAAAQNVSFAcAAAAAAACgaimKAwAAAAAAAFC1rCkOAAAAAAAAQNUyUhwAAAAAAACAqqUoDgAAAAAAAEDV2iztAIph/fr18eqrr8ZWW20VdXV1aYcDABUpSZJYsWJFDB06NOrre95vTn4GgJ6RmwEge+RnAMiWrubmqiiKv/rqq7HDDjukHQYAVIVFixbFBz7wgR6/jvwMAMUhNwNA9sjPAJAtneXmqiiKb7XVVhHx14Pt169fytEAQGVavnx57LDDDm15tafkZwDoGbkZALJHfgaAbOlqbq6KoviGaWX69evnxAEAeqhY07XJzwBQHHIzAGSP/AwA2dJZbu75oicAAAAAAAAAkFGK4gAAAAAAAABULUVxAAAAAAAAAKpWVawpDgAAAAAAQHVat25dvPfee2mHkYrevXtHfb0xrtBTiuIAAAAAAABkTpIk0draGkuXLk07lNTU19fHsGHDonfv3mmHAhVNURwAAAAAAIDM2VAQHzRoUDQ2NkZdXV3aIZXV+vXr49VXX43XXnstdtxxx5o7figmRXEAAAAAAAAyZd26dW0F8W233TbtcFKz3Xbbxauvvhpr166NzTffPO1woGIVvAjBI488Escff3wMHTo06urq4t57793k/qeddlrU1dV1uO2zzz5t+3z1q1/t8Piee+5Z8MEAAAAAAABQ+TasId7Y2JhyJOnaMG36unXrUo4EKlvBRfFVq1bFyJEj47rrruvS/tdcc0289tprbbdFixbFgAED4tOf/nTOfvvss0/Ofr/+9a8LDQ0AAAAAAIAqUutThtf68UOxFDx9+rhx42LcuHFd3r9///7Rv3//tvv33ntv/OUvf4nTTz89N5DNNoumpqZCwwEAAAAAAACAjSr7muI33XRTjBkzJnbaaaec7c8991wMHTo0+vTpEwcffHBMnTo1dtxxx7yvsWbNmlizZk3b/eXLl5c0ZgCgc/IzAGSL3AwA2SM/A0A6Cp4+vSdeffXV+PnPfx6f/exnc7Y3NzfHjBkzYubMmXH99dfHggUL4vDDD48VK1bkfZ2pU6e2jUDv379/7LDDDuUIHwDYBPkZALJFbgaA7JGfASAdZS2K33rrrbH11lvHCSeckLN93Lhx8elPfzr222+/GDt2bNx///2xdOnS+PGPf5z3dSZPnhzLli1ruy1atKgM0QMAmyI/A0C2yM0AkD3yM9SWGTNmxN577x2NjY2x1157xX333Zd2SFCzyjZ9epIkcfPNN8cpp5wSvXv33uS+W2+9dey+++7x/PPP5328oaEhGhoaShEmANBN8jMAZIvcDADZIz9Dccx5cUnZ3qt5+Lbdet5dd90VEydOjBtvvDGam5vj2muvjXPOOUdnGEhJ2UaK/8///E88//zzccYZZ3S678qVK+OFF16IIUOGlCEyAAAAAAAAKJ5p06bFeeedF+PHj4/hw4fHcccdt9Flg4HSK7govnLlymhpaYmWlpaIiFiwYEG0tLTESy+9FBF/nf7l1FNP7fC8m266KZqbm2Pfffft8Nj5558f//M//xMLFy6M3/72t/H3f//30atXrxg/fnyh4QEAAAAAAEBqVqxYEY8++mgce+yxbdseeOCBOOCAA1KMCmpbwdOnP/744/HhD3+47f6kSZMiImLChAkxY8aMeO2119oK5BssW7Ys7rrrrrjmmmvyvubLL78c48ePjyVLlsR2220Xhx12WDz66KOx3XbbFRoeAAAAAAAApGbevHlRX18fI0eOjNWrV8dtt90W1157bdxzzz1phwY1q+Ci+FFHHRVJkmz08RkzZnTY1r9//1i9evVGn3P77bcXGgYAAAAAAABkTktLS+y5554xd+7cOOywwyIi4sQTT4xx48ZFRMR///d/x3nnnRfr16+Pr3zlK/HZz342zXChJpRtTXEAYOPmvLgk5wYAAAAAVKaWlpY48MADY8SIETFnzpyYNm1azJw5M6ZMmRJr166NSZMmxa9+9at44okn4qqrroolS7QHQqkpigMAAAAAAECRbCiK9+vXL0aPHh3//M//HKecckrMmTMnHnvssdhnn31i++23jy233DLGjRsXDz74YNohQ9VTFAcAAAAAAIAiWLt2bTz11FOx11575WyfN29eHHbYYfHqq6/G9ttv37Z9++23j1deeaXcYULNKXhNcQAAAAAAAKCj+fPnxzvvvBNTpkyJ7bbbLhobG+P666+PhQsXxhlnnBG//vWv0w4RapKR4gAAAAAAAFAELS0tMWTIkOjbt28cfvjhccQRR8SiRYvi4Ycfjqamphg6dGjOyPBXXnklhg4dmmLEUBuMFAeADJrz4pKc+83Dt00pEgAAAADIliy3lbW0tERzc3Pcc889eR8fPXp0PPnkk/HKK69E//794+c//3lcfPHFZY4Sao+iOACkoH3RGwAAAACofC0tLXHooYdu9PHNNtssvvWtb8WHP/zhWL9+fVxwwQWx7bbZLfJDtVAUBwAAAAAAgCKYN29efO5zn9vkPh//+Mfj4x//eJkiAiIUxQEAAAAAAKAo3njjjbRDAPJQFAeAEjNVemnl+/tmeV0pAAAAAADKqz7tAAAAAAAAAACgVIwUBwAAAACAKtZ+ljUzrAFQaxTFAQAAAADIFEtlAQDFZPp0AAAAAAAAAKqWkeIAAABAUZmiFQAAgCxRFAcAao6GegAAAACA2qEoDgAAAAAANcSa7QDUGkVxAKCi5Ltw72wfF/YAAAAAALWrPu0AAAAAAAAAAKBUFMUBoALMeXFJhxsAAAAAkF0zZsyIvffeOxobG2OvvfaK++67L+2QoGaZPh0AAAAAAIDKsfA35XuvnQ/t1tPuuuuumDhxYtx4443R3Nwc1157bZxzzjmxaNGiIgcIdIWiOAAUmVHcAAAAAFDbpk2bFuedd16MHz8+IiKOO+64mDFjRrpBQQ1TFAcAAABKqn2nwebh227y8Xz7AABApVixYkU8+uijMW3atLZtDzzwQBxwwAEpRgWFK/RaLsvXcYriAEDV62z0voZ4ALqiki72AQCA9MybNy/q6+tj5MiRsXr16rjtttvi2muvjXvuuSft0KBm1acdAAAAAAAAAFSLlpaW2HPPPWPu3LmxxRZbxJlnnhnHH398jBs3LiIi/v7v/z622Wab+NSnPpVypFA7FMUBAAAAAACgSFpaWuLAAw+MESNGxJw5c2LatGkxc+bMmDJlSkREfOlLX4of/vCHKUcJtcX06QAAAAAAZJ6lTErL3xeKp6WlJU455ZTo169fjB49OkaPHh3PPvtszJkzJyIijjrqqJg1a1a6QUI3dLZMZZYpigMAAAAAADkUyaF71q5dG0899VTstddeOdvnzZsXH/vYx1KKClAUBwAAALqtOyMFKnl0AQAAbMr8+fPjnXfeiSlTpsR2220XjY2Ncf3118fChQvjjDPOSDs8qFmK4gAAAAAAlJVRyEC1amlpiSFDhkTfvn3j8MMPjy222CIOO+ywePjhh6OpqSnt8KBmKYoDAAAAAABQOXY+NO0INqqlpSWam5vjnnvuSTsUyISsdIRTFAcAAIB2TO8NAFQy5zKQnpaWljj00E0X7ceMGRPz5s2LVatWxQc+8IG488474+CDDy5ThFCbFMUBAAAAAACgCObNmxef+9znNrnPL3/5yzJFA2ygKF5k+XrgWQ8HAAAAAACg+r3xxhtphwDkoSgOAAAAAEDFy8qapQBQiap96Q1FcQAAAGqORvPK4zMDoL1qb7wHAIpHURwAAAC6QZEWAAAAKoOiOAAAAAAAAAA9kuVZXOrTDgAAAAAAAAAASsVI8QKZHg8AAAAAAACoZFke1V0KiuIAAABUvc4u9mutMQAAskYuBgBKSVEcAHrIhTsAAAAAAGSXojgAQB6WTAGAdOl4CAAAQLHUpx0AAAAAAAAAAJSKojgAAAAAAAAU2YwZM2LvvfeOxsbG2GuvveK+++5LOySoWaZPBwAAAAAAoGI83vp42d5rVNOobj3vrrvuiokTJ8aNN94Yzc3Nce2118Y555wTixYtKnKEQFcoineiszXMrHEGAAAAAADA+02bNi3OO++8GD9+fEREHHfccTFjxox0g4IaVvD06Y888kgcf/zxMXTo0Kirq4t77713k/vPmjUr6urqOtxaW1tz9rvuuuti5513jj59+kRzc3M89thjhYYGAAAAAAAAqVqxYkU8+uijceyxx7Zte+CBB+KAAw5IMSqobQUXxVetWhUjR46M6667rqDnPfvss/Haa6+13QYNGtT22B133BGTJk2KSy+9NH7/+9/HyJEjY+zYsfH6668XGh4AAAAAAACkZt68eVFfXx8jR46M1atXxw9+8IO49tpr4/zzz087NKhZBU+fPm7cuBg3blzBbzRo0KDYeuut8z42bdq0OPPMM+P000+PiIgbbrgh7rvvvrj55pvjwgsvLPi9AAAAAAAAIA0tLS2x5557xty5c+Owww6LiIgTTzwxxo0bF4sWLYpTTjklXn/99dhss83i4osvjk9/+tMpRwzVr+CR4t21//77x5AhQ+Lv/u7v4je/+U3b9nfffTfmzp0bY8aM+VtQ9fUxZsyYmD17drnCAwAAAACgROa8uCTnRuXxGULXtbS0xIEHHhgjRoyIOXPmxLRp02LmzJkxZcqU2GyzzWL69Onx9NNPx4MPPhjnnnturFq1Ku2QoeoVPFK8UEOGDIkbbrghRo0aFWvWrIkf/OAHcdRRR8WcOXPiwAMPjDfffDPWrVsXgwcPznne4MGDY/78+Xlfc82aNbFmzZq2+8uXLy/pMQBAqW21+LGc+ysGj04pku6TnwEgW+RmAMge+RlqQ0tLS5xyyinRr1+/GD16dIwePTqeffbZmDNnTnz1q1+NIUOGREREU1NTDBw4MN56663YYostUo4aqlvJR4rvsccecfbZZ8dBBx0UhxxySNx8881xyCGHxLe//e1uv+bUqVOjf//+bbcddtihiBGXnh51AFSjSs/PAFBt5GYAyB75Garf2rVr46mnnoq99torZ/u8efPaplLfYO7cubFu3Tq/BVAGZZs+/f1Gjx4dzz//fEREDBw4MHr16hWLFy/O2Wfx4sXR1NSU9/mTJ0+OZcuWtd0WLVpU8ph7QhEcgFpQafkZAKqd3AwA2SM/Q/WbP39+vPPOOzFlypR44okn4tlnn41zzz03Fi5cGGeccUbbfm+99Vaceuqp8f3vfz/FaKF2lHz69HxaWlrapobo3bt3HHTQQfHQQw/FCSecEBER69evj4ceeigmTpyY9/kNDQ3R0NBQrnABgC6QnwEgW+RmAMge+Rmq34YaWN++fePwww+PLbbYIg477LB4+OGH2waDrlmzJk444YS48MIL45BDDkk5YqgNBRfFV65c2TbKOyJiwYIF0dLSEgMGDIgdd9wxJk+eHK+88kr88Ic/jIiI6dOnx7Bhw2KfffaJd955J37wgx/Er371q3jwwQfbXmPSpEkxYcKEGDVqVIwePTqmT58eq1atitNPP70IhwgAxWXWDwAAAKDW5WsfaR6+bQqRUItGNY1KO4SNamlpiebm5rjnnnvyPp4kSZx22mnxkY98JE455ZQyRwe1q+Ci+OOPPx4f/vCH2+5PmjQpIiImTJgQM2bMiNdeey1eeumltsfffffdOO+88+KVV16JxsbG2G+//eKXv/xlzmucdNJJ8cYbb8Qll1wSra2tsf/++8fMmTNj8ODBPTk2AICiaX+x70IfILt0YAMAANLS0tIShx566EYf/81vfhN33HFH7LfffnHvvfdGRMS///u/x4gRI8oUIdSmgoviRx11VCRJstHHZ8yYkXP/ggsuiAsuuKDT1504ceJGp0sHAAAAAACArJs3b1587nOf2+jjhx12WKxfv76MEQERKa0pTi4jzwCq21aLH0s7BAAAACgbs7YAteyNN95IOwQgj/q0AwAAAAAAAACAUjFSHACoeu1H668YPDqlSAAAAAAAKDcjxQEAAAAAAACoWkaKAwBVxzruAAAAAABsYKQ4AAAAAAAAAFVLURwAAAAAAACAqmX6dACg5uWbbn3F4NEpRAJAVsgNAAAAUD0UxQGAmmPNcQAAAACA2mH6dAAAAAAAAACqlpHiAAAAAAAAAFVszotL0g4hVUaKAwAAAAAAPTbnxSU5N6h1M2bMiL333jsaGxtjr732ivvuuy/tkKBmGSkOAFSU9uuBrxg8OqVIAMgKDa61qf3n3jx825QiAQCg3Fb/7ndle6/GD36wW8+76667YuLEiXHjjTdGc3NzXHvttXHOOefEokWLihwhFEf7dteI6mp7VRQHgAxS+O26fCdr5aAhHgA2rrP87NwGAIBqN23atDjvvPNi/PjxERFx3HHHxYwZM9INCmqY6dMBAAAAAACgSFasWBGPPvpoHHvssW3bHnjggTjggANSjApqm5HiAAAAAAAAUCTz5s2L+vr6GDlyZKxevTpuu+22uPbaa+Oee+5JOzSoWYriANBDpjoHAACA4qr2dU3T5u8LpdXS0hJ77rlnzJ07Nw477LCIiDjxxBNj3LhxsXTp0hgzZkysXbs21q5dG1/60pfizDPPTDliqH6K4gAAANS8ztbABgAA6KqWlpY48MADY8SIETFnzpz4zW9+ExdddFFMmTIlLr744njkkUeisbExVq1aFfvuu2+ceOKJse2226YdNhSskgaMKYpn0JwXl3TY1jzcjyEAAAAAAOmopMIHpK2lpSVOOeWU6NevX4wePTpGjx4dzz77bMyZMyd69eoVjY2NERGxZs2aSJIkkiRJOWKoforiAAAAAABQ4xS9oTjWrl0bTz31VOy111452+fNmxcf+9jHIiJi6dKlceSRR8Zzzz0XV111VQwcODCNUKHosrw8h6I4AGxCvtk7OmP6VQCAwmS54QQAAAoxf/78eOedd2LKlCmx3XbbRWNjY1x//fWxcOHCOOOMMyIiYuutt4558+bF4sWL48QTT4xPfepTMXjw4JQjh85Vctu3ojgAQB56yANUt0q+kAcAALKrpaUlhgwZEn379o3DDz88tthiizjssMPi4Ycfjqamppx9Bw8eHCNHjoz//d//jU996lMpRQy1QVG8yPRuBwAAAADYtGJ0RNaZGWpX4wc/mHYIG9XS0hLNzc1xzz335H188eLF0djYGFtttVUsW7YsHnnkkfjc5z5X5iih9iiKA0AFyDuabfi48gcCAAAA3WCWlnR15+9fjM+s/bJ0zcO37fFrQta1tLTEoYceutHH//znP8dZZ50VSZJEkiTxhS98IUaMGFHGCKkV3VkatJopigMAAAAAkHkK60AlmDdv3iZHfo8ePTpaWlrKFxBsRK3lVUVxAAAAAABSVWsN80D1euONN9IOAchDURwAAAAoK4UPAAAAyklRPAXtL/5XDB6dUiQAAABQGTotpA8fV55AAAAAakQ1dWhWFAeAAlXTiQAA1Ar5GwBKZ86LS9IOISLy5HsdpgAgdVnJz4riAAAA0A1ZubBPmw4HAAAAZJ2ieIFKMfW56dQBAAAgl2I7AAAAxVKfdgAAAAAAAAAAUCpGipeB3u0AUP3yraHXPHzbFCIBoFRc2wEAAEBlUhQHAAAAKp4OagAAAGyMongPGSkAALWhfc5fMXh0SpEAAABkmzZTACBrFMU74QQOoLbkG2EEAAAAQK58bec6kAOQVfVpBwAAAACFmPPikpwbAABAFs2YMSP23nvvaGxsjL322ivuu+++tEOihm21+LGcW60xUhwAAACoeHkbdYaPK38gAACU3KvP/aVs7zV0t2269by77rorJk6cGDfeeGM0NzfHtddeG+ecc04sWrSoyBECXaEoXiHaj35oHr5tSpEAAABAZXAtDQBAWqZNmxbnnXdejB8/PiIijjvuuJgxY0a6QUENUxQHAOiG9qPRrJsGAADUqlqcghVgU1asWBGPPvpoTJs2rW3bAw88EAcccECKUUFtUxQHAAAAAIAqklZHBR3I4a/mzZsX9fX1MXLkyFi9enXcdtttce2118Y999yTdmhQs+rTDgAAAAAAAACqRUtLS+y5554xd+7c2GKLLeLMM8+M448/PsaNG9e2z+rVq2OnnXaK888/P8VIoXYoigMAAAAAAECRtLS0xIEHHhgjRoyIOXPmxLRp02LmzJkxZcqUtn2+/vWvx4c+9KEUo4TaYvp0AAAAAAAAKJKWlpY45ZRTol+/fjF69OgYPXp0PPvsszFnzpyIiHjuuedi/vz5cfzxx8eTTz6ZcrRUgzkvLkk7hMxTFAcAAKCipbVmJgAAQHtr166Np556Kvbaa6+c7fPmzYuPfexjERFx/vnnx1VXXRW//e1v0wgRapKiOABsgkZ2AAAAAKCr5s+fH++8805MmTIltttuu2hsbIzrr78+Fi5cGGeccUb89Kc/jd133z123313RXEoI0VxAAAAAACg6DoMNhg+Lp1AoIxaWlpiyJAh0bdv3zj88MNjiy22iMMOOywefvjhaGpqikcffTRuv/32uPPOO2PlypXx3nvvRb9+/eKSSy5JO3SoaoriAAAAAAAAVIyhu22Tdggb1dLSEs3NzXHPPffkfXzq1KkxderUiIiYMWNGPPnkkwriUAb1hT7hkUceieOPPz6GDh0adXV1ce+9925y/7vvvjv+7u/+Lrbbbrvo169fHHzwwfHAAw/k7PPVr3416urqcm577rlnoaEBAAAAAABAalpaWmK//fZLOwyIrRY/lnOrdQUXxVetWhUjR46M6667rkv7P/LII/F3f/d3cf/998fcuXPjwx/+cBx//PHxxBNP5Oy3zz77xGuvvdZ2+/Wvf11oaAAAAAAAAJCaefPmdbkoftppp8XVV19d4oiAiG5Mnz5u3LgYN67r635Mnz495/43vvGN+OlPfxo/+9nP4oADDvhbIJttFk1NTYWGU5Xy9dZYMXh0CpEAAAAAAADQVW+88UbaIQB5lH1N8fXr18eKFStiwIABOdufe+65GDp0aPTp0ycOPvjgmDp1auy44455X2PNmjWxZs2atvvLly8vacxZ0KFQPrzrHRMAoBxqMT8DQJbJzQCQPfIzAKSj4OnTe+rqq6+OlStXxmc+85m2bc3NzTFjxoyYOXNmXH/99bFgwYI4/PDDY8WKFXlfY+rUqdG/f/+22w477FCu8AGAjZCfASBb5GYAyJ5az89zXlyScwOAcilrUfy2226Lyy67LH784x/HoEGD2raPGzcuPv3pT8d+++0XY8eOjfvvvz+WLl0aP/7xj/O+zuTJk2PZsmVtt0WLFpXrEACAjZCfASgVjafdIzcDUCpyc/fJzwCQjrJNn3777bfHZz/72bjzzjtjzJgxm9x36623jt133z2ef/75vI83NDREQ0NDKcIEALpJfgagXDosL0VecjMAZI/8DIVLkiTtEFJV68cPxVKWovh//ud/xj/90z/F7bffHscdd1yn+69cuTJeeOGFOOWUU8oQHQC1TI92AAAAAMiezTffPCIiVq9eHX379k05mvS8++67ERHRq1evlCMhy3Qe71zBRfGVK1fmjOBesGBBtLS0xIABA2LHHXeMyZMnxyuvvBI//OEPI+KvU6ZPmDAhrrnmmmhubo7W1taIiOjbt2/0798/IiLOP//8OP7442OnnXaKV199NS699NLo1atXjB8/vhjHCABQcnlPPIePK38gAAAAAFWgV69esfXWW8frr78eERGNjY1RV1eXclTltX79+njjjTeisbExNtusbJM/Q1Uq+Bv0+OOPx4c//OG2+5MmTYqIiAkTJsSMGTPitddei5deeqnt8e9///uxdu3a+PznPx+f//zn27Zv2D8i4uWXX47x48fHkiVLYrvttovDDjssHn300dhuu+26e1wA0C161AEAAABANjQ1NUVEtBXGa1F9fX3suOOONdchAIqt4KL4UUcdtcn1CzYUujeYNWtWp695++23FxoGAAAAAABQQToMRjDDGp2oq6uLIUOGxKBBg+K9995LO5xU9O7dO+rr69MOAyqeuRYAAAAAAADIrF69ellTG+gRXUsAAAAAAAAAqFpGigMAVefp1S/m3N+7cXhKkQAAAAAAkDZF8Qo158UlOfebh2+bUiQAAAAAAAAA2WX6dAAAAAAAAACqlpHiAJAC03sDAJTeVosfy90wfFw6gQBUmQ6/rwBAxWvfZh1RXe3WiuIAQM3RKQEAAAAAqFTtl1neKqU4KomiOABAibQ/OW0evm1KkQBUNqPRAAAAgJ5QFE+B0WkAAAAAAGRJGu3W7TuTR+hQDkBpKIoDAAAAJdVZI3u1r10HAABQC7J8bacoDgA9ZAaQypflkzUAoHgsbQLQPdYtBQAiKrstvD7tAAAAAAAAAACgVIwUB4AiyzfqGADomnzrShqNBgAAAPSEoniF2mrxY7kbho9LJxAAiqKSp50BgEog1wJA93VoiywBuRoAKCVFcQAAAAAAoOzydrgwAAyAElAUBwAoETO7AAAAAACkT1G8SrRfd695+LYpRQIA2WOddwAAgPLpyjVYMaZHN+U6AJRWNbWrKooDAOTR2QmfxhaA6qJRHQCyJd81mfwMAH+Vd/kJNklRHAAqQN7GgBTiyKJq6q0IAJSWpU0ASsN1GQCQdYriAJBBGhQAAHIZzQ8AxaPdAYBaywWK4gAAAEBZ1VrjCwCkrRS51xT3AFQSRfEC6ZkOAAAAAFB+OlWVl7ZwAKqJongZdHay5uQCAACgeDSYVwefIwAAAMWiKA4AAADtmA4UAEpL5ycAoJwUxQGoaVstfiztEAAAAAAyT0cGgMpWit/xSsoNiuI9VEkfNgDF4bcfANJlCarsc74EULnmvLikw7atUogDAKCYFMUBqCntL+47u7DXoAsA1UmOBwAAgNqhKF4lOkz/O3xcOoEAQIkpYgBAtqSVm80YAAAA1IqnZ/887RAqnqJ4Jyq14T3fNEfNw7dNIRIAqE7daYhvn5/lZoDiqNTrNgAAACgV18q5FMUzqCv/SPWABwAAAAAAAOicojgAAAB0gem6sy1fB3OfEQCUVzHOl8yyBlBdOuSGlOJQFAcAAIBuyMqFPT2w8De593c+NJ04ACqM6VjpKp0KAcgKRXEAaspWix9LOwQAAAAAAKCMFMUBAADIDB3YAAAAoDBmcemcojgAvI+TB0qpQ6Fn+Lh0AgEAAKCqad8AgFz1aQcAAAAAAAAAAKVipDgAQIm075m/d+PwlCIBgNrUWS5++rXluY/vXOqIACilDr/7KcUBAGSPojgAQBGYmg4AAAAAIJsUxSuUkWcAAAAAQLFttfixtEOgxnX4Nzh8XDqBAFBVFMUBqFpzXlzSYdtWKcQBAAAAWdX+2tl1MwBkn1krC6coDgBQJta3A+hIQzwAAABQaoriAAAAZJoe8AAAAEBP1KcdAAAAAAAAAACUipHiAAAApGarxY+lHQI1pMNSJo3DU4oEAACAclIUB6CmmY6VNLVfR7d5+LYpRQIAANQqHdQoJx3UALpHO3bPKYpXqbwns8PHlT8QAAAAAAAAgBRZUxwAAAAAAACAqlVwUfyRRx6J448/PoYOHRp1dXVx7733dvqcWbNmxYEHHhgNDQ2x6667xowZMzrsc91118XOO+8cffr0iebm5njsMdP2ANAzWy1+rMMNAAAAKMzTq1/MuQEAVJqCi+KrVq2KkSNHxnXXXdel/RcsWBDHHXdcfPjDH46WlpY499xz47Of/Ww88MADbfvccccdMWnSpLj00kvj97//fYwcOTLGjh0br7/+eqHhAQAAAHTPwt90vAEAAFDxCl5TfNy4cTFuXNfXpr7hhhti2LBh8a1vfSsiIvbaa6/49a9/Hd/+9rdj7NixERExbdq0OPPMM+P0009ve859990XN998c1x44YWFhliT2vfQ3LtxeEqRAAAAAAAAAGRHwUXxQs2ePTvGjBmTs23s2LFx7rnnRkTEu+++G3Pnzo3Jkye3PV5fXx9jxoyJ2bNnlzo8AIDUdJjSf3jXOx4CVCtTsgIAAFDzzFpVdCUvire2tsbgwYNztg0ePDiWL18eb7/9dvzlL3+JdevW5d1n/vz5eV9zzZo1sWbNmrb7y5cvL37gAEBB5GcAyBa5GQCyR37Ola9DpFlQASiFgtcUz4KpU6dG//7922477LBD2iEBQM2TnwEgW+RmALrj6dUvdrhRPPIzAKSj5EXxpqamWLx4cc62xYsXR79+/aJv374xcODA6NWrV959mpqa8r7m5MmTY9myZW23RYsWlSx+AKqLC/vSkZ8BIFvkZgDIHvkZANJR8unTDz744Lj//vtztv3iF7+Igw8+OCIievfuHQcddFA89NBDccIJJ0RExPr16+Ohhx6KiRMn5n3NhoaGaGhoKGncAEBh5GcAusS6aGUjNwNA9sjP3dD+/HHnQ9OJA4CKVnBRfOXKlfH888+33V+wYEG0tLTEgAEDYscdd4zJkyfHK6+8Ej/84Q8jIuKcc86J7373u3HBBRfEP/3TP8WvfvWr+PGPfxz33Xdf22tMmjQpJkyYEKNGjYrRo0fH9OnTY9WqVXH66acX4RABAAAAAIBK9PRrueuu771zOnEAUNkKLoo//vjj8eEPf7jt/qRJkyIiYsKECTFjxox47bXX4qWXXmp7fNiwYXHffffFP//zP8c111wTH/jAB+IHP/hBjB07tm2fk046Kd5444245JJLorW1Nfbff/+YOXNmDB48uCfHBgAAALBRHZbTeW14h300vANAebXPz3s3dszPAFCogoviRx11VCRJstHHZ8yYkfc5TzzxxCZfd+LEiRudLh0AAAAAAACgFrSfJYOeK/ma4gAA5Neh93tKcQAAAAAAVDNF8SrVYQq40NAOAAAAAAAA1B5F8Vqy8De593c+NJ04AACAmvX40udy7jfG4JQigY50MAcAAKhOiuIAVK18jZoAAAAAAJBl2raLT1EcgOrRfkYMqDT5/g2b2QWoci70AQAAgFJTFAcAAAAAqFE6qAFAygz2KgtF8Rry9GvLc+7vvXM6cQAAALVr4ZJVaYcAAAAA1Jj6tAMAAAAAAAAAgFJRFAcAAAAAAACgapk+HQAAAAAAyKT2697vnVIcAFQ2RfEa4uQBAAAAAAAAsuPxpc+lHUJNUBQHAMiIfCfAo+LQFCIBADaY8+KSnPvNw7dNKRIAICLi6dk/z7m/98HjUooEgEpiTXEAAAAAAAAAqpaR4gBUjadfW552CAAAVJmtFj+Wu2G40WgAAACVRlEcAAAAAKBGPN76eNohAACUnaI4AAAAwEY8vfrFnPt7pxQHQNG0Ppl2BAAAZacoDgAAAAAAVAQd1gDoDkVxAAAASsYUrQCQLQuXrEo7BACAslMUB6BqtO8pDABkgClaAQAAgJQpigMAAFAyRqMBAAAAaVMUr2HtpzEc1TQqpUgAgIj8hSPZGQAAAACql87k5VGfdgAAAAAAAAAAUCpGitey9mv7GSkOAJljZhcAAKBHFv4m7QgAgPeTm1NhpDgAAAAAAAAAVctIcQAossbnXsm5v3q37VOKBAAAAAAAOte+XTuiutq2FcUBoECK3tUnyyd8z7y2POf+qKaUAgHogvZLPgAA6Xt86XNphwAAkDpFcQAqVqU0vGe54Er29X3rmXZbPpJKHAAAAABAbammAWKK4gBUrtYn044AAIASqabGFwA6yteB/P387gNQrbIyi0tnubjaKIrTpv2Iy1FNo1KKBABKqysnfBriAQAAoAIs/E3HbTsfWv44ADJG+2YuRXEAeB8nCgAAAACVI9+Iy1GhKA5ALkXxGrZwyaqc+zs3pRQIAAAAdCIznRfbj0YzEg3IuPZtgOWQb3Yunc4BoHS6MzNmMd6nkvK7ojgApKCSTx4AoBLJvdni8wAAAKCcFMUBqFhp9HbPpxg97Kg8XfncNfADAAC1oFyj0wCAbKmkDs+K4gAAAEDFSavxpf26pdYsBbLk8dbH0w6BGldJxRGAtHRnsJfOZT2nKJ4BmTlRaH0y937TqHTiAMgQJxv0RGZyPAAAQBG5VgYAKo2iOAAV4/GWGWmHQA3T6APQBe072kIXFSPP5nsNHdIAAACIUBTPpKxcyOebbmmU0eMAZWGEcXFVU0G7fX6WmwGoBNWUiwGqXVbaJrMSB3+jrQKgeyxvkg2K4mXgZAEAAKhGWbmw7841l+u02uRzB4CNkycBqGaK4rRZuGRVzv2dm1IKBKBI9ConazQwAAAAQHG1b9eOiDCnGpAplhrLBEXxAmnMBoDsM0UrAAAAlJZrbwAqiaJ4D2U18RsdCdA1Wfkd7yyOvI8fXaJgACDjupU3qUml+LfQYTRau2UERjUZmwZUH7mVStN+mR/5Gag0lZp7uxR3Su3aiuIVwgh1AACA8tHRGICKZHpWMka7NgBZoSheQ5yAAEAVaN/Ipbc7UGUUo6tDpY5qACDbOssvzhmKSz4HoJooiqegFCcTTlCAatN+mivYQCcvgMrmdxwAIH1ptSc7FwQgLYriANQUnYiqj88UIDs6rLVcJOX4rZdPAIBqltVzHbMEAVAuiuJsXL41iEzRCgBFo4c8AKSrKLnY0iYAkD3yMwDtKIoDUDWy2usZukoPeYBs05kJALKtwzXV0enEAQBkj6J4hVL4AQAAAACg2uiICNQaNb/yqO/Ok6677rrYeeedo0+fPtHc3ByPPfbYRvc96qijoq6ursPtuOOOa9vntNNO6/D4Mccc053QAABqyuOtj+fcAIDua3zulQ43gGrnd49a4NoZgIJHit9xxx0xadKkuOGGG6K5uTmmT58eY8eOjWeffTYGDRrUYf+777473n333bb7S5YsiZEjR8anP/3pnP2OOeaYuOWWW9ruNzQ0FBoaAAAAFUbje+XzGQIAAJB1BRfFp02bFmeeeWacfvrpERFxww03xH333Rc333xzXHjhhR32HzBgQM7922+/PRobGzsUxRsaGqKpqanQcACoIQuXrEo7BAAAaowpXAGoRDqtAUCugori7777bsydOzcmT57ctq2+vj7GjBkTs2fP7tJr3HTTTXHyySfHFltskbN91qxZMWjQoNhmm23iIx/5SHzta1+LbbfdtpDwAAAAKKfWJ3Pu5mt8VUAEgNLpzjTQiqVUGh3UgGojF6ejoKL4m2++GevWrYvBgwfnbB88eHDMnz+/0+c/9thj8eSTT8ZNN92Us/2YY46JE088MYYNGxYvvPBC/Mu//EuMGzcuZs+eHb169erwOmvWrIk1a9a03V++fHkhhwEAlID8DADZIjcDVB9rIVc++RkA0lHw9Ok9cdNNN8WIESNi9OjROdtPPvnktv8fMWJE7LfffrHLLrvErFmz4uijj+7wOlOnTo3LLrus5PFG6K0BkGV+o7OlnPkZgMohX6dHbk5HvoLVqKZRKUQCQBbJzwCQjvpCdh44cGD06tUrFi9enLN98eLFna4HvmrVqrj99tvjjDPO6PR9hg8fHgMHDoznn38+7+OTJ0+OZcuWtd0WLVrU9YMAAEpCfi6NxudeybkBpOnx1sdzbmSb3AwA2SM/95zrZCDrXDtnU0EjxXv37h0HHXRQPPTQQ3HCCSdERMT69evjoYceiokTJ27yuXfeeWesWbMm/vEf/7HT93n55ZdjyZIlMWTIkLyPNzQ0RENDQyGhAwAlJj8DUGs6NMJ2nOgsVXJzedQ/k1vMWN+0b0qRAHT8TYqIiEEDyh8IGyU/A0A6Cp4+fdKkSTFhwoQYNWpUjB49OqZPnx6rVq2K008/PSIiTj311Nh+++1j6tSpOc+76aab4oQTTohtt902Z/vKlSvjsssui09+8pPR1NQUL7zwQlxwwQWx6667xtixY3twaABUG71/AQAAAACoJHk7rVF2BRfFTzrppHjjjTfikksuidbW1th///1j5syZMXjw4IiIeOmll6K+PndW9meffTZ+/etfx4MPPtjh9Xr16hV/+MMf4tZbb42lS5fG0KFD46Mf/WhcfvnlesxlUPtpHqyLBpRM65NpRwCZ9Pqjf8q5P+gEo9EAoJw6dNTs3zedQACAbtPODZSUtu1MKrgoHhExceLEjU6XPmvWrA7b9thjj0iSJO/+ffv2jQceeKA7YdBD7S/kV++2fc79hUtWdXjOzpteOh6g26ytAgCVR293AAAAoBJ0qyhOdeqsSA4AZFD7nqd6twNFpNMaFIfRaECaLEUGAKAozibkPWHep/xxAABdl6+ApeEdKBpTwAEAUAHaz2i0vsnSYwC1TlEcAAAAAKAS6bAGeQd3tba7P6g8oQBEhKXGskpRHAAAAKCbWpe9nXN/UL4CldFpQIlodAcA6Jr6tAMAAAAAAAAAgFIxUpzCtO/xbo1SoEjq583Pub9+cEqBQIXpMDJkcZ7TO/kaAACAWtahXdssLkBxPN76eIdtRiRnk6I4BWnf8P54U+6XfZRGdwAAqBntp40G8mvfAVSHNaC7OvyeAACZ59o5GxTFAcgk66IBAABQ69qPPuts5JlGdwCA/BTFAUhFoRf2QBe9tSDtCAAAAAAAMkVRnB4xBRwAAAD8Td4ZjwYMK38gAAAAtFEUByAbjG6Fkmk/M8MondgAAAAAoGDt29moHIriAGSSddAAIHvyjoAFAAAAyDhLuAIAAAAAAABQtYwUBwAAAAAAAOhE/bz5OffXD16bUiQUSlGcolr9u9/l3G/84AdTigQAAAAAKlv7hnegNPKtETyqaVQKkQCVxjJjlUNRnKKa/9azOfcPDEVxAACoBL9/4EcdN47cs/yBQA1o/307cOw/phQJANSmDiM9nfcCVD1FcQrSuuztnPtNA1IKBACIiDy5uX/flCIBACI65uYI184AAABpUxSnZ95akHt/wLB04gAqjingoHw6fN9MAQcA5eXaGSiRfB1xAIB0yc/ZpCgOAABAXjqxAQBQDTpb87d+Xp6NOpQDVBVFcQDKYvXvfpd2CMD/af99bPzgB1OKBEiT3AwAAADUivq0AwAAAAAAAACAUjFSHICymP/Ws2mHAPyf9t/HA8NIcQAASFuXZnF5a0HpA4Eq1H5936b+fTt9zuOtj+fcH2U6dYCKpigOQNGZjhUAALrO0iYAAAClpShOSeUrjLm4BwCACtFuNFr7ETYAAABQzQwAqx7WFAcAAAAAAACgahkpDgAAAAAAANQUsx3XFkVxSmr+W8922HZg+EEBAAAAACC76ufNz93QNCqdQAAoCkVxAAAAgHJ6a0Hu/QF7pBMHAACwSfkGf1KZFMUpu/bTUZiKAipfvmlmAAAAgK7T6A4AUDqK4qROkRwAiqd12dsdtjUN2PRz5GKoTRreASBdXepg3n5mCQCgpAwAq16K4gAAAAAAAACd0WGtYimKAwBUuw7rlg5LJw6grAqeBcKFPWSamV0AoITynQu7dgaoKoriAAAAAAAAm6CDGkBlUxQHAAAAyBDrGAIAABSXojjFZXpWoKtM0QqZMf+tZ3PuHxh6u0Ol6UoBTZENAAAAqFWK4gAUXfsCW1e0Lnu7BJEA3ZGvcGZaOAAAKIyplgEAskNRHICCGWkGtUeDHlS+7nRaA4qjfQfQpv59c+63/37uOWCPkscEABSXDuaQfdq1a5uiOAAAAABAiWmIB4Ds06G8eimKAwAAAGRIvoY4o8cBIFt0dAEiLAtaSRTFAQAAACqMpU2g8nVnJJqGdyiOzpY2ycdyJwCVTVEcgPJ4a0HaEQAAAAAAADVIUZyy06MOAAAAAACADV597i8594futk1KkVCtFMUB6JQ1kqDKtJ+5YcCwdOIAAAAAACgDRXEAAACAjDPrGmRL+9FsWy99PqVIAICNMdiL91MUp+yWv9k7d8OA3Lvtf6QaP/jBEkcEABRKvoZsaX1tbYdtTUN6frnXuuztHr8GAFSjd56Zn7uhC3m3fecWACBdcnNtURQn8/L15NHwDgClk++CoNDRaO1HzkRYCwqKSW93AKBU3l65dYdtfbdcusl92j8OQG2p5fXAu5I3yQZFcagBtZyQKA896qB0stLY1H4UajFGoEK5pXVO5FwMoDL5/c6V1b9HVzpjZjV2AKD7yjYg460FBT8lK+155NKaSY8sX9ov536/rZcX/T3yTQU5vOjvAuVX6xfl7X8/os70rFAM+XqnRv81OXc75O92S5lUslr/ba1mWflssxJHKXTIzRHyM5BZhf4eV9MsNsU4lmr6e5RKKXJ++yUF+w18N3eHbjS6d4eGesiOfL/H7VXq73N3ck01X29VMp9L8eVtv2tHfi6+bhXFr7vuurjqqquitbU1Ro4cGd/5zndi9OjRefedMWNGnH766TnbGhoa4p133mm7nyRJXHrppXHjjTfG0qVL49BDD43rr78+dtttt+6ERwm1/6I2btHJE/KezOdOv5rVkWdpnZB09r7FeM9yXfx2J1mWIsF252/aWRxdiTMrJwtZiaNQpp2h0tT6v9mu5M1qVozf2qz8XmclDjYtX8fRUtBoDmXS/tp5wLBOn9L+d6B953FLkVGpCr0e78prVJPOcnOtX5dQG8pWTCowP+c7R29ql4/lYjallq7Hs5Kr23dYi8jTaY2qUXD18Y477ohJkybFDTfcEM3NzTF9+vQYO3ZsPPvsszFo0KC8z+nXr188++zfptatq6vLefyb3/xmXHvttXHrrbfGsGHD4uKLL46xY8fG008/HX369Ck0RCrMSysW5dxvGtL5xX+hyvUDa0pOakWHk/yG8ryvhnmoru9BOTqCUbmycoHcHZ3FXkmd64DKka8ITvaUqgN+peaNSs73WdGV4mAar1HJ1ylUp1L8G22/nODWsUuPX7PWlCIPlGKgVpaX6qjUc4DucN5ATxVcFJ82bVqceeaZbaO/b7jhhrjvvvvi5ptvjgsvvDDvc+rq6qKpqSnvY0mSxPTp0+Oiiy6KT3ziExER8cMf/jAGDx4c9957b5x88smFhliVCk3axegV2pUT4iWr1rTb0vPp1Isx6iWtInhar1EKxTgRKPQ1u6NU08ZlVTEKSKUoQuXrUddxpzxTskKN6c6ojnLodBrHPDrk69fm59zts9eePY6rK9K6gC606JlPOaY5LcWxlOoCu9DP0tSxneu0E1uHWZ06z9UavKG61HrjbTXN9FIM5ZiSvjvnbuU43yvGc7rTntXhPLy+4JfILCPW6YlinHOW4/q6/Xd46606f05a7ZCVmq+6U4zu7uum8RpZUY5jqaa/V8elx9rXyUqjs981ebZzBRXF33333Zg7d25Mnjy5bVt9fX2MGTMmZs+evdHnrVy5MnbaaadYv359HHjggfGNb3wj9tlnn4iIWLBgQbS2tsaYMWPa9u/fv380NzfH7Nmz8xbF16xZE2vW/O0f2fLlxV/HuhppwEpflpNtMRr3S6EYPftKEUdWdSfOd57JLWStXtoxNZRritZy8FtYOvJz9cjX0aUrF/edycpvaaU2BnRFVvN5JUnj30e5poHNSsO7XFw+cnNH3fn3V6n/ZvOdw7efYr0z5VpWrJqK01nNteUqGGTlcyBXWp1vya/W83Ol5tWILg4Kaad9Pu6zdeHvm9Xf1qzGVQxZbV+utSVEinEu8k4Rlu/tWASnUhX06b/55puxbt26GDx4cM72wYMHx/z58/M+Z4899oibb7459ttvv1i2bFlcffXVccghh8RTTz0VH/jAB6K1tbXtNdq/5obH2ps6dWpcdtllhYROF5TiBDnvj0UN/X5kpZCcVVnpIV4M1XRCUk0FcMqr2vJzMUZ1d7ZPVqY67Ir2y53suNUOJX/PWpsNpL1Kyi2FzjCTldHmxXjNrHwGaSnFKJ5KahjNumrLzVnVlX/D7fdZnuRO3dCdWdZKIau/aVmNq1Sy3KG+UFmJA7Kk1vJzVq97O7xm/44jPTu0bWekXbscS5H5/e45f8PSaz/Yqyse+1P7GdM66k4HmM5kpYNarV9/F94lokAHH3xwHHzwwW33DznkkNhrr73ie9/7Xlx++eXdes3JkyfHpEmT2u4vX748dtih9A2zpZJWT/WsfAm7o5p7oVUyiZ5yS2vKmFo/ediYasvPQH6VVKxn0yq5t7tc3DVyc3XpzhSt7WVl3c5CXzOt14By61puLs8UreVQq/lcfs6mcp0bty/ilWvpsXKopCVN6RmfQXmUo1NRreXigoriAwcOjF69esXixYtzti9evHija4a3t/nmm8cBBxwQzz//fERE2/MWL14cQ4YMyXnN/fffP+9rNDQ0REND+wXxqkd3/qFXcoG7GPwIUw38O05HrSX+Uqr2/FzJKrnw1Z6RwNB17Wd4SEtWRwfVgmrPzWYqyKZKzpvOCciiJauyUQSXz4un2vMzhSlFkbw7HZnLNegsjVwqf6fD37385NnOFVQU7927dxx00EHx0EMPxQknnBAREevXr4+HHnooJk6c2KXXWLduXfzxj3+MY489NiIihg0bFk1NTfHQQw+1FcGXL18ec+bMic997nOFhEcNqeYedUC6NIzWJieNpMEFIuTnN5k0pVX0ad+Brd+AHocBACWXL+dpR6lMrk+BiOpvGy94+vRJkybFhAkTYtSoUTF69OiYPn16rFq1Kk4//fSIiDj11FNj++23j6lTp0ZExJQpU+JDH/pQ7LrrrrF06dK46qqr4s9//nN89rOfjYiIurq6OPfcc+NrX/ta7LbbbjFs2LC4+OKLY+jQoW2FdwAqzPLXOt2lfW/3vluWKpie0TBfnczKAlDdNNCyQVnWIO2G9ufC226RO2KwO+sYtr62tkcxAaSt2hviqX75ZmracStT40O1ycosLhSu4KL4SSedFG+88UZccskl0draGvvvv3/MnDkzBg8eHBERL730UtTX17ft/5e//CXOPPPMaG1tjW222SYOOuig+O1vfxt777132z4XXHBBrFq1Ks4666xYunRpHHbYYTFz5szo06dPEQ6RWmDkOFS+rDZYQjXorOEdyL7257ulUIwLe7mXUvLvq+dcO0P3lSMXA+nJdy7s2pm0GcEPxVVwUTwiYuLEiRudLn3WrFk597/97W/Ht7/97U2+Xl1dXUyZMiWmTJnSnXAqTs1fyLcfQdpvSM7d9j3q9KaD4irGhXxW1ijNKr3bIf9vjYZ3NnBhD7RndD9Qk7owyxpkSVfatQtt+675tnKgy3RQK79qu07rVlGcwkjsAABAIVzsF6YYy2JU8oU9ANmkQzkAVCCd1qqWojgA9JDOT9A1pmwF0iRfVx4jzYA06aCWvmobnVaJ5NbyM+saQOkoilNU1iwFAACA0rP0GAAAQNcpigNAiendnj692zvRflqofkPSiQNqlJFoAABdY/kTAAphKZPSq6TcrCgOAFQdRfDKYDp1gNolV3ci3zqGOq0BQLp0KAeoaIriRebCHgAAKDe933uuknq3Q7HooAYAQFaZVY1iUxQvkKJ3ZejKj6WLfQAA0uDCnlrkWhrIMrm5esg3AMDGKIqTeflGvey41Q4pRAJQPkarAUC68jWqy8cAkC2K4F2nnQGAYugs92Y5NyuKA1S5YvR4NyVr8RV6cpDlk4lycyEP0E351iimIHJQ1zl3AQDo2KZmsBdAehTFO+FCHgCgIxf2AK4XN/B3qB7WGAcAAKqVojgAPWfkGSnSEA8AAKTNDGtUGtfSQNYUY8bTUliyak3aIVAkiuIAlIWTB8iO9t/HbbdoyN0hX0eXfkNKGBFQKA3vAABsiqJ39TCTC0BxKIpTs5xMAEC2yM1Uq6z2dgeyo9MOaxEdO62167BmaROofDqTAwCUjqI4QBXR6A4AAETkn1FCoRwAAKhViuLwf4xOAwAAAACgSzqZxaVUtGPDxhW81Fi+JQSpWorilFS+aZ/yTgNXINPCAQAAQGnlm4lKwzsAAFCJFMVJX0o96gAAAACoPZYeA7qi/YAvg70AKpuiOAAAAABQtTJTBDdFKwBAahTFAQDIJFO2AkD2WMcUAIBCOYckCxTFAcjRfhonAKBCGY0GAABVT7GRSpSZWVyoKYriAAAARVQpF/ft10gEALpOh3IAgMqiKA5Q41zIA8XQ/rdkx612SCkSAAAAAIBciuIAFaxSRqIBAD2jExsAAABA9ymKAwAAdJMOakCts44pWSQ/AwB5LX8t7QhIkaI4QAVxYQ+kpv1FQ78h6cQBAAApMGsLAEBlUxQHAAAAAAAos/Ydbnbcaocev2a+QTVmcqHcyjG4S4c1CqUoDhthCjjYBNPMAABAusziAgAA0GWK4mRPNy7sS9Gjrj097EiD6dIBIFvkZgDIlkzn5k46lC9ZtaZMgQAAXSU/Vy9FcQCAGtf+ZH/bLRp6/Jr5prAq1TRw76fDGjXLLC5AF6TRoVxuBqBm5DsnN5MLQGYoigPUmHKstaI3HQAAdF++8+lidFoDAGqTTmsAiuIAAACZU4pObDqtAQAAALVKUZyyK8UUrQAAEXq/U3yZXqcUAGqQ3AwAQHcoilOVyrFOGgAAAFQy186wceVYegwAgPJRFAcAoCw0vFMJjD4DAADyMQMqbJxraSqBojjZt/y1jtv6DSl/HFBiThwAoDYZiQYAAJSTpceoRK6d6SlFcQAAAAAgk3QgB4DapAhOsSmKQw/oUUdPuLAHKlr7mVzM4gIAdJFraQAAyiLfTMTULEVxADrXyclD+zWVAAAAInQGBoBC5BsZu+NWO6QQCUD1URSHItLbHQC6rv3Ffiku9PM1xMvPAFQls7gAAABslKI4AAAAAJAJZhcAqkqFdFoz2AuoBYriAADkyLckwrZbNKQQCQAAlF6+6YoLZs1SAIBMUxQHAKCm6AEPAPmVY2kTAAAqm1ldqFSK4lBCGt3JgqL0eAeAKlTRF/JdGI2Wb9YHAACAzuS7VtK2Xbsq+tq5E66ba4uiOAAAAEDGtW+ws7QJ1aCaG9kBisVMLgDFoSgOQC7roAEAAEDRGY0G1UWHNYDKoihO6spx8pBv+mg96ig3PeABIF1ZycWWNgHIZekxAABcK1NqiuIAAPRc+1km+g1JJw4AAACg5HRqAypNfXeedN1118XOO+8cffr0iebm5njsscc2uu+NN94Yhx9+eGyzzTaxzTbbxJgxYzrsf9ppp0VdXV3O7ZhjjulOaAAAAAAApbX8tdwbAACZVvBI8TvuuCMmTZoUN9xwQzQ3N8f06dNj7Nix8eyzz8agQYM67D9r1qwYP358HHLIIdGnT5+48sor46Mf/Wg89dRTsf3227ftd8wxx8Qtt9zSdr+hwfobbILRaFSArEzRClAp2k+TVa6lTjr7vdbbHQDKy8gzSs30rACF09ZZvZx7USsKLopPmzYtzjzzzDj99NMjIuKGG26I++67L26++ea48MILO+z/H//xHzn3f/CDH8Rdd90VDz30UJx66qlt2xsaGqKpqanQcAAAAABoL9/I1QI7lHelcFiuTmwAUBUM9oLSMnsLm1BQUfzdd9+NuXPnxuTJk9u21dfXx5gxY2L27Nldeo3Vq1fHe++9FwMGDMjZPmvWrBg0aFBss8028ZGPfCS+9rWvxbbbbltIeJB5+XrT6XVFJVqyak3aIQBAp4xkAIBskZsBIPvka6pVQUXxN998M9atWxeDBw/O2T548OCYP79rX5KvfOUrMXTo0BgzZkzbtmOOOSZOPPHEGDZsWLzwwgvxL//yLzFu3LiYPXt29OrVq8NrrFmzJtas+VtBZvny5YUcBkCPZXVKGVPAkSb5GUhbVi/ci5Kf9XanG+Tm6ta+o+q2W1iGDqASyM8AkI6Cp0/viSuuuCJuv/32mDVrVvTp06dt+8knn9z2/yNGjIj99tsvdtlll5g1a1YcffTRHV5n6tSpcdlll5UlZsqvXBf2aa1bClCt5GeA9JjFhXzkZgDIHvkZANJRX8jOAwcOjF69esXixYtzti9evLjT9cCvvvrquOKKK+LBBx+M/fbbb5P7Dh8+PAYOHBjPP/983scnT54cy5Yta7stWmRkJJXrnWfm59wAKpX8DADZIjdTLdpfN7t2BiqZ/Fy9lqxa0+FWCi+tWJRzywq5mrJb/lrHWyfK8R0luwoaKd67d+846KCD4qGHHooTTjghIiLWr18fDz30UEycOHGjz/vmN78ZX//61+OBBx6IUaNGdfo+L7/8cixZsiSGDBmS9/GGhoZoaDAtGJAdaZ3omY6VLJGfq5spWskajSzQObmZapbVJa1qnfwMnZOfKTYzolIIuZpaVvD06ZMmTYoJEybEqFGjYvTo0TF9+vRYtWpVnH766RERceqpp8b2228fU6dOjYiIK6+8Mi655JK47bbbYuedd47W1taIiNhyyy1jyy23jJUrV8Zll10Wn/zkJ6OpqSleeOGFuOCCC2LXXXeNsWPHFvFQqWrti3r98neoqAQu7LPJyQJA7cqXA+RnNiVLozUAclTRtTNUAiPQAACyo+Ci+EknnRRvvPFGXHLJJdHa2hr7779/zJw5MwYPHhwRES+99FLU1/9tVvbrr78+3n333fjUpz6V8zqXXnppfPWrX41evXrFH/7wh7j11ltj6dKlMXTo0PjoRz8al19+uR5zlJUedQAAZIaZXAAAgAqmgzmbojM5aSi4KB4RMXHixI1Olz5r1qyc+wsXLtzka/Xt2zceeOCB7oQBAEBW5SvoVdFotM5mEHGhDwBQZXRYAwCoaN0qigPlYzp1AGpVvl7DZnKBvzIdKwB0ndFoQGZZ2oQSq/VlQV07836K4lBhTDtTfP6mAJXDcie1o9Yv3AEAgMrXviC37RaWjKW4DKqDrlMUB8hDQ3wuPeoAAIBKooG4thkZDkCt0q4NG6coDgBAwfR2B4Bs6VZuNmUrAABQIxTFybx8I1TL0fBeSdOz6gEPAABAGirp2pniMhINgEJpx65dZnEhCxTFASqIkweAytCVRmIX/7k0rAPQE/IIabPsGACUWPtZjqBAiuJUJ1PAwcY5eQDSIj8DACnJVzTXQa0wNdfxwLUzkIZ8vz0FXjvnG1RTKTO5GEkOlJKiOAAAQBeZtQWAWlFzRXAAyq6zXKNTG1BMiuJQhfSoA4Dsq+V8rZG9HSPRAAAA6IRraegZRXHookqedqbWdHZyUEtFBwAqRy0XyQEAAIDqYZY1skhRHGqARvZcldSjzskDwKa1/53UYQ0A6CrXyrkq6VoZAKqRXNxzS1atSTsEMkxRHGqQtViySQEcAGqE6dIBKkI1F82rqdE9rWtpje4A6ajm/AyUlqI4tSFfw2O/IeWPAyqAC3ugO9r/dmy7RUPnT2qfn+XmqlFNDe0AFKaSZ3HpLH/pYA4AdJfr5G7QoZwiUxSnInWr4Z2C6HEHAOmSiwGArnLeAFAZdCgvPjkwG8yCSiVQFAcqipOcbtCjDoAS0+MdIHvyzQClQzlshOtmoMpV8kwuhdJ+DGyMojgQEYU3Zps2DoCiK0Lv92q+0O9Kri5GLu6sAaHWCuB6uwOwKaVoeNeYv2lyMwAQYRlQCqcoDnRJdxrANaoDQOWTrwFqSFdGyxbYaS1fAbOaOq1llfxdfBregWpSzR3K2+vK4K5ydHID0qcoDj1QSycPWeXkAgA2risX9nJpZdAQDwAAkF2urSH7FMWBopH4AQAA4G+KsfyJ6dQBqlf7zq/bbtGQUiS1xVKiUJsUxakK3Tp5KMK6pRRG0bwMujLdIQBVq7NcKxdXDiPDASpTMZYeK8V7AABd15VcKx+XgLZtSkxRHCAl+dbWS4NGd6CaWccUALLHUmQAAJUrK+3aUChFcdjAyHEAKJp8HW5MA0fN0tsdqCaunakGcjNADh3WgFqgKA5F5OQBAIooX2NlERre5Ws2RY93oJpZt5RKJDcDtKODGkC3KIoDlIkLeQAAgI50WAMAAEpNUZyqlJXe7y7sKTlTvgEAAFBFStKh3LUzUCG6tRRZCUaO5/st1rZduwz2olooigMAUBZZ6bQGG2Tlwj5fwxcAAAAAxaMoDlBjNLwD5DKzC0VnNBoAdElWOqgBAH9Tlvzchetm7dgUm6I4AAAAAJmhwxoAAFBsiuKwMdZiAQAAoIpZ2gQA6Cqd1ugRM6qRAYriACWSlWngTDMDVJUSdFqjdmQlNwNArSpbLi5Cw7traaBi5fsNLMG1syI5UGkUxQEqiR51QBUxOg0AqoAOawAAFU0HcmqFojg1IV/v3oIb3svUww4AgAzTQQ2oYlntsGYpMkrNqHAAKDLXzmSQojikzDQzAJAtcjOlpuEdoOfkawAogJldABTFAbqjktZBA6gURqdRbKaAA0iBRnfeRy4GgGyRm6lliuIAAFQuDe+8j4t7AKhBOpMDVS6rHcgBKo2iODUrqycTnTXmGplWxYpwIW86VgBqhVlbAEov3/VFVq6doSdcOwOVLKvt2tQ4bdtUAEVx6IkURqeZwhUA0mcdUwrhwh6oafkaSEtw7Sw3sylyMUA7JWjXlouBrFMUh/+jhx3vV5bRZ0aeAQAAUKFSW7bEtTQAdJllxuBvFMUBKpje7kA169aUrV1pJLXueEVyIQ8AANBRUQZ7lWlGVKPJS8+1M2ycojhQ85woAAAAdE2lzLKm0b34jAwHoBDaXKuY9cOpUIriUEwprDEOAKSvs4t9DfGFq5QGFBfyAN2QwrVzvrwiPwNAunRiqxCK4FQJRXEoJdPOAAAbIX9ngJFnAD1WSVO2AgBdlO9aSX7OpErpUA5ZoCgO1JyynSjoQQdQdBreq1dZ8rMiOAB0menSASBbFMChZxTFYSO6UozM6tpp7dX6SDQnCwBsUhmK5HIRAGSP5U9ypXK+UqICuA7mQC0rWbt2Rq6dqzk/Z6btwGAvqpSiOJRTmaad6Sx51vrJBQBkUS11YsvMhX43uLAHoBpUci4GoOeKMgsb1cMsLdQIRXFImylceyQzF/JOHABSUZIL+Yzm5nw5L6uF88zkZwDKriij07pzfZWRfE3p6aAGkBEpXTvXUodyoLgUxYG8TCNXfi7sAXquLEXyiLJc7CssF4FOawCZVCmj06pplrXMnFfIzQAUWWdF8krqYF4Wli6hhimKQw+U5UK+K0kqhR7x3TmZ6EovvsxcqHfGuioAFSHfb201jybvTKkaA6ppDdL25GuACtFZXkgpV5drNFtmrqVLkJ/lYoDSK0o7d0auk7uSE0uRn7tTnK9UcjOVqltF8euuuy6uuuqqaG1tjZEjR8Z3vvOdGD169Eb3v/POO+Piiy+OhQsXxm677RZXXnllHHvssW2PJ0kSl156adx4442xdOnSOPTQQ+P666+P3XbbrTvhQWo6SwZdOpmo4F7ThZ5MZPpEwIU8ABWgHDO7ZCZfy80AVSu1DudFaKwvRZ7MTO4FoGaVrIN5hepObs5MPje4C9oUXBS/4447YtKkSXHDDTdEc3NzTJ8+PcaOHRvPPvtsDBo0qMP+v/3tb2P8+PExderU+NjHPha33XZbnHDCCfH73/8+9t1334iI+OY3vxnXXntt3HrrrTFs2LC4+OKLY+zYsfH0009Hnz59en6UUO0y0guvvYpJ/CX6ezlZAMimonRia69CRqd1d59qITcDkLaKnsJVBzWAmlbwtXR38kaGr51TkdIAOvmZalVwUXzatGlx5plnxumnnx4RETfccEPcd999cfPNN8eFF17YYf9rrrkmjjnmmPjyl78cERGXX355/OIXv4jvfve7ccMNN0SSJDF9+vS46KKL4hOf+ERERPzwhz+MwYMHx7333hsnn3xyT44PMqUryaRka592JiOF9KIo9Pj1lgPgfcq2Lnl71ZSLO2PqcwA2IbVr5/a5uEwd0Ku5IV4uBqgd1TQFezWTm6llBRXF33333Zg7d25Mnjy5bVt9fX2MGTMmZs+enfc5s2fPjkmTJuVsGzt2bNx7770REbFgwYJobW2NMWPGtD3ev3//aG5ujtmzZ+ctiq9ZsybWrPnbF3fZsmUREbF8+fJCDqdLVr/zTtFfEzalO//mBjQWoTHg7QWF7b/V4J6/54rFnb9uvn3K4K3VTg7IvmLnvQ2vlyRJt55frvwsN1NuXfk3l0oujihOPi6FEuRvuZlKIDdD+RTj312H/N1ZLs73eKG5uCvXwd3R/nVLcG3dPhfnO/+Rr8maUrQXy8/QNUVp5+7OdXJ7aV03l+m6uLP2CLmZLErr2rmgovibb74Z69ati8GDc39EBg8eHPPnz8/7nNbW1rz7t7a2tj2+YdvG9mlv6tSpcdlll3XYvsMOFTLdFAAUw8WXluRlV6xYEf379y/4efIzADVPbgaAbClRbo6QnwGg21K6di54+vQsmDx5cs7o8/Xr18dbb70V2267bdTV1RXtfZYvXx477LBDLFq0KPr161e0102TY6oMjqkyOKbKUY3HVYpjSpIkVqxYEUOHDu3W88uRn32WlcExVY5qPC7HVBkcU9fIzempxuNyTJXBMVUGx1QZSnVM8nM6HFNlcEyVoRqPKaI6j8sxdU1Xc3NBRfGBAwdGr169YvHi3GkfFi9eHE1NTXmf09TUtMn9N/x38eLFMWTIkJx99t9//7yv2dDQEA0NuVNCbL311oUcSkH69etXNf/YNnBMlcExVQbHVDmq8biKfUzd6eW+QTnzs8+yMjimylGNx+WYKoNj6pzcnK5qPC7HVBkcU2VwTJWhFMckP6fHMVUGx1QZqvGYIqrzuBxT57qSm+sLecHevXvHQQcdFA899FDbtvXr18dDDz0UBx98cN7nHHzwwTn7R0T84he/aNt/2LBh0dTUlLPP8uXLY86cORt9TQAAAAAAAADoioKnT580aVJMmDAhRo0aFaNHj47p06fHqlWr4vTTT4+IiFNPPTW23377mDp1akREfOlLX4ojjzwyvvWtb8Vxxx0Xt99+ezz++OPx/e9/PyIi6urq4txzz42vfe1rsdtuu8WwYcPi4osvjqFDh8YJJ5xQvCMFAAAAAAAAoOYUXBQ/6aST4o033ohLLrkkWltbY//994+ZM2fG4MGDIyLipZdeivr6vw1AP+SQQ+K2226Liy66KP7lX/4ldtttt7j33ntj3333bdvnggsuiFWrVsVZZ50VS5cujcMOOyxmzpwZffr0KcIhdl9DQ0NceumlHaazqWSOqTI4psrgmCpHNR5XNR5TV1TjcTumylCNxxRRncflmCqDY6oe1Xrc1XhcjqkyOKbK4JgqQzUeU1dV47E7psrgmCpDNR5TRHUel2MqrrokSZKyvysAAAAAAAAAlEFBa4oDAAAAAAAAQCVRFAcAAAAAAACgaimKAwAAAAAAAFC1FMUBAAAAAAAAqFo1XRT/+te/Hoccckg0NjbG1ltv3aXnJEkSl1xySQwZMiT69u0bY8aMieeeey5nn7feeiv+3//7f9GvX7/Yeuut44wzzoiVK1eW4Ag6KvS9Fy5cGHV1dXlvd955Z9t++R6//fbby3FI3fp7HnXUUR3iPeecc3L2eemll+K4446LxsbGGDRoUHz5y1+OtWvXlvJQ2hR6TG+99VZ84QtfiD322CP69u0bO+64Y3zxi1+MZcuW5exX7s/puuuui5133jn69OkTzc3N8dhjj21y/zvvvDP23HPP6NOnT4wYMSLuv//+nMe78v0qtUKO6cYbb4zDDz88ttlmm9hmm21izJgxHfY/7bTTOnwmxxxzTKkPI0chxzRjxowO8fbp0ydnn0r7nPL9HtTV1cVxxx3Xtk/an9MjjzwSxx9/fAwdOjTq6uri3nvv7fQ5s2bNigMPPDAaGhpi1113jRkzZnTYp9DvaBZUY27uzvvLz/Jzd8nNcnMlfE5yc2Xl5ojqzM9y81/JzXJzd8nP8rP8nC65WW6Wm3umGvOz3Cw3y81dkNSwSy65JJk2bVoyadKkpH///l16zhVXXJH0798/uffee5N58+YlH//4x5Nhw4Ylb7/9dts+xxxzTDJy5Mjk0UcfTf73f/832XXXXZPx48eX6ChyFfrea9euTV577bWc22WXXZZsueWWyYoVK9r2i4jklltuydnv/cdcSt35ex555JHJmWeemRPvsmXL2h5fu3Ztsu+++yZjxoxJnnjiieT+++9PBg4cmEyePLnUh5MkSeHH9Mc//jE58cQTk//6r/9Knn/++eShhx5Kdtttt+STn/xkzn7l/Jxuv/32pHfv3snNN9+cPPXUU8mZZ56ZbL311snixYvz7v+b3/wm6dWrV/LNb34zefrpp5OLLroo2XzzzZM//vGPbft05ftVSoUe0z/8wz8k1113XfLEE08kzzzzTHLaaacl/fv3T15++eW2fSZMmJAcc8wxOZ/JW2+9VZbjSZLCj+mWW25J+vXrlxNva2trzj6V9jktWbIk53iefPLJpFevXsktt9zStk/an9P999+f/Ou//mty9913JxGR3HPPPZvc/8UXX0waGxuTSZMmJU8//XTyne98J+nVq1cyc+bMtn0K/TtlRTXm5u68v/wsP3eH3Cw3V8rnJDdXVm5OkurMz3LzX8nNcnN3yM/ys/ycPrlZbpabu68a87PcLDfLzV1T00XxDW655ZYunTysX78+aWpqSq666qq2bUuXLk0aGhqS//zP/0ySJEmefvrpJCKS3/3ud237/PznP0/q6uqSV155peixv1+x3nv//fdP/umf/ilnW1f+MZdCd4/pyCOPTL70pS9t9PH7778/qa+vz/lRvP7665N+/fola9asKUrsG1Osz+nHP/5x0rt37+S9995r21bOz2n06NHJ5z//+bb769atS4YOHZpMnTo17/6f+cxnkuOOOy5nW3Nzc3L22WcnSdK171epFXpM7a1duzbZaqutkltvvbVt24QJE5JPfOITxQ61ywo9ps5+D6vhc/r2t7+dbLXVVsnKlSvbtqX9Ob1fV77HF1xwQbLPPvvkbDvppJOSsWPHtt3v6d8pbdWSm4v5/vKz/NwZubkjubkyPie5uTJyc5JUT36Wm/9Gbi6taszNSSI/J4n8nAXy81/JzbnkZrm5K6oxP8vNcnMWVEJurunp0wu1YMGCaG1tjTFjxrRt69+/fzQ3N8fs2bMjImL27Nmx9dZbx6hRo9r2GTNmTNTX18ecOXNKGl8x3nvu3LnR0tISZ5xxRofHPv/5z8fAgQNj9OjRcfPNN8df/42XVk+O6T/+4z9i4MCBse+++8bkyZNj9erVOa87YsSIGDx4cNu2sWPHxvLly+Opp54q/oG8T7H+jSxbtiz69esXm222Wc72cnxO7777bsydOzfnu1BfXx9jxoxp+y60N3v27Jz9I/76N9+wf1e+X6XUnWNqb/Xq1fHee+/FgAEDcrbPmjUrBg0aFHvssUd87nOfiyVLlhQ19o3p7jGtXLkydtppp9hhhx3iE5/4RM53oho+p5tuuilOPvnk2GKLLXK2p/U5dUdn36di/J0qRdZzc7HeX36WnzsjN+cnN1fG5yQ3V1dujsh+fpabc8nNcnMh5Oe/kZ/l50oiN8vNxVDpuTmiOvOz3Pw3crPc3JnNOt+FDVpbWyMichLOhvsbHmttbY1BgwblPL7ZZpvFgAED2vYpZXw9fe+bbrop9tprrzjkkENytk+ZMiU+8pGPRGNjYzz44IPx//1//1+sXLkyvvjFLxYt/ny6e0z/8A//EDvttFMMHTo0/vCHP8RXvvKVePbZZ+Puu+9ue918n+OGx0qpGJ/Tm2++GZdffnmcddZZOdvL9Tm9+eabsW7durx/w/nz5+d9zsb+5u//7mzYtrF9Sqk7x9TeV77ylRg6dGjOD/YxxxwTJ554YgwbNixeeOGF+Jd/+ZcYN25czJ49O3r16lXUY2ivO8e0xx57xM033xz77bdfLFu2LK6++uo45JBD4qmnnooPfOADFf85PfbYY/Hkk0/GTTfdlLM9zc+pOzb2fVq+fHm8/fbb8Ze//KXH/54rRdZzc7HeX36Wn7vy3nJzR3Jz9j8nubn6cnNE9vOz3Pw3crPcXCj5+a/kZ/m50sjNcnMxVHpu3vD+1Zaf5ea/kpvl5q6ouqL4hRdeGFdeeeUm93nmmWdizz33LFNEPdfVY+qpt99+O2677ba4+OKLOzz2/m0HHHBArFq1Kq666qpuJ6VSH9P7k+qIESNiyJAhcfTRR8cLL7wQu+yyS7dfd1PK9TktX748jjvuuNh7773jq1/9as5jxf6c6Lorrrgibr/99pg1a1b06dOnbfvJJ5/c9v8jRoyI/fbbL3bZZZeYNWtWHH300WmEukkHH3xwHHzwwW33DznkkNhrr73ie9/7Xlx++eUpRlYcN910U4wYMSJGjx6ds73SPqdKU425OUJ+7g75WX4uJ7m5MsjN6anG/Cw3F05ulpvLTX6uDPJzOuTm7pObe0Zurm1yc2WQm4uj6ori5513Xpx22mmb3Gf48OHdeu2mpqaIiFi8eHEMGTKkbfvixYtj//33b9vn9ddfz3ne2rVr46233mp7fqG6ekw9fe+f/OQnsXr16jj11FM73be5uTkuv/zyWLNmTTQ0NHS6f3vlOqb3xxsR8fzzz8cuu+wSTU1N8dhjj+Xss3jx4oiITH9OK1asiGOOOSa22mqruOeee2LzzTff5P49/Zw2ZuDAgdGrV6+2v9kGixcv3ugxNDU1bXL/rny/Sqk7x7TB1VdfHVdccUX88pe/jP3222+T+w4fPjwGDhwYzz//fMmTUk+OaYPNN988DjjggHj++ecjorI/p1WrVsXtt98eU6ZM6fR9yvk5dcfGvk/9+vWLvn37Rq9evXr82RdTNebmCPlZfv6bLORnuTmX3FwZn5PcnF5ujqjO/Cw3y80byM2lIz/nJz/Lz8UgN8vNG4s3Qm4ulmrMz3JzfnKz3JxXj1clrwK33HJL0r9//073W79+fdLU1JRcffXVbduWLVuWNDQ0JP/5n/+ZJEmSPP3000lEJI8//njbPg888EBSV1eXvPLKK0WP/f16+t5HHnlk8slPfrJL7/W1r30t2Wabbboda1cV6+/561//OomIZN68eUmSJMn999+f1NfXJ4sXL27b53vf+17Sr1+/5J133ineAeTR3WNatmxZ8qEPfSg58sgjk1WrVnXpvUr5OY0ePTqZOHFi2/1169Yl22+/fTJ16tS8+3/mM59JPvaxj+VsO/jgg5Ozzz47SZKufb9KrdBjSpIkufLKK5N+/fols2fP7tJ7LFq0KKmrq0t++tOf9jjerujOMb3f2rVrkz322CP553/+5yRJKvdzSpK//tY3NDQkb775ZqfvUe7P6f0iIrnnnns2uc8FF1yQ7Lvvvjnbxo8fn4wdO7btfk8/+7RVS24uxvvLz/JzV8nNfyU3V8bnlCRyc6Xl5iSpnvwsN2+c3Fxc1Zibk0R+zkd+lp/TIjf/ldwsNxeiGvOz3NyR3Cw3542xx69Qwf785z8nTzzxRHLZZZclW265ZfLEE08kTzzxRLJixYq2ffbYY4/k7rvvbrt/xRVXJFtvvXXy05/+NPnDH/6QfOITn0iGDRuWvP322237HHPMMckBBxyQzJkzJ/n1r3+d7Lbbbsn48ePLckydvffLL7+c7LHHHsmcOXNynvfcc88ldXV1yc9//vMOr/lf//VfyY033pj88Y9/TJ577rnk3/7t35LGxsbkkksuKfnxJEnhx/T8888nU6ZMSR5//PFkwYIFyU9/+tNk+PDhyRFHHNH2nLVr1yb77rtv8tGPfjRpaWlJZs6cmWy33XbJ5MmTM3lMy5YtS5qbm5MRI0Ykzz//fPLaa6+13dauXZskSfk/p9tvvz1paGhIZsyYkTz99NPJWWedlWy99dZJa2trkiRJcsoppyQXXnhh2/6/+c1vks022yy5+uqrk2eeeSa59NJLk8033zz54x//2LZPV75fpVToMV1xxRVJ7969k5/85Cc5n8mG35AVK1Yk559/fjJ79uxkwYIFyS9/+cvkwAMPTHbbbbeSn6R295guu+yy5IEHHkheeOGFZO7cucnJJ5+c9OnTJ3nqqadyjruSPqcNDjvssOSkk07qsD0Ln9OKFSvaclBEJNOmTUueeOKJ5M9//nOSJEly4YUXJqecckrb/i+++GLS2NiYfPnLX06eeeaZ5Lrrrkt69eqVzJw5s22fzv5OWVWNubkr7y8//5X83DNys9xcKZ/TBnJzZeTmJKnO/Cw3y81yc/mOS36Wn4tFfv4buflv5Ga5uVDVmJ/lZrlZbu6ami6KT5gwIYmIDreHH364bZ+ISG655Za2++vXr08uvvjiZPDgwUlDQ0Ny9NFHJ88++2zO6y5ZsiQZP358suWWWyb9+vVLTj/99JwTklLq7L0XLFjQ4RiTJEkmT56c7LDDDsm6des6vObPf/7zZP/990+23HLLZIsttkhGjhyZ3HDDDXn3LYVCj+mll15KjjjiiGTAgAFJQ0NDsuuuuyZf/vKXk2XLluW87sKFC5Nx48Ylffv2TQYOHJicd955yXvvvZfJY3r44Yfz/luNiGTBggVJkqTzOX3nO99Jdtxxx6R3797J6NGjk0cffbTtsSOPPDKZMGFCzv4//vGPk9133z3p3bt3ss8++yT33XdfzuNd+X6VWiHHtNNOO+X9TC699NIkSZJk9erVyUc/+tFku+22SzbffPNkp512Ss4888yyX1gVckznnntu276DBw9Ojj322OT3v/99zutV2ueUJEkyf/78JCKSBx98sMNrZeFz2th3fMNxTJgwITnyyCM7PGf//fdPevfunQwfPjwnV22wqb9TVlVjbu7K+8vPfyM/94zcLDdXwueUJHJzJeXmJKnO/Cw3y81yc8/Iz/Kz/Jwuuflv5Ga5uTuqMT/LzXKz3Ny5uiRJkgAAAAAAAACAKlSfdgAAAAAAAAAAUCqK4gAAAAAAAABULUVxAAAAAAAAAKqWojgAAAAAAAAAVUtRHAAAAAAAAICqpSgOAAAAAAAAQNVSFAcAAAAAAACgaimKAwAAAAAAAFC1FMUBAAAAAAAAqFqK4gAAAAAAAABULUVxAAAAAAAAAKqWojgAAAAAAAAAVev/B5RZP7edBsQnAAAAAElFTkSuQmCC", "text/plain": [ "<Figure size 2000x500 with 4 Axes>" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axs = plt.subplots(1, n_replicas, sharey=True, figsize=(20, 5))\n", "for i in range(n_replicas):\n", " ax = axs[i]\n", " ax.set_title(f\"replica {i+1}\")\n", " for j in range(n_temps):\n", " ax.hist(result[i][j][:, :, 0].flatten(), density=True, alpha=0.245, label=\"$\\\\beta_\"+str(j)+\"$\", bins=100)\n", "plt.tight_layout()\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "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.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }