From aae63ac65d955c3651aa42d51d36c2b603ce4e67 Mon Sep 17 00:00:00 2001 From: Dipa C Date: Tue, 23 Dec 2025 02:12:38 -0500 Subject: [PATCH] Sampling A1 --- .../a1_sampling_and_reproducibility.ipynb | 72 +++++++++++++++---- 1 file changed, 60 insertions(+), 12 deletions(-) diff --git a/02_activities/assignments/a1_sampling_and_reproducibility.ipynb b/02_activities/assignments/a1_sampling_and_reproducibility.ipynb index 264a448a..e2f23768 100644 --- a/02_activities/assignments/a1_sampling_and_reproducibility.ipynb +++ b/02_activities/assignments/a1_sampling_and_reproducibility.ipynb @@ -8,43 +8,78 @@ "# Assignment 1: Sampling and Reproducibility\n", "\n", "The code at the end of this file explores contact tracing data about an outbreak of the flu, and demonstrates the dangers of incomplete and non-random samples. This assignment is modified from [Contact tracing can give a biased sample of COVID-19 cases](https://andrewwhitby.com/2020/11/24/contact-tracing-biased/) by Andrew Whitby.\n", - "\n", - "Examine the code below. Identify all stages at which sampling is occurring in the model. Describe in words the sampling procedure, referencing the functions used, sample size, sampling frame, any underlying distributions involved. \n" + " \n", + "> Examine the code below. Identify all stages at which sampling is occurring in the model. Describe in words the sampling procedure, referencing the functions used, sample size, sampling frame, any underlying distributions involved." ] }, { "cell_type": "markdown", "id": "4ea73db3", "metadata": {}, - "source": [] + "source": [ + "**Code Block 1 - Infect random subset of people:** \n", + "Function --> \n", + "   infected_indices = np.random.choice(ppl.index, size=int(len(ppl) * ATTACK_RATE), replace=False) \n", + "   ppl.loc[infected_indices, 'infected'] = True \n", + "Sample size = 100 (1000 * 10%) \n", + "Sampling frame = 1000 \n", + "Sampling type = simple random sampling (without replacement) \n", + "Distribution = binomial \n", + " \n", + "**Code Block 2 - Primary contact tracing:** \n", + "Function --> \n", + "   ppl.loc[ppl['infected'], 'traced'] = np.random.rand(sum(ppl['infected'])) < TRACE_SUCCESS \n", + "Sample size = 20 (100 * 20%) \n", + "Sampling frame = 100 \n", + "Sampling type = simple random sampling (with replacement, although doesn't really matter for this) \n", + "Distribution = binomial \n", + "\n", + " \n", + "**Code Block 3 - Secondary contact tracing:** \n", + "Function --> \n", + "   event_trace_counts = ppl[ppl['traced'] == True]['event'].value_counts() \n", + "   events_traced = event_trace_counts[event_trace_counts >= SECONDARY_TRACE_THRESHOLD].index \n", + "   ppl.loc[ppl['event'].isin(events_traced) & ppl['infected'], 'traced'] = True \n", + "Sample size = anywhere from nil (0) to all primary contact traced (20) \n", + "Sampling frame = 20 \n", + "Sampling type = convenience \n", + "Distribution = conditional\n", + " \n" + ] }, { "cell_type": "markdown", "id": "3d9b2ccc", "metadata": {}, "source": [ - "Modify the number of repetitions in the simulation to 10 and 100 (from the original 1000). Run the script multiple times and observe the outputted graphs. Comment on the reproducibility of the results." + "> Modify the number of repetitions in the simulation to 10 and 100 (from the original 1000). Run the script multiple times and observe the outputted graphs. Comment on the reproducibility of the results." ] }, { "cell_type": "markdown", "id": "4cf5d993", "metadata": {}, - "source": [] + "source": [ + "The 1000 repetition is quite overlapping between the infections vs traced, and over every re-run of the simulation. This is expected as the 1000 repetitions if basically averaging the results over the population many times (actually, as many times as the sample population of this simulation). \n", + "The 10 repetitions look very noisy, both between the re-runs and in overlap between the infections vs traced. Again, to be expected because it's a very small sample of the random simulation, thus any random individual variations can greatly skew the plot at each run. \n", + "Interestingly (for me), the 100 repetitions is not as noisy as I expected. Yes it is noiser than the 1000 repetition, but definitely better than the 10 repetitions...almost like the accuracy is logarithmic (probably not the best word) - in the sense that as you increase your repetitions/n, your noise decreases exponentially. I am not too familiar with stats, but I believe the law of large numbers shows something similar, where as your number of repetitions/n increases, the standard error/noise decreases exponentially (sorry I'm not explaining this well)." + ] }, { "cell_type": "markdown", "id": "32603ce7", "metadata": {}, "source": [ - "Alter the code so that it is reproducible. Describe the changes you made to the code and how they affected the reproducibility of the script. The output needs to produce the same output when run multiple times." + "> Alter the code so that it is reproducible. Describe the changes you made to the code and how they affected the reproducibility of the script. The output needs to produce the same output when run multiple times." ] }, { "cell_type": "markdown", "id": "77613cc3", "metadata": {}, - "source": [] + "source": [ + "I set a random seed (np.random.seed(50)) before the simulation so that all random generators (e.g. np.random.choice & np.random.rand) would have the same start point at each repetition, and carry out their defined \"randomness\" from there (which is usually systematic leaps). This allows the simulation to run the same sequence and produce the same outcomes every time, removing variability that would come with true randomness (and in turn, improving reproducibility)." + ] }, { "cell_type": "markdown", @@ -56,10 +91,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 30, "id": "ab8587a0", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAJOCAYAAACqS2TfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYrFJREFUeJzt3QeYVOXdN+CHIiAIoqggFjQWxK6YKNEYO5b4ajSJSSxojDHEXmL01dhjwUSNPTEqGk2MGjWJHWsSu9hxgz1gEBQbZaXPd/2fvLPf7LrALuxh2d37vq5hmJkzZ545Z87s/M7T2pVKpVICAAAAmlz7pl8lAAAAEIRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG6Agv3+979P66yzTlpiiSVSz549m7s4LKB27dql008/vbmLQR0HHnhgWm211Vr8Z6I530d9Zs2alU444YS0yiqrpPbt26c999xzkb32u+++m7ft8OHDF9lrAhRJ6AYWWPwgih9Gzz33XGrp7rnnnkIC1b/+9a/8Y3qNNdZIV199dfrtb3873+e8+OKLab/99ss/djt37pyWXXbZtMMOO6TrrrsuzZ49OxXliSeeyNvg008/TUX6wx/+kC6++OL5Lhdlic/X/C7bbLNNoeVtKx599NFa2zVOEn3pS19KBxxwQHr77bebu3gtVjlA1nfZYostCt+ft9122wI9/9prr00XXHBB+ta3vpWuv/76dMwxxzTbd8Gi9re//S19/etfTyussELq2rVrPg6+853vpPvuu6+5iwa0UB2buwAAi0vovvzyy5s8eMcP3zlz5qRf//rXac0115zv8r/73e/Sj3/849S7d++0//77p7XWWitNnjw5PfTQQ+nggw9O77//fvrf//3fVFToPuOMM/JJgiJr5OOH9quvvpqOPvroeS6311571dpmU6ZMSUOHDk3f/OY382Nlsa0Whc8//zx17Nj6/2weeeSR6ctf/nKaOXNmev755/OJorvvvju98sorqW/fvmlxEyez4hhb3H3ve99Lu+66a637ll9++bS4evjhh9NKK62ULrrookX+XdCvX798vMWJn0Xtl7/8ZfrpT3+aQ/dJJ52UQ/ebb76ZHnzwwXTzzTennXfeeZGXCWj5Wv+vB4Bm9MEHH+TrhoTYp556KgfuQYMG5ZMA3bt3r3ksfpRGi4L4gdpWbLjhhvlSNnHixBy6475oCTA306ZNS506dcpNYptSly5dUlvwta99LdduhoMOOiitvfbaOYhHbWeEkIUxderU1K1bt9SUmiOYLYhNN910np/bxfG7q7m6w0QNfXMcb9Gk/qyzzko77rhjeuCBB+b6fb4oFXHMAIue5uVAk4pa0qWWWiqNGTMmfeMb38j/j9qSqEUOUVu23Xbb5R8RUZsRNR31NVn/+9//ng499NDUq1ev1KNHj9zE9ZNPPqm17F/+8pe022675dq3aIYdTbjjB1N9TbCffvrpXMu0zDLL5NeO4Ba1z+Uyl8tX2fRzfq644oq03nrr5deOMhx22GG1mmZH/8zTTjutpkZrfv0/o5Y5lrnppptqBe6yzTbbLJe18sfYcccdV9MMvX///rmWplQq1XperPPwww9Pd955Z1p//fXzslHuyqaSUa6o3Qmrr756zTaIprEhmrbHfovmlvH8ddddN1155ZX1vo9777031xLFe4h9F7Wm5f0cTcGj1vTf//53zWssTD/WchPaqIE65ZRT8mctaqYmTZqUPv7443T88cenDTbYIH8Ooyy77LJLeumll+oN6rENImDGj/0VV1wx16a/9dZbtbZj5f4rN3+PWrBy64Cll146B9Xq6upa649auwiuyy23XN4u//M//5P+85//NLhPcPzYj5YOUasf5dtoo41yCK6vGXN8BqJ2Oo6H2Fex/Z999tkF3MIp7/fwzjvv1NrHEc7jWIr3E8fhqFGj6v0uiG0Yx14st+++++bHyq0/Yt/E+4njI2oQ63ZVufHGG9PAgQPTkksumbtZfPe7301jx479wuuUP0NROx/LxT6oKz4T8VrxmSibPn16PkajRUVsqziWoh9z3F8pbkfz6ihnef+99957qUjlz3ZcN1V/54Z8Zsvrf+SRR/I+LR+n5XLEvosm4fEdEtszPpPxXV33+3lhvgvm9h6j9r38uYuy77HHHqmqqqrR73Fu4sRefE623HLLeh+P77/GHpeN2Y+L6ph544030t5775369OmT17Xyyivn5T777LN5bh9gwanpBppchN4IN1tvvXUaNmxYDpER+uKH0sknn5x/RESgueqqq3KYjprdCHqVYvn4sRQ/oEaPHp0DXvw4K/+ACfFjJX6gHHvssfk6fpCdeuqp+UdT9EUsGzFiRD4BEEHqqKOOyj804ofaXXfdlW/HD8Zx48bl5WLQs4aIckVIjr7WUftaLmOEm8cffzzXvsUP0xtuuCHdcccd+bEoY2XNbaX4MRhNyGObrbrqqvN9/QjW8cM/fhjHj76NN9443X///Tk4R5ir2yT0n//8Z7r99tvTT37yk/xD7pJLLsk/uuLkSJzYiP3x+uuvpz/+8Y/5uREOK5u/RvnjR3a8ZjSxjj6Psa74IRgnG8pin/zgBz/Iy0ataOzDF154IQf873//+3n/xw+7CCzlMsZ2WVhxsiVqtyNQRUCK/7/22mv5RMO3v/3t/PmaMGFC+s1vfpNDQDxWbiodn9f4fMT2jx+e8ZmIJv3xeYiWBRFe5yX6esb6zz333NwcO7oIxI/z888/v9aP6VtuuSV3GYh+vI899lgOqg0RgT0CSoSIOC7itW699da8zjjJE+WtFKEmyh+f6zhW4hiM/Rv9shekVrh84iE+JyGOkSFDhqTBgwfn9xif3fh8bLXVVnlfV55EiZrDWC4ei5MBcUIkxGc2PivxPfHDH/4wL/ePf/wjt/aIk0vhF7/4Rfr5z3+et28s8+GHH6ZLL700HyPxOvXVwsb7i+4H8VmPfR2fg7L4LMRnI/ZxiM9ufJ7j2PjRj36UBgwYkE8KxucyjoVYvixeP8JMfIa/+tWv5u+ahu6/sthOEeoqRRhsjpr6eX1m45iPfRzbP7p0xDIhtk+Iz1XsuwixcSIpTsZcdtlleZ+Uv/uK+C6I5t3xeYn+1fH9G8dFfB4iIMd7qHvyriHHZV3xeITV+H474ogjcmhtquOyoYo+ZmbMmJHXH8dCvMf4exh/M+LvYZQ7PpNAAUoAC+i6666LKtXSs88+W3PfkCFD8n3nnHNOzX2ffPJJackllyy1a9eudPPNN9fc/69//Ssve9ppp31hnQMHDizNmDGj5v5hw4bl+//yl7/U3FddXf2FMh166KGlrl27lqZNm5Zvz5o1q7T66quX+vXrl8tRac6cOTX/P+yww/L6G+KDDz4oderUqbTTTjuVZs+eXXP/ZZddltdx7bXX1twX7y3u+/DDD+e5zpdeeikvd9RRRzWoDHfeeWde/uyzz651/7e+9a28nd98882a+2K5KG/lfeXXu/TSS2vuu+CCC/J977zzzhder75tPXjw4NKXvvSlmtuffvppqXv37qXNN9+89Pnnn891W++22255fzRWbMO6n5dHHnkk3xflqFvG+AxU7p8Q761z586lM888s+a+2F+xjgsvvPALr1lZ7rqvXd63P/jBD2o955vf/GapV69eNbdHjhyZlzv66KNrLXfggQd+YZ31ufjii/NyN954Y819cWwMGjSotNRSS5UmTZpU895iuXjtjz/+uGbZOGbi/r/97W/zfJ3ytoztEdt63Lhxpbvvvru02mqr5c9UHOeTJ08u9ezZs3TIIYfUeu748eNLSy+9dK37y98FJ554Yq1lH3744Xz/kUceOdft/e6775Y6dOhQ+sUvflHr8VdeeaXUsWPHWvfH61R+nu6///563++uu+5a6/P6+9//vtS+ffvSP/7xj1rLXXXVVfn5jz/+eL794osv5ts/+clPai33/e9/v0H7r7xf6rvENp/b+yjvj8plKtcX35XzUn7+rbfe2ujPbPj6179eWm+99WrdF9sqnn/TTTfVuv++++6rdf/CfhfU9x433njj0gorrFD66KOPan2PxT484IADFug91ufUU0/Nz+/WrVtpl112yZ+1OIYX9LhszH5cFMfMCy+88IXPBVA8zcuBQsRZ9rI4ux5Nn6OmO87Al8V98Vh9IyNHzVNlDVDUJkcNa/R1LosaibKo2YtapGh6GDVKMWp4iLP7URMTfaLr1ow1pAn53GpcorYg1lnZb/iQQw7JTSijyWRjRe18qK9ZeX1iO3To0CHXNFWK5uaRD6NZZ6Woka+ssY0a9yhrQ0elrtzWUTsV2zpqjOP55SaJUTMc++HEE0/8Qn/MBd3WDRU1r5VlDNFcuLx/ojb7o48+yjVp8bmLmq+yP//5z7lmP2p96mpIuaMffqX4DMZrlfdpuRl/tAyoVN/rzW1fR21UDMRVFsdG7PuoiYxa80r77LNP7kZRWZ7Q0H0dtZNR2xktAaI2N7oxRJPZqE2LfRy1YVGW+AyUL/FZ3HzzzXPLi7ri2K0U2zu2a7nrRX3bO2qqoyY6vi8qXye2QwwuWN/rVDaHj/35pz/9qea+aPocZY9tUxa1klF7G9P5Vb5GuTl9+TXK3zl1j7X5DQRY33dalKHyEs2Rm8P8PrNzE9ssakKjz3PlNovmzHFslbdZU38XxACSMatD1CJX1j7H91iUpfLvwsK+x2jBFK1FNtlkk9x6KGrk4/1Fn/zKpuyNPS4bo8hjplyTHe9tfs3tgaajeTnQ5Mr9zSrFH/roN1b3B1fcX19fwPiRUCl+0EXz8HIf4xD9DaMfbzT1rPtDqhwEy01joy9zU4lm7iHCW6VoyhpNH8uPN0YE4BA/VBtahghFdUN6uQlo3TLU12Q9gll9274+0Ww0fvA9+eSTX/ihFts69mMR27qh6nZPqOwDGX3v48RLZV//clPpEOWOfbmgI5PX3bblwBvbNvZr7IsI/3XL2JDR7EM8P46HugPDNXRfV5anIaKLRgSUCNIRXuN1ytsm+oKGcjCd2+e4LJ4Xx32l2N7x2Z1X0914nTh5VPd7oGxeTbLjNaPrRASnaEIbJ18ikER/78rQHa8RIWpuI4iXB80q77+63QzqHv/zE+8lTn4tDub3mZ2b2GZxvNft21x3mzX1d8HcvnNDfD4jQNYdcGxB32OIIB2X+LsS44FEs+74PO2+++65y0n8jWvscdlQRR8z8T0UXbIuvPDC3PUrjvXoZhGD/GlaDsURuoEmFz/WG3N/3YG/GiJq26KmNX48nXnmmfkHcfwQihrMn/3sZy1iCqG6ASx+bEWf0iIszLaPH3zbb799rhGMH2ox2FScYIianuiLuThs67q13OGcc87J/Ruj5jb6fMcP1viBHDWUTVnmpvxcLw7liYGa5hYOy9st+vxGDVpddU9cVLY2aIx4nThBFy026ns/8xsHIPptR5/ueP6ee+6Z+9PH57eyZjleI95rfKbrE5/z5jS3GuH6BopcVJ+R2GYRuCOs1WdxmgKtKY7L+PsSNelxidAaLT4ihMffnqL246I4Zn71q1/lVgMxGGmM0h6189H3PfqH1w38QNMQuoHFUpy133bbbWtuR3O9aGJYnuc2BlSLpoJRgxWDxJRVjrAcyrVTUTsxr1qmxjR5jFHXQwyeFjXbZdHkPF5/QWqzYrCcqD2MWvsYaXZ+P/ijDNHMPWrGK2u7y83qy2VsjLltgxhUKGoM//rXv9aqParbxLdyW8+rFrfopuZlt912W/4MXXPNNV84YVMeKK5c7vghHTWhRQxqFfsifhDHZ6OyFioGYGro819++eW8jsof4wuzrxdUeR9H8FrQWttYR9RMxujyc6u5i2UiHEWtXIwo31jxnRAtY6KJeQxIFcdVNBOu+xoxkn2cUJrXZ7K8/8otIsri+C9SuWa2ckaEhalBbQqxzeJ7JwYvq+9EV+VyTfldUPmdW1ccB3E8Fz2tVnSviNAdf4cac1w2xX4s4piJE05xidZiTzzxRN6nMbjp2Wef3eByAQ2nTzewWIopjyIElcXoyDFaa4zcGspn8itrLSL0RlPiStEPL36ExEjidX/0VD63/IOt7jL1ibARNb0xAnjlOiLcRdPLxo5qXBbNt2N9McJ1nGSoa+TIkTXT0cTJh6gpiVGDK0XNc/yQLW+nxpjbNqhvW8f7jGnEKu200075BEDUmMQUXPPa1otiapood91areiTGiP1VoqmyNH3se62bKra6hgpONT9bMaowg0R+3r8+PG1+ijHsRDPj9qrxtS6NcV7idq/aEVQeXyWxWjJ8xPbO7Zr9J2d2/aO0dZj/8UydfdB3I4TbvMSISjmGo8TRlErH9ursml5iL6v8Vm4+uqr6x2ZOporh/KxFMd7pfhOKVKEttgGMX1ipbqfo0Uptll870TLkbpiG5e/O5r6uyBOoMQMDfH9V/n9FKE+amrLJ2MXVnSdiS409SmPk1E+8dLQ47Ip9mNTHjPRZD7KWSnCdxwzdafKA5qOmm5gsRQBOmqg4kde1G7ED5SosYq+ZyGm7YkahBhAK5rGRdCMH9d1f2zED4kI7NEXL360xTQ38QMuaiOiT3jUHoQYKCfEuiJYxI+X8tRC9TWhjClw4sdNzJMaZSqXMeahjb5xCyLeU8wXHgNuRVPYCN9RMxq12VGzHzXN5VqIeD9Rixu1d9HPPZrNxo/PaC4YzafnN81VfcrbINYZ7z1qfeN14gd0nGSI/8d0QXFCIIJK1HaWa31ChLEI/TGIXmyHmBYo9lHUJsaP2fIJg3id+KEa/QpjufiBGutuajENWHQ9iH0e2zaa7kez2MrWCSGmrYup3aI8zzzzTO7jGIEravRiX8RcwAsj3m/8aI6QFj98y1OGxbRUDantiwG4oql0NAeNEy8xNVLU4kc/+1hnQwffawqxj+N4is9mnNCKz0kcDzH1XAwgGLVl9Z28qBSf23h+hNho0RLHUNQWxvRH8VhMvxSf3/isx3EWn+9oIh7vM1oLxBR8sU0q59uuT4TsCEBxMitCRbmvbVmUIZqdx4Bb0Wojyh6BMr4b4v74bojazfjeiP69cXxHQIzPUkwv19CWCgsq+tfGdHfxHuIzEtskpnUq95tuDhEk4zsgwnQMbBbfDfE9EfsxTmjFGApxsqOI74KYBjJOgMQUkzF9VnnKsNhODZnrviGibLF/4xiNz2W0OIqQH9PHxeczPocxwFpjjsum2I9NecxEq49YPsoUNeIRwONvZ/zNi+8poCCLYIR0oI1NGRZTrdRV3/QzIaaLiWlj6q7zscceK/3oRz8qLbPMMnn6lX333bfWVDEhpvTZYost8nRkffv2LZ1wwgk10wXVnZ7ln//8Z2nHHXfM09hE+TbccMNa02XF1GJHHHFEafnll8/TIzXk6zGmCFtnnXVKSyyxRKl3796loUOHfmFasoZOGVYppqeJ6YjiPcW6Yxtsv/32peuvv77WFFgxfdMxxxxTs9xaa62Vp/2qnJInxOvHlGj1bfvYX5XOOuus0korrZSn4amcPuyvf/1r3mZdunTJU0idf/75NVNt1Z1iLJb96le/mvdLjx49Sl/5yldKf/zjH2senzJlSn5/MfVUPL+h04fNa8qw+qa/iSnDjjvuuNKKK66Yy7LllluWnnzyyfxZjEulmG7s5JNPztPLxbbs06dPnn7trbfemu+UYXX3bfkzXLldpk6dmvfBsssumz/Pe+65Z2n06NF5ufPOO2++733ChAmlgw46qLTccsvl6d822GCDL0wZVZ6CKD4DdTVkaqt5bcv6lo0p42KasPhMrLHGGnkKtOeee26+3wXl4y3KGcdPvJ847mJ6prpTM/35z38ubbXVVnk9cYnlYzvGtqt8nfo+Q3EcrLLKKvVOrVc5xVN8luO7KaaSi2Mtpis844wzSp999lnNcjHtVUzXFFNORTl233330tixYxs1ZVh9+6VSfe8jPlt77713ngYxyhZTIr766qsLPWVYQz6zc/vODr/97W/zdorjKr5T4/MY378xzVxTfBfMbVq0Bx98MB/H5fXFfnjttddqLdOY91jXzJkzS1dffXU+PqMs8ZmIbb/JJpvk/Td9+vRGH5eN2Y+L4ph5++2383RqcczGsRvfSdtuu23etkBx2sU/RQV6gMaKUWKjZvLZZ5/NtUzQWkVNYdSa3XjjjWnfffdt7uIAAAXRpxsAChZNYeuKJqjR/aFyIEAAoPXRpxsACjZs2LDc7zP6X8a0WjEoU1yin2VzT00FABRL6AaAgsXgTCNGjMijPsdAdDH1Wgz+VHcaKwCg9dGnGwAAAAqiTzcAAAAUROgGAACAgujTnVKaM2dOGjduXOrevXtq165dcxcHAACAxVj00p48eXLq27dvno1kXoTulHLgNnosAAAAjTF27Ni08sorz3MZoTulXMNd3mA9evRo7uIAAACwGJs0aVKuuC1nyXkRumMI9/9rUh6BW+gGAACgIRrSPdlAagAAAFAQoRsAAAAKInQDAABAQfTpBgCAxdTs2bPTzJkzm7sY0OYsscQSqUOHDk2yLqEbAAAWwzmAx48fnz799NPmLgq0WT179kx9+vRp0GBp8yJ0AwDAYqYcuFdYYYXUtWvXhf7RDzTupFd1dXX64IMP8u0VV1wxLQyhGwAAFrMm5eXA3atXr+YuDrRJSy65ZL6O4B3H4sI0NTeQGgAALEbKfbijhhtoPuVjcGHHVRC6AQBgMaRJObSOY1DoBgAAgILo0w0AAC3EmDFj0sSJExfZ6y233HJp1VVXLXzQuP333z898cQTeZqmokZsHz58eDr66KObdUT4GKDr0EMPTbfddlv65JNP0gsvvJA23njj1BYMb8D2P/3009Odd96ZXnzxxXz7wAMPzMvHfS2Z0A0AAC0kcPfvPyBNm1a9yF6zS5euafToqgYH7wUJSRdddFF6//33c9BaeumlU1NYbbXVcsCLS9k+++yTdt1119Sc7rvvvhw+H3300fSlL30pn9RoTieeeGLeV//6179q7ov/DxgwIA0ZMiSXtSz+HycMYv+WBxkr2q9//et8oqKlE7oBAKAFiBruCNwDBtyYunYdUPjrVVdXpaqq/fLrFlnb/dZbb6WBAwemtdZaKxUpguKiCovzeq8x/dRXv/rVuS4zY8aM1KlTp0VSnm233Tadf/75ubVBzEcdHnnkkbTKKqvkEwOV4v4ttthikW7DpZvoJExz06cbAABakAjc3btvWvilKYL9Nttsk4488sh0wgknpGWXXTYHu2hCXFkj/ec//zndcMMNedCqqCkPUZv6wx/+MC2//PKpR48eabvttksvvfRSrXX/7W9/S1/+8pdTly5dco3xN7/5zZrX/Pe//52OOeaYvM7yYFhRU9uzZ89a67jyyivTGmuskUNu//790+9///taj8dzf/e73+V1x0jWcWLgr3/9a83j0UR83333zeWMMBqPX3fddfVui3hvRxxxRG6xEOuN914u7+GHH55r5eN9DB48ON//2GOPpa985Supc+fOOahHrfSsWbNqbdtYXzxvmWWWSb17905XX311mjp1ajrooINS9+7d05prrpnuvffeue6frbbaKjfprwzY8f/DDjssffzxx+ndd9+tdX+E9DB9+vR0/PHHp5VWWil169Ytbb755l8I6bG942RN165d8/b76KOPvvD65513Xi53lPXggw9O06ZN+8I223PPPRv8eSrX1Mf7is/Fuuuumx588MG8vcutL+KkRmzv2KaxTL9+/dK5556biiR0AwAAhbn++utzMHv66afTsGHD0plnnplGjBiRH3v22WfTzjvvnL7zne/kJubRnDh8+9vfzvMjR2AcOXJk2nTTTdP222+fg2C4++67c5CL5uLRL/qhhx7KATXcfvvtaeWVV86vE+uMS33uuOOOdNRRR6Xjjjsuvfrqq7npdITVqNGtdMYZZ+Tyvfzyy/n1ImSXy/Hzn/88vfbaa7mcVVVVOcTPrcl4vLcoU5QtyhTvvXIbRfB//PHH01VXXZX+85//5NeKkwpxsiHWe80116Szzz77C9s2Xu+ZZ57JAXzo0KF520VN+vPPP5922mmn3F++urr+LgmxX+I1Kt9zhOfY1ltuuWXN/W+//XY+WVAO3RFan3zyyXTzzTfn7RKvGfvxjTfeyI/Hvo4QHcu9+OKL+Xl1y37LLbfkwHzOOeek5557LofgK664ot5yNvTzFHPcR0iPoB+P//a3v00nn3xyredfcskl+cRJvP7o0aPTTTfdVHMCpCialwMAAIXZcMMN02mnnZb/HzXBl112WQ7JO+64Y64hjprcqCUuN2/+5z//mUNkhO54LPzyl7/MNZUxANmPfvSj9Itf/CJ997vfzYG4bKONNsrXUQPaoUOHXHtaXmd9Yp1Rk/qTn/wk3z722GPTU089le8vh8sQy3zve9/L/4+AGKEtyhchM4LoJptskjbbbLP8+LzCWzSVjjJF2eqWK7ZLBMiyCIrRxDu2VdTSrrPOOmncuHHpZz/7WTr11FNT+/bta97zKaeckv9/0kkn5ZrjCOGHHHJIvi+WjcAewTiahtcn3uutt96a/x8nEKK2Od7T1ltvnQN4nIiI66gVjnXEe47a/Lju27dvfl7Uekd/9bg/tlGcYIjtEzXSYe21184D5cUyZRdffHEO5nEJEcqjVrpubXdjPk8RvqMJf5S3vI3jsxKPlUW543lRGx7bNmq6i6amGwAAKEyEpEpRoxmBem6iZnfKlCmpV69eaamllqq5vPPOOzlQhag9jdrYhRE101GbWylux/1zK3/UsEZz93L5o2Y5antjBPIImBEsF0T0aa9btkGDBtWaJzrKFtvlvffeq7dsEeZjm22wwQY190XT7TCv7R1Ntl9//fVc+x5hNcJorOvrX/96TZPxuI7a8zgJ8sorr+Qa5QjSlfsnmsOX90+UP5qcVxo0aNAX3uP8lmns5ylqruNkReVJjXILiMqTKPH5ie4E0VT9gQceSEVT0w0AABQm+gxXiiA5Z86cuS4fwTKCVN0+wqHcJ3tRDuY1r/Lvsssuuf/4Pffck2tZ40RA9IeO2vLGiDDfVGWrvK8c2ue1vSPMR9P2aEoelwjbIZqdxyB60bQ89kU0vy/vnwjl0ew/ritF+F7cPk91RVeFOIETXQKiZj26Duywww65FUVR1HQDAACLjQhFMZp2x44d80BglZdyf+mo7YwmxXMTITJqY+clpsWKPtSV4nYMvtUY0UQ+pte68cYbc5Pp6Ee8sKJs0We6crqsKFs0T48+4U0pTmCUB0KL2uqo+S6H22hOHn3Jx44dW9PkPpqex7aN2uW6+6dcwxzljz7VlZ566qkvvMf5LdNYUXsdZZ0wYULNfZV958uitUJMIRcDz/3pT3/Kg/mV++kXQU03AACw2Ihax2hmHANiRT/naMYc/ZnLg6dF/+no0xu1yjHyePTtjlG9o7Y5+jyX+1b//e9/z49Fk+j6Bjf76U9/mms5I0TGa8Zo6DEIW9R+NlT0mY6m4eutt14e0fuuu+7KYXJhRT/zCPAxOFoMRhbNpuM9R7/zcn/uphSBOuZLL5/0KIta76i1Lw+4FmJ/xGByBxxwQPrVr36Vt9+HH36YT4LEyZDddtstN9uOGvR47h577JHuv//+Wv25QwxiF029Y3/GsjGg2ahRo/L85Qsq+m7HZyJOgsRnZ/LkyTV93su1/hdeeGFuSRHljm0Z/dnjZEHdke2bktANALRqMWhONJGkdYtQVeRc0ouTmD+7Nb1OXRGOIkDHYGIxiFcEughFMbBXuY9y1MZGWDrrrLPy4GFRcxmPl8WI1tEcOgJYhOHKGuOyCPUx4FcEwwiAq6++eh4IrFzT2xBRox4DmMXUWlFj/LWvfS338V5YMRVXbIM4MRCDpcXgcDHgWDlAFhG6Y5vF4GfRwqAydEfYj2nMKpt1x3aKgc9i5PcYaT2Ov6gV/8Y3vpEfj/9HLXI899RTT80nNaLssb/KoqY5+oBHX/gYPG3vvffOfeQjoC+oaO4eA+7FdHNxkiAC/AUXXJB23333PBBciNYCEchjpPVYPpaLbV3EyYyydqX6PoFtzKRJk/Jogp999lk+YAGA1hO4+/cfkKZNq3+6HFqPLl26ptGjq1pF8I4AEn1OIwSWg0JzfZ5b03albXr88cfz4HBvvvlmPgnTFMdiYzOkmm4AoNWKGu4IKAMG3Ji6dl34Jp8snqJGtqpqv7y/W3M4jPcWAXhRttxoSy0IaB3uuOOOPKBbTAsWQTtaMUTz9cYG7qYkdAMArV4E7u7d/38/RWipIgALwTB30Y87+vZHy5A4aRRN26PveXMSugEAAGgVDjjggHxZnJgyDAAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABTElGEAANBCxNzDEydOXGSvF/Mct8Z5wbfZZpu08cYbp4svvji1FgceeGD69NNP05133tng973aaqulo48+Ol8ojtANAAAtJHAP6N8/VU+btshes2uXLqlq9OgGBe927drN8/HTTjstnX766akliHJGeH3xxRcXaj1bbLFFDrlXXXVVzX3x/6FDh6brrrsuB+Wy+P9bb72V/vGPf6RF5dlnn03dunVbZK/XVgndAADQAkQNdwTuGwcMSAO6di389aqqq9N+VVX5dRsSut9///2a///pT39Kp556aho9enTNfUsttVTN/0ulUpo9e3bq2LF1x5Ftt9023XHHHbXue+SRR9Iqq6ySHn300VqhO24PGTJkkZZv+eWXX6Sv11bp0w0AAC1IBO5Nu3cv/NLYYN+nT5+ay9JLL51rvsu3//Wvf6Xu3bune++9Nw0cODB17tw5/fOf/8w1u3vssUfq3bt3DuVf/vKX04MPPlhrvdOnT08/+9nPclCN56255prpmmuuqXn81VdfTbvsskt+fqxn//33r9UEf+rUqemAAw7Ij6+44orpV7/61Tzfx/Dhw9MZZ5yRXnrppfwe4hL3lVsbRHljXT169Ejf+c530oQJE+YZuuPEw/jx42vue+yxx9KJJ56YQ3bZO++8k/7973/n5cPYsWPzunv27JmWXXbZ/JrvvvtuzfJxwuLYY4/Nj/fq1SudcMIJ+URGpYa872heXtnEPt7r7373u/TNb34zde3aNa211lrpr3/9a63nxO24v0uXLrm8119/fX5eNG0P8T523333tMwyy+Ra9PXWWy/dc889qS0TugEAgEUiwuZ5552Xqqqq0oYbbpimTJmSdt111/TQQw+lF154Ie288845sEW4LYvg+Mc//jFdcskl+Xm/+c1vamrNI+htt912aZNNNknPPfdcuu+++3IIjsBa9tOf/jQH3b/85S/pgQceyGH3+eefn2sZ99lnn3TcccflsBi193GJ++bMmZPD78cff5zXN2LEiPT222/nx+Zmyy23TEsssUSu3Q6vvfZa+vzzz9PBBx+cPvrooxy2QzweIXbQoEFp5syZafDgwfkkRTQ1f/zxx/P7jW0zY8aMvHwE6DgRcO211+aTF1GmujXqjX3fZXHCIbbfyy+/nPfNvvvum9cforzf+ta30p577plPShx66KHp5JNPrvX8ww47LJ8o+fvf/55eeeWVdP7559dq5dAWte72HAAAwGLjzDPPTDvuuGPN7ajF3WijjWpun3XWWTk8Rm3q4Ycfnl5//fV0yy235IC7ww475GW+9KUv1Sx/2WWX5cB9zjnn1NwXQTRqxeO5ffv2zbXiN954Y9p+++3z41Ezu/LKK8+1jEsuuWQOidH0PWrpy6IMESIjeMb6ww033JDDefSNjlr6uqKm9ytf+UoOvN/73vfy9VZbbZVr7L/61a/m26uvvnq+jsAd90dZI+BHjXO5n3z0/45a7Vhup512yrXTJ510Utprr71q+onff//9Na8bJzMa+77Losl7lDXEdo2THc8880wO/XHCo3///umCCy7Ij8f/o6XBL37xi5rnxwmTvffeO22wwQZf2F9tlZpuAABgkdhss81q3Y5wePzxx6cBAwbkUBlhN2qzyzXdMZBZhw4d0te//vV61xe1rVFLHM8rX9ZZZ538WDRdj0vUDm+++ea1gn6ExcaKckXYLgfusO666+Zyx2PzGjG83JQ8ruN2iPdUeX+5aXm8pzfffDPXdJffU5R52rRp+f189tlnufa98j3FCYLKbbsw7ztaIFSeNIhm9B988EG+HU3l655ciJMKlY488sh09tln51r+GDzv5ZdfTm2d0A0AACwSdUfKjsAdNdtRoxpNqSNkRw1puRl11DrPS4T2aI4ez6u8vPHGG2nrrbdOi4MI01Hr/p///CeH6/IJhHLojoAcfbijmXz5PUW/97rvKdbx/e9/v/DyRnP4SlHbHjXvDfXDH/4wN7uPvvXRMmCzzTZLl156aWrLhG4AAKBZRH/laM4cA3dF2I7m3JUDhsV9Efiib3J9Nt100zRq1Kg8IFgMsFZ5iYC/xhpr5BD59NNP1zznk08+yQF2Xjp16pQHK6sUtfERjuNSFn20o1951HjPTTQjj/VdccUVubY6AnWIGuMPP/wwN4cvN0Mvv6c4abDCCit84T3FAHVxiYHRKt/TrFmz0siRI2tuL+j7np+oKY++85WiaX1d0Rrgxz/+cbr99ttz//irr746tWVCNwAA0CxiFOwIZlGTG82qoya3slY1wnRMo/WDH/wgz5sd/amjdjj6eZcH7YpBvqIPcoS/qDWOvs0HHXRQDs3RNDsGLYtBxR5++OHc/zhCfvv2845B8brxWlGuGAk9BgaLPuVxEiAGFosByaKfcwzyFjXWdZvNV4ra+pivO2p7o8l1NJcPEcQr7y/XMMf6l1tuuTxoW9T+l99zNNt+77338jJHHXVUHpAutkmMDP+Tn/ykZvTwsKDve35i4LR4vRhNvtzfvjyye7n/+dFHH533QZQ7ttMjjzyST1i0ZQZSAwCAFiTmz24tr3PhhRfmQB21wRE0I8xNmjSp1jJXXnll+t///d8cLGPE75gzPG6HGCgtasvjeTHAWITjfv365UG/ygEzBv0qN0OPftJR8xr9ouclBgKLkwHRNDzCbAxkFqE1RgI/4ogjctP1WH+8TkOaTsd6YjTvcn/usgjsEUrL/blDTNUVy8Z7ioHSJk+enFZaaaU8IFr0rw7xHqJfd5yQiHLENozWApXva0He9/zEoG+33XZbXtevf/3rPPhbjF4+dOjQPAhciJMdcTIkThBEeXfeeed00UUXpbasXanuhG5tUBzY0UwjPoTlDzIA0PJFLUs05Rw4cGTq3n3T5i4OBZk8+fk0cuTA3Lw2mua2dNEEOWoJI+DENFJlMbjYgP79U/W0aYusLF27dElVo0fnoAv1iZHLY/T0ymb3rcW0uRyLjc2QaroBAKAFiOAbATiaOy8qUfsscFMp+qZHf/RevXrlVgZRox7TuzF3QjcAALQQEYCFYJpTDPIWU4JFX/r4LEZT85gznLkTugEAAGiQ6J/d1vtoN5bRywEAAKAgQjcAAAC0xtB9+umn5/ncKi/rrLNOrdHiYrj56KQfc83F0P0TJkyotY4YxXG33XbLQ+vHBPIxF11MDg8AAC1Z5XzVQMs9Bpu9T/d6662XHnzwwZrbHTv+/yIdc8wx6e6770633nprHo49RsWLuepilLzyHHARuPv06ZOeeOKJPFddTFAfE8ufc845zfJ+AABgYXTq1CnPvTxu3Li0/PLL59tROQUsGjGr9owZM9KHH36Yj8U4Blt06I6QHaG5rpjv7Jprrkl/+MMf0nbbbZfvi0npBwwYkJ566qm0xRZbpAceeCC99tprObT37t07bbzxxumss87KE8lHLfrCbhwAAFjU4kd+zAscFUoRvIHmEa2pY4T2OCZbdOiOIef79u2bJxsfNGhQOvfcc/MbGzlyZJo5c2baYYcdapaNpufx2JNPPplDd1xvsMEGOXCXDR48OA0dOjSNGjUqbbLJJs30rgAAYMFF5VH87o1uk9G6E1i0OnTokCuIm6KVSbOG7s033zwNHz489e/fP5/JO+OMM9LXvva19Oqrr6bx48fnL5uePXvWek4E7HgsxHVl4C4/Xn5sbqZPn54vZZMmTWridwYAAAsnfuxHt8m4AC1Xs4buXXbZpeb/G264YQ7h/fr1S7fccktacsklC3vdqE2PgA8AAABtZsqwqNVee+2105tvvpn7eUfn9U8//bTWMjF6ebkPeFzXHc28fLu+fuJlJ510Uu4zXr6MHTu2kPcDAABA27ZYhe4pU6akt956K6244opp4MCBuSnNQw89VPP46NGj8xRh0fc7xPUrr7ySPvjgg5plRowYkXr06JHWXXfdub5O586d8zKVFwAAAGhVzcuPP/74tPvuu+cm5TEy42mnnZY7rH/ve9/LU4QdfPDB6dhjj03LLrtsDsZHHHFEDtoxiFrYaaedcrjef//907Bhw3I/7lNOOSXP7R3BGgAAANps6H7vvfdywP7oo4/yHIRbbbVVng4s/h8uuuiiPDz73nvvnQc+i5HJr7jiiprnR0C/66678mjlEca7deuWhgwZks4888xmfFcAAACwGITum2++eZ6PxzRil19+eb7MTdSS33PPPQWUDgAAAFpRn24AAABoTYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgHYtaMQBASzNt2pg0c+bE5i4GjVRdXZWvq6r+e13Xcsstl1ZdddVFXCqA/xK6AQD+L3A/90z/NGvOtOYuCgtov/32q/f+rl26pKrRowVvoFkI3QAAKeUa7gjcZ3cdkFZv37W5i0MjzJldnaZ+XpXWHTAgdetae99VVVen/aqq0sSJE4VuoG2H7vPOOy+ddNJJ6aijjkoXX3xxvm/atGnpuOOOSzfffHOaPn16Gjx4cLriiitS7969a543ZsyYNHTo0PTII4+kpZZaKg0ZMiSde+65qWPHxeatAQAtSATuAR27N3cxaITZKaXJKaWNu3ZN3bvbd8DiZbEYSO3ZZ59Nv/nNb9KGG25Y6/5jjjkm/e1vf0u33npreuyxx9K4cePSXnvtVfP47Nmz02677ZZmzJiRnnjiiXT99den4cOHp1NPPbUZ3gUAAAAsZqF7ypQpad99901XX311WmaZZWru/+yzz9I111yTLrzwwrTddtulgQMHpuuuuy6H66eeeiov88ADD6TXXnst3XjjjWnjjTdOu+yySzrrrLPS5ZdfnoM4AAAAtOnQfdhhh+Xa6h122KHW/SNHjkwzZ86sdf8666yT++I8+eST+XZcb7DBBrWam0cT9EmTJqVRo0YtwncBAAAAX9SsHZ+jr/bzzz+fm5fXNX78+NSpU6fUs2fPWvdHwI7HystUBu7y4+XH5ib6h8elLEI6AAAAtJqa7rFjx+ZB02666abUpUuXRfraMdDa0ksvXXNZZZVVFunrAwAA0DY0W+iO5uMffPBB2nTTTfNI43GJwdIuueSS/P+osY5+2Z9++mmt502YMCH16dMn/z+u43bdx8uPzU2Mkh59xsuXOAEAAAAArSZ0b7/99umVV15JL774Ys1ls802y4Oqlf+/xBJLpIceeqjmOaNHj85ThA0aNCjfjutYR4T3shEjRqQePXqkddddd66v3blz57xM5QUAAABaTZ/umENx/fXXr3Vft27dUq9evWruP/jgg9Oxxx6bll122RyMjzjiiBy0t9hii/z4TjvtlMP1/vvvn4YNG5b7cZ9yyil5cLYI1gAAANBmB1Kbn4suuii1b98+7b333nngsxiZ/Iorrqh5vEOHDumuu+5KQ4cOzWE8QvuQIUPSmWee2azlBgAAgMUudD/66KO1bscAazHndlzmpl+/fumee+5ZBKUDAACAFjZPNwAAALRWQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIB2LWjEAtFZjxoxJEydObO5i0ABVVVX5urr6v9fz0pBlAKCxhG4AaGTgHtC/f6qeNq25i0IjVFXt1+BlS6UZhZYFgLZF6AaARoga7gjcNw4YkAZ07drcxWE+plZXp9eqqlK3JQek9h3mvb8en/VRumLau6k0Z9YiKx8ArZ/QDQALIAL3pt27N3cxmI/JKaWI0N07dE0dOs57f70zu3qRlQuAtsNAagAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAABgcQrdb7/9dtOXBAAAAFqZBQrda665Ztp2223TjTfemKZNm9b0pQIAAIC2Grqff/75tOGGG6Zjjz029enTJx166KHpmWeeafrSAQAAQFsL3RtvvHH69a9/ncaNG5euvfba9P7776etttoqrb/++unCCy9MH374YdOXFAAAANrSQGodO3ZMe+21V7r11lvT+eefn9588810/PHHp1VWWSUdcMABOYwDAABAW7VQofu5555LP/nJT9KKK66Ya7gjcL/11ltpxIgRuRZ8jz32aLqSAgAAQAvTcUGeFAH7uuuuS6NHj0677rpruuGGG/J1+/b/zfCrr756Gj58eFpttdWaurwAAADQukP3lVdemX7wgx+kAw88MNdy12eFFVZI11xzzcKWDwAAANpW6H7jjTfmu0ynTp3SkCFDFmT1AAAA0Hb7dEfT8hg8ra647/rrr2+KcgEAAEDbDN3nnntuWm655eptUn7OOec0RbkAAACgbTYvHzNmTB4sra5+/frlxwCgJYi/WRMnTmzUc6qqqvL11OrqNLmgctF0Yj8BQIsL3VGj/fLLL39hdPKXXnop9erVq1EDssXl3XffzbfXW2+9dOqpp6Zddtkl3542bVo67rjj0s0335ymT5+eBg8enK644orUu3fvWj+Yhg4dmh555JG01FJL5X7kURMfc4gDwNzE34/+/QekadMWLJS9VlWVZjV5qSjKnFIpdWjuQgDQJi1QMv3e976XjjzyyNS9e/e09dZb5/see+yxdNRRR6Xvfve7DV7PyiuvnM4777y01lprpVKplPuDx9zeL7zwQg7gxxxzTLr77rtzX/Gll146HX744WmvvfZKjz/+eH7+7Nmz02677Zb69OmTnnjiifT++++nAw44IC2xxBKauQMwT1HDHYF7wIAbU9euAxr8vOrqqlRVtV/qtuSA1L1D10LLyMKbOevjNG3aO/l3BgC0mNB91lln5drp7bffvqZGec6cOTnwNibs7r777rVu/+IXv8g130899VQO5DHl2B/+8Ie03Xbb1QzgNmDAgPz4FltskR544IH02muvpQcffDDXfm+88ca5bD/72c/S6aefnkdQB4B5icDdvfumjX5e+w5dU4eO3QspE01n9hzNywFogQOpRZj905/+lP71r3+lm266Kd1+++3prbfeStdee+0CB92otY5m5FOnTk2DBg1KI0eOTDNnzkw77LBDzTLrrLNOWnXVVdOTTz6Zb8f1BhtsUKu5eTRBnzRpUho1atRcXyuaqscylRcAAABoagvV8XnttdfOl4Xxyiuv5JAd/bejT/Ydd9yR1l133fTiiy/mAN+zZ89ay0fAHj9+fP5/XFcG7vLj5cfmJvp8n3HGGQtVbgAAACgkdEet9PDhw9NDDz2UPvjgg9y0vNLDDz/c4HX1798/B+zPPvss3XbbbXkgtOgfXqSTTjopHXvssTW3o6Z7lVVWKfQ1AQAAaHsWKHTHgGkRumMQs/XXXz+1a9dugQsQtdlrrrlm/v/AgQPTs88+m37961+nffbZJ82YMSN9+umntWq7J0yYkAdOC3H9zDPP1FpfPF5+bG46d+6cLwAAALDYhe7oe33LLbekXXfdtckLFLXm0ec6AniMQh616XvvvXd+bPTo0XmKl2iOHuI6Bl+L2vaYxiyMGDEi9ejRIzdRBwAAgBYXuitrpxe2mXfMyR2Do02ePDmPVP7oo4+m+++/P08RdvDBB+dm4Msuu2wO0kcccUQO2jFyedhpp51yuN5///3TsGHDcj/uU045JR122GFqsgEAAGiZofu4447LTcAvu+yyhWpaHjXUMc1YzK8dIXvDDTfMgXvHHXfMj1900UWpffv2uaY7ar9jZPIrrrii5vkdOnRId911Vxo6dGgO4926dct9ws8888wFLhMAAAA0a+j+5z//mR555JF07733pvXWWy83A68UU4g1RMzDPS9dunRJl19+eb7MTb9+/dI999zTwJIDAADAYh66Y2Czb37zm01fGgAAAGjrofu6665r+pIAAABAK9N+QZ84a9as9OCDD6bf/OY3eRC0MG7cuDRlypSmLB8AAAC0rZruf//732nnnXfO03fFAGcx8Fn37t3T+eefn29fddVVTV9SAAAAaAs13UcddVTabLPN0ieffJKWXHLJmvujn3fMqw0AAAAsYE33P/7xj/TEE0/k+borrbbaauk///lPU5UNAAAA2l5N95w5c9Ls2bO/cP97772Xm5kDAAAACxi6d9ppp3TxxRfX3G7Xrl0eQO20005Lu+66a1OWDwAAANpW8/Jf/epXafDgwWnddddN06ZNS9///vfTG2+8kZZbbrn0xz/+selLCQAAAG0ldK+88srppZdeSjfffHN6+eWXcy33wQcfnPbdd99aA6sBAABAW9ZxgZ/YsWPab7/9mrY0AAAA0NZD9w033DDPxw844IAFLQ8AAAC07dAd83RXmjlzZqqurs5TiHXt2lXoBgAAgAUdvfyTTz6pdYk+3aNHj05bbbWVgdQAAABgYUJ3fdZaa6103nnnfaEWHAAAANqqJgvd5cHVxo0b15SrBAAAgLbVp/uvf/1rrdulUim9//776bLLLktbbrllU5UNAAAA2l7o3nPPPWvdbteuXVp++eXTdtttl371q181VdkAAACg7YXuOXPmNH1JAAAAoJVp0j7dAAAAwELWdB977LENXvbCCy9ckJcAAACAthm6X3jhhXyZOXNm6t+/f77v9ddfTx06dEibbrpprb7eAAAA0FYtUOjefffdU/fu3dP111+flllmmXzfJ598kg466KD0ta99LR133HFNXU4AAABoG326Y4Tyc889tyZwh/j/2WefbfRyAAAAWJjQPWnSpPThhx9+4f64b/LkyQuySgAAAGh1Fih0f/Ob38xNyW+//fb03nvv5cuf//zndPDBB6e99tqr6UsJAAAAbaVP91VXXZWOP/749P3vfz8PppZX1LFjDt0XXHBBU5cRAAAA2k7o7tq1a7riiitywH7rrbfyfWussUbq1q1bU5cPAAAA2lbz8rL3338/X9Zaa60cuEulUtOVDAAAANpi6P7oo4/S9ttvn9Zee+2066675uAdonm56cIAAABgIUL3Mccck5ZYYok0ZsyY3NS8bJ999kn33XffgqwSAAAAWp0F6tP9wAMPpPvvvz+tvPLKte6PZub//ve/m6psAAAA0PZquqdOnVqrhrvs448/Tp07d26KcgEAAEDbDN1f+9rX0g033FBzu127dmnOnDlp2LBhadttt23K8gEAAEDbal4e4ToGUnvuuefSjBkz0gknnJBGjRqVa7off/zxpi8lAAAAtJWa7vXXXz+9/vrraauttkp77LFHbm6+1157pRdeeCHP1w0AAAAsQE33zJkz084775yuuuqqdPLJJxdTKgAAAGiLNd0xVdjLL79cTGkAAACgrTcv32+//dI111zT9KUBAACAtj6Q2qxZs9K1116bHnzwwTRw4MDUrVu3Wo9feOGFTVU+AAAAaBuh++23306rrbZaevXVV9Omm26a74sB1SrF9GEAAABAI0P3Wmutld5///30yCOP5Nv77LNPuuSSS1Lv3r2LKh8AAAC0jT7dpVKp1u177703TxcGAAAANNFAanML4QAAAMAChu7or123z7Y+3AAAANAEfbqjZvvAAw9MnTt3zrenTZuWfvzjH39h9PLbb7+9MasFAACAVqlRoXvIkCFfmK8bAAAAaILQfd111zVmcQAAAGjTFmogNQAAAGDuhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAADQGkP3ueeem7785S+n7t27pxVWWCHtueeeafTo0bWWmTZtWjrssMNSr1690lJLLZX23nvvNGHChFrLjBkzJu22226pa9eueT0//elP06xZsxbxuwEAAIDFKHQ/9thjOVA/9dRTacSIEWnmzJlpp512SlOnTq1Z5phjjkl/+9vf0q233pqXHzduXNprr71qHp89e3YO3DNmzEhPPPFEuv7669Pw4cPTqaee2kzvCgAAAP6rY2pG9913X63bEZajpnrkyJFp6623Tp999lm65ppr0h/+8Ie03Xbb5WWuu+66NGDAgBzUt9hii/TAAw+k1157LT344IOpd+/eaeONN05nnXVW+tnPfpZOP/301KlTp2Z6dwAAALR1i1Wf7gjZYdlll83XEb6j9nuHHXaoWWadddZJq666anryySfz7bjeYIMNcuAuGzx4cJo0aVIaNWpUva8zffr0/HjlBQAAAFpt6J4zZ046+uij05ZbbpnWX3/9fN/48eNzTXXPnj1rLRsBOx4rL1MZuMuPlx+bW1/ypZdeuuayyiqrFPSuAAAAaMsWm9AdfbtfffXVdPPNNxf+WieddFKuVS9fxo4dW/hrAgAA0PY0a5/ussMPPzzddddd6e9//3taeeWVa+7v06dPHiDt008/rVXbHaOXx2PlZZ555pla6yuPbl5epq7OnTvnCwAAALTamu5SqZQD9x133JEefvjhtPrqq9d6fODAgWmJJZZIDz30UM19MaVYTBE2aNCgfDuuX3nllfTBBx/ULBMjoffo0SOtu+66i/DdAAAAwGJU0x1NymNk8r/85S95ru5yH+zoZ73kkkvm64MPPjgde+yxeXC1CNJHHHFEDtoxcnmIKcYiXO+///5p2LBheR2nnHJKXrfabAAAANps6L7yyivz9TbbbFPr/pgW7MADD8z/v+iii1L79u3T3nvvnUcdj5HJr7jiipplO3TokJumDx06NIfxbt26pSFDhqQzzzxzEb8bAAAAWIxCdzQvn58uXbqkyy+/PF/mpl+/fumee+5p4tIBAABAKxm9HAAAAFoboRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIB2LWjEAACxKU6ur53pfVVVVM5SIIiy33HJp1VVXbe5iQIMJ3QAAtGhzSjPmGqxf/7/r/fbbbxGXiqJ06dI1jR5dJXjTYgjdAAC0aKXSrHzdpUv/tETHpWo91m12dUqfV6UBA25MXbsOaKYS0lSqq6tSVdV+aeLEiUI3LYbQDQBAq9C+fdfUoWP32vc1W2koUrlVg6bmtARCNwAArdbEOTNy8I7aUVqPcneBrl26pKrRowVvFmtCNwAArdbk0qw0J6V0ZufV0hpL9Gru4rCQ5syuTlM/r0rrDhiQxkT4rqrS1JzFntANAECrt3r7JdOAOk3PaXlmx4mUlNLGXbumbs1dGGgg3VwAAACgNYbuv//972n33XdPffv2Te3atUt33nlnrcdLpVI69dRT04orrpiWXHLJtMMOO6Q33nij1jIff/xx2nfffVOPHj1Sz54908EHH5ymTJmyiN8JAAAALGahe+rUqWmjjTZKl19+eb2PDxs2LF1yySXpqquuSk8//XTq1q1bGjx4cJo2bVrNMhG4R40alUaMGJHuuuuuHOR/9KMfLcJ3AQAAAIthn+5ddtklX+oTtdwXX3xxOuWUU9Iee+yR77vhhhtS7969c434d7/73TxVwH333ZeeffbZtNlmm+VlLr300rTrrrumX/7yl7kGHQAAAJrLYtun+5133knjx4/PTcrLll566bT55punJ598Mt+O62hSXg7cIZZv3759rhkHAACA5rTYjl4egTtEzXaluF1+LK5XWGGFWo937NgxLbvssjXL1Gf69On5UjZp0qQmLj0AAAAsxjXdRTr33HNzrXn5ssoqqzR3kQAAAGiFFtvQ3adPn3w9YcKEWvfH7fJjcf3BBx/UenzWrFl5RPPyMvU56aST0meffVZzGTt2bCHvAQAAgLZtsQ3dq6++eg7ODz30UK1m4NFXe9CgQfl2XH/66adp5MiRNcs8/PDDac6cObnv99x07tw5TzFWeQEAAIBW1ac75tN+8803aw2e9uKLL+Y+2auuumo6+uij09lnn53WWmutHMJ//vOf5xHJ99xzz7z8gAED0s4775wOOeSQPK3YzJkz0+GHH55HNjdyOQAAAG06dD/33HNp2223rbl97LHH5ushQ4ak4cOHpxNOOCHP5R3zbkeN9lZbbZWnCOvSpUvNc2666aYctLfffvs8avnee++d5/YGAACANh26t9lmmzwf99y0a9cunXnmmfkyN1Er/oc//KGgEgIAAEAr7NMNAAAALZ3QDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACtKxqBXT9MaMGZMmTpzY3MWgAMstt1xaddVVm7sYAABAExO6W1DgHtC/f6qeNq25i0IBunbpkqpGjxa8AQCglRG6W4io4Y7AfeOAAWlA167NXRyaUFV1ddqvqirvY6EbAABaF6G7hYnAvWn37s1dDAAAABrAQGoAAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAgnQsasUAbd2YMWPSxIkTm7sYzEVVVVW+rq7+73VDNXZ5AKBtE7oBCgrcA/r3T9XTpjV3UZiPqqr9Fuh5pdKMJi8LAND6CN0ABYga7gjcNw4YkAZ07drcxaEeU6ur02tVVanbkgNS+w4N30ePz/ooXTHt3VSaM6vQ8gEArYPQDVCgCNybdu/e3MWgHpNTShGbu3fomjp0bPg+emd2daHlAgBaFwOpAQAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQToWtWKAlmzMmDFp4sSJC/z8qqqqfD21ujpNbsJy0XRi3wAAFE3oBqgncPfvPyBNm7bwoey1qqo0q0lKRVHmlEqpQ3MXAgBotYRugDqihjsC94ABN6auXQcs0Dqqq6tSVdV+qduSA1L3Dl2bvIwsvJmzPk7Tpr2TSqVScxcFAGjFhG6AuYjA3b37pgu1jvYduqYOHbs3WZloOrPnaF4OABTPQGoAAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAACgIEI3AAAAFEToBgAAgIII3QAAAFAQoRsAAAAKInQDAABAQYRuAAAAKIjQDQAAAAURugEAAKAgQjcAAAAUROgGAACAggjdAAAAUBChGwAAAAoidAMAAEBBhG4AAAAoiNANAAAABRG6AQAAoCBCNwAAABRE6AYAAICCdCxqxQAAAEWYWl2dpv7f/++5555UVVXVzCWiqX3pS19KgwYNSq2B0A0AALQIc0oz8nWE7Kf+r9nuz3/+8+YuFgVon1L65xNPtIrgLXQDAAAtQqk0K1936dI/zS5NT3Omv5tO79Q3rdFh6eYuGk3ordmfpdNnjEtvv/220L04ufzyy9MFF1yQxo8fnzbaaKN06aWXpq985SvNXSwAAKCJtW/fNbUr/Xd4qgjc63Xu3dxFoilNj3/GpdaiVQyk9qc//Skde+yx6bTTTkvPP/98Dt2DBw9OH3zwQXMXDQAAgDasVdR0X3jhhemQQw5JBx10UL591VVXpbvvvjtde+216cQTT2zu4kGDGABk8dsX1dULvk8W5rkAALQeLT50z5gxI40cOTKddNJJNfe1b98+7bDDDunJJ59s1rJBQ7w/Y0ZucrLffvs1d1Goo6pq4fdJ6f8GfAEAoG1q8aF74sSJafbs2al379r9OOL2v/71r3qfM3369Hwp++yzz/L1pEmT0uJqypQp+Xrk5MlpyuzZzV0cmtCTn32W5qSUjlxhhbRm9+7NXRxSStNmzEhjxo5NnZZYIbVv13mB1jFq9pR07+xP0qszPkzVs///9w2Lj1mzJ+UuY51nTkwd55Qnnpm/t2f992/GqFmfpM9Lvo9b0362b1vnfrZfW+++fvv/Tmzbt63P27P/m8uqq6sX24xWLlepVJrvsu1KDVlqMTZu3Li00korpSfqDCd/wgknpMceeyw9/fTTX3jO6aefns4444xFXFIAAABak7Fjx6aVV165ddd0L7fccqlDhw5pwoQJte6P23369Kn3OdEUPQZeK5szZ076+OOPU69evVK7du3S4nomZZVVVsk7tUePHs1dHApiP7cN9nPbYV+3DfZz22A/tw32c9swqQn2c9RdT548OfXt23e+y7b40N2pU6c0cODA9NBDD6U999yzJkTH7cMPP7ze53Tu3DlfKvXs2TO1BPGh8AXQ+tnPbYP93HbY122D/dw22M9tg/3cNvRYyP289NINmx++xYfuELXWQ4YMSZtttlmem/viiy9OU6dOrRnNHAAAAJpDqwjd++yzT/rwww/TqaeemsaPH5823njjdN99931hcDUAAABYlFpF6A7RlHxuzclbg2gOf9ppp32hWTyti/3cNtjPbYd93TbYz22D/dw22M9tQ+dFvJ9b/OjlAAAAsLhq39wFAAAAgNZK6AYAAICCCN0AAABQEKG7GV1++eVptdVWS126dEmbb755euaZZ+a5/K233prWWWedvPwGG2yQ7rnnnlqPR/f8GMF9xRVXTEsuuWTaYYcd0htvvFHwu2BR7+cDDzwwtWvXrtZl5513Lvhd0JT7edSoUWnvvffOy8f+i2kOF3adtMz9fPrpp3/heI7jn5azn6+++ur0ta99LS2zzDL5En976y7v73Pb2M/+PreOfX377bfnaYh79uyZunXrlmdF+v3vf19rGcd029jPBzblMR0DqbHo3XzzzaVOnTqVrr322tKoUaNKhxxySKlnz56lCRMm1Lv8448/XurQoUNp2LBhpddee610yimnlJZYYonSK6+8UrPMeeedV1p66aVLd955Z+mll14q/c///E9p9dVXL33++eeL8J1R9H4eMmRIaeeddy69//77NZePP/54Eb4rFnY/P/PMM6Xjjz++9Mc//rHUp0+f0kUXXbTQ66Rl7ufTTjuttN5669U6nj/88MNF8G5oqv38/e9/v3T55ZeXXnjhhVJVVVXpwAMPzH+L33vvvZpl/H1uG/vZ3+fWsa8feeSR0u23355/h7355puliy++OP82u++++2qWcUy3jf08pAmPaaG7mXzlK18pHXbYYTW3Z8+eXerbt2/p3HPPrXf573znO6Xddtut1n2bb7556dBDD83/nzNnTv5Rd8EFF9Q8/umnn5Y6d+6cf/DROvZz+Qtgjz32KLDUFL2fK/Xr16/eMLYw66Tl7OcI3RtttFGTl5UFt7DH3qxZs0rdu3cvXX/99fm2v89tYz8Hf58XT03x93STTTbJFSHBMd029nNTH9OalzeDGTNmpJEjR+amKGXt27fPt5988sl6nxP3Vy4fBg8eXLP8O++8k8aPH19rmaWXXjo3rZjbOml5+7ns0UcfTSussELq379/Gjp0aProo48KehcUsZ+bY50snCL3STRJ7Nu3b/rSl76U9t133zRmzJgmKDHNtZ+rq6vTzJkz07LLLptv+/vcNvZzmb/PrWtfRwXlQw89lEaPHp223nrrfJ9jum3s56Y+poXuZjBx4sQ0e/bs1Lt371r3x+04iOsT989r+fJ1Y9ZJy9vPIfqS3HDDDfnL4fzzz0+PPfZY2mWXXfJr0TL2c3Osk4VT1D6JH2nDhw9P9913X7ryyivzj7noNzp58uQmKDXNsZ9/9rOf5ZMo5R9//j63jf0c/H1uPfv6s88+S0sttVTq1KlT2m233dKll16adtxxx/yYY7pt7OemPqY7NvoZQLP67ne/W/P/GGhtww03TGussUY+E7f99ts3a9mAxok/3mVxLEcI79evX7rlllvSwQcf3Kxlo/HOO++8dPPNN+fv4xjIh7a1n/19bj26d++eXnzxxTRlypQcuI499tjcGmmbbbZp7qKxCPdzUx7TarqbwXLLLZc6dOiQJkyYUOv+uN2nT596nxP3z2v58nVj1knL28/1iS+HeK0333yziUpO0fu5OdbJwllU+yRGUV177bUdzy1wP//yl7/MYeyBBx7IP8zK/H1uG/u5Pv4+t9x9HU2T11xzzTyi9XHHHZe+9a1vpXPPPTc/5phuG/u5qY9pobsZRBOGgQMH5jMqZXPmzMm3Bw0aVO9z4v7K5cOIESNqll999dXzh6pymUmTJqWnn356ruuk5e3n+rz33nu5f0lMW0HL2M/NsU4WzqLaJ3G2/a233nI8t7D9PGzYsHTWWWflbgIxBU0lf5/bxn6uj7/Pree7O54zffr0/H/HdNvYz01+TDfJcGws0LD2Mcrh8OHD81D1P/rRj/Kw9uPHj8+P77///qUTTzyx1lRSHTt2LP3yl7/MU1XEiLf1TRkW6/jLX/5Sevnll/Noe6YvaF37efLkyXkKoieffLL0zjvvlB588MHSpptuWlprrbVK06ZNa7b32dY1dj9Pnz49TzsTlxVXXDHv0/j/G2+80eB10jr283HHHVd69NFH8/Ecx/8OO+xQWm655UoffPBBs7xHGr+f429vTFNz22231ZpWJr6vK5fx97l172d/n1vPvj7nnHNKDzzwQOmtt97Ky8dvsvhtdvXVV9cs45hu/ft5chMf00J3M7r00ktLq666av4Sj2Hun3rqqZrHvv71r+dh6ivdcsstpbXXXjsvH/O63n333bUejykMfv7zn5d69+6dP3Tbb799afTo0Yvs/VD8fq6uri7ttNNOpeWXXz6H8ZiGKOYhFMRa1n6OL+8451n3Ess1dJ20jv28zz775EAe61tppZXy7ZgvlJazn+N7uL79HCdNy/x9bv372d/n1rOvTz755NKaa65Z6tKlS2mZZZYpDRo0KAe6So7p1r+fq5v4mG4X/zS+fhwAAACYH326AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBoA1o165duvPOO5u7GGn8+PFpxx13TN26dUs9e/Zs7uIAQOGEbgCYiwMPPDCH1bh06tQprbnmmunMM89Ms2bNSour008/PW288cZfuP/9999Pu+yyS2puF110US7Liy++mF5//fXmLg4AFK5j8S8BAC3XzjvvnK677ro0ffr0dM8996TDDjssLbHEEumkk076wrIzZszI4bw5lEqlNHv27Lk+3qdPn7Q4eOutt9LAgQPTWmut1dxFAYBFQk03AMxD586dc2Dt169fGjp0aNphhx3SX//615qa8D333DP94he/SH379k39+/fP97/yyitpu+22S0suuWTq1atX+tGPfpSmTJlSs87y884444y0/PLLpx49eqQf//jHObSXRcg/8sgj0worrJC6dOmSttpqq/Tss8/WPP7oo4/mGvh77703h9go54033pjX+dJLL9XU0A8fPrze5uUNLeMvf/nLtOKKK+Zl4oTDzJkz57m9rrzyyrTGGmvkkw+xPX7/+9/XPLbaaqulP//5z+mGG27I5YnXmJtrr702rbfeevl9xesffvjhNY9deOGFaYMNNshN1FdZZZX0k5/8pFbZ//3vf6fdd989LbPMMnmZWE+cMCl79dVXc63/UkstlXr37p3233//NHHixJrHb7vttrz+8raJfT516tR5vm8AmBuhGwAaIYJYZTh+6KGH0ujRo9OIESPSXXfdlcPZ4MGDc+CLkHzrrbemBx98sFZoLD+vqqoqh+c//vGP6fbbb8+BueyEE07IAfX6669Pzz//fG7aHuv9+OOPa63nxBNPTOedd15eV/SVPu6443LIjCbccdlnn32+8B4aWsZHHnkk10zHdZQjAnw5xNfnjjvuSEcddVQuQwTbQw89NB100EH5+SFeK1oOfOc738ll+/Wvfz3X4B4BP04ExMmBOMkR77+sffv26ZJLLkmjRo3K5Xr44Yfz9iqL58ZJi7///e/5+eeff34O2OHTTz/NJxs22WST9Nxzz6X77rsvTZgwIZcpRLm+973vpR/84Ac1+2evvfbKLQkAYIGUAIB6DRkypLTHHnvk/8+ZM6c0YsSIUufOnUvHH398zeO9e/cuTZ8+veY5v/3tb0vLLLNMacqUKTX33X333aX27duXxo8fX/O8ZZddtjR16tSaZa688srSUkstVZo9e3Z+7hJLLFG66aabah6fMWNGqW/fvqVhw4bl24888kikwNKdd95Zq8ynnXZaaaONNvrCe4ll77jjjkaVsV+/fqVZs2bVLPPtb3+7tM8++8x1e331q18tHXLIIbXui+fsuuuuNbdje8a65yXe58knn1xqqFtvvbXUq1evmtsbbLBB6fTTT6932bPOOqu000471bpv7NixefuMHj26NHLkyPz/d999t8GvDwDzoqYbAOYhaq+jljSaeEeT5Kg5jsHKyqIZcmU/7qgd3WijjXKz5rItt9wyzZkzJ9eIl8UyXbt2rbk9aNCg3ER67NixuXY5mnHH88qiH/lXvvKVvP5Km222WaPfU0PLGDXmHTp0qLkdzbw/+OCDea63sszl9dYt87zE+seNG5e23377uS4TtfLx+EorrZS6d++em4d/9NFHqbq6Oj8ezfLPPvvs/NqnnXZaevnll2ueG03vo+Y99mn5ss466+THYrvHdol1x3799re/na6++ur0ySefNLj8AFCX0A0A87DtttvmkbbfeOON9Pnnn+fmzJVhtfL/zaHI14+gXyn6YUcwL7r5/ry8++676Rvf+EbacMMNc/P7kSNHpssvvzw/Vm72/8Mf/jC9/fbbOYxH8/I4MXHppZfmx+LERvT3jn1aeYn9u/XWW+eTDNFVIPrKr7vuuvl50Tf9nXfeKfR9A9B6Cd0AMJ9QG/2JV1111dSx4/wn/RgwYECuTa0ceOvxxx/P/ZDLA62FWCZCfNlTTz2Va11jYLDyQGTxvLKo+Y4+0REE5yWeN69RzBtTxsaK9VaWubze+ZW5UtRcx4Br0ee9PhGyI/j/6le/SltssUVae+21c814XbEdY3C66CsffcyjxjpsuummuS94vEbs18pL+QRGnFyIWvLoY//CCy/kbRr91QFgQQjdANCE9t1339wUfciQIXkwsWjKfMQRR+Ra1xgpuyxqZQ8++OD02muv5ZG1oxl0DGQWwTfCX4yU/tOf/jQP9BXLHHLIIbn5dDxnXiJMRq1s1N7GiNwxoNiClrGxorwx0FoMhBY1xzHKeITe448/vlHrieb7EapjsLRYTwwkV66pjnAcJyDidtRmx+joV111Va3nH3300en+++/P2yGeG+8vTgiUB1mLwehisLQ4iRFNymPZGPAtTlY8/fTT6ZxzzsmDrI0ZMyaX/8MPP6x5PgA0ltANAE0o+mlHiItg9+Uvfzl961vfyn2EL7vsslrLxX0xV3U0aY5+4v/zP/9Tq694jEi+99575yActbNvvvlmXm+MOD4v8ZwYITyaxcd0ZDEy+oKWsbFiirEYkTymGYv+4L/5zW/yHOfbbLNNo9YTJwMuvvjidMUVV+T1RHPyCN8h+lxHmI8Ryddff/100003pXPPPbfW8yM8R7iOoBzbImrDY10hpnaL2vdYZqeddsp9tyOk9+zZM5/wiOnbYtTzXXfdNT/vlFNOyScAoj8/ACyIdjGa2gI9EwBYIDE/dUxdVTlvNgDQOqnpBgAAgIII3QAAAFAQzcsBAACgIGq6AQAAoCBCNwAAABRE6AYAAICCCN0AAABQEKEbAAAACiJ0AwAAQEGEbgAAACiI0A0AAAAFEboBAAAgFeP/AUpiH5VV6k/DAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "# Import necessary libraries\n", "import pandas as pd\n", @@ -80,6 +126,8 @@ "TRACE_SUCCESS = 0.20\n", "SECONDARY_TRACE_THRESHOLD = 2\n", "\n", + "np.random.seed(50) # set random seed to make code reproducible\n", + "\n", "def simulate_event(m):\n", " \"\"\"\n", " Simulates the infection and tracing process for a series of events.\n", @@ -130,7 +178,7 @@ "\n", " return p_wedding_infections, p_wedding_traces\n", "\n", - "# Run the simulation 1000 times\n", + "# Run the simulation 10 / 100 / 1000 times\n", "results = [simulate_event(m) for m in range(1000)]\n", "props_df = pd.DataFrame(results, columns=[\"Infections\", \"Traces\"])\n", "\n", @@ -193,7 +241,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "sampling-env", "language": "python", "name": "python3" }, @@ -207,7 +255,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.0" + "version": "3.11.13" } }, "nbformat": 4,