diff --git a/env_household_contact_tracing.yml b/env_household_contact_tracing.yml index 86886fd..c9ecc36 100644 --- a/env_household_contact_tracing.yml +++ b/env_household_contact_tracing.yml @@ -13,3 +13,6 @@ dependencies: - pandas - pyvis - beautifulsoup4 + - statsmodels + - torch + - ax diff --git a/examples/calibration.ipynb b/examples/calibration.ipynb new file mode 100644 index 0000000..0b07397 --- /dev/null +++ b/examples/calibration.ipynb @@ -0,0 +1,635 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "source": [ + "import household_contact_tracing.branching_process_models as bpm\r\n", + "from household_contact_tracing.branching_process_controller import BranchingProcessController\r\n", + "from household_contact_tracing.calibration import StandardCalibrationHouseholdLevelTracing\r\n", + "import seaborn as sns\r\n", + "import matplotlib.pyplot as plt" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "# Calibration tutorial" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "This is a quick tutorial on how to calibrate the epidemic models. Household structured branching processes are the underlying mathematical framework, and these have several key metrics that are commonly calibrated:\r\n", + "* The household secondary attack rate (probability that a non-index case in a household gets infected)\r\n", + "* The growth rate of the epidemic\r\n", + "\r\n", + "The household secondary attack rate typically only depends upon the parameter 'household_pairwise_survival_prob', however this could always be varied. The growth rate of the epidemic is more complicated, and this could depends on every parameter in the model. Typically, we calibrate the epidemic model to a growth rate when there is no contact tracing - an unconstrained epidemic model. Without contact tracing, the epidemic growth rate will depend only upon the household secondary attack rate, the outside household infectivity scaling, and the probability that a case reports their infection and isolates (which depends on the symptomatic proportion of cases)." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Evaluating metrics\r\n", + "\r\n", + "By default, models will be initialised with a statistics view, which is able to estimate several metrics in a model, inclduing the growth rate and the household secondary attack rate.\r\n", + "\r\n", + "We configure some model parameters, set up a model, run the simulation for 25 days, and then use the statistics view to print some summaries of the simulated epidemic." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 3, + "source": [ + "params = {'outside_household_infectivity_scaling': 0.1,\r\n", + " 'contact_tracing_success_prob': 0.0, # doesn't matter, no tracing\r\n", + " 'overdispersion': 0.32,\r\n", + " 'asymptomatic_prob': 0.2,\r\n", + " 'asymptomatic_relative_infectivity': 0.35,\r\n", + " 'infection_reporting_prob': 0,\r\n", + " 'contact_trace': False,\r\n", + " 'test_delay': 2,\r\n", + " 'contact_trace_delay': 1,\r\n", + " 'incubation_period_delay': 5,\r\n", + " 'symptom_reporting_delay': 1,\r\n", + " 'household_pairwise_survival_prob': 0.2,\r\n", + " 'do_2_step': False, # doesn't matter, no tracing\r\n", + " 'reduce_contacts_by': 0.3,\r\n", + " 'prob_has_trace_app': 0, # doesn't matter, no tracing\r\n", + " 'hh_propensity_to_use_trace_app': 1, # doesn't matter, no tracing\r\n", + " 'test_before_propagate_tracing': True, # doesn't matter, no tracing\r\n", + " 'starting_infections': 100, \r\n", + " 'node_will_uptake_isolation_prob': 1, # doesn't matter, no tracing\r\n", + " 'self_isolation_duration': 0, # doesn't matter, no tracing\r\n", + " 'quarantine_duration': 0, # doesn't matter, no tracing\r\n", + " 'transmission_probability_multiplier': 1,\r\n", + " 'propensity_imperfect_quarantine': 0, # doesn't matter no tracing\r\n", + " 'global_contact_reduction_imperfect_quarantine': 0, # doesn't matter, no tracing\r\n", + "\r\n", + " }" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 4, + "source": [ + "controller = BranchingProcessController(bpm.HouseholdLevelTracing(params))\r\n", + "controller.graph_view.set_display(False) # don't draw the network\r\n", + "controller.csv_view.set_display(False) # don't save csv info\r\n", + "controller.run_simulation(25, max_active_infections=10000)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "When estimating the growth rate, it is good practice to discard the first 10 or so days. It can take some time for the distribution of infectious ages to mix properly, and for contact tracing to be initialised." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 5, + "source": [ + "controller.statistics_view.growth_rate_summary(discard_first_n_days=10)" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "15 time periods were used to estimate the growth rate.\n", + "The estimated growth rate was 8.57% (95.0% CI: 7.5-9.64%) per day.\n", + "The estimated doubling time is 8.43 (95.0% CI: 7.53-9.59) days.\n" + ] + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 6, + "source": [ + "controller.statistics_view.household_secondary_attack_rate_summary(use_first_generation_only=False)" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Household secondary attack rate summary:\n", + "173 local household epidemics were eligible to be included.\n", + "All households with completed local epidemics were included. This may lead to a biased sample, as it is possible that local epidemics with a long duration were not included.\n", + "There were 312 non-index susceptible individuals exposed, of which 193 were infected.\n", + "This yields a household secondary attack rate of 62% (95% CI: 56-67%).\n" + ] + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Household secondary attack rate simulations" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Only households with completed local epidemics are included when estimating the household secondary attack rate. This however can lead to biased samples as households with long lasting epidemics might be excluded, and these may be likely to have larger final sizes.\r\n", + "\r\n", + "To get around this a second simulation has been implemented that simulates only the local epidemics, and does not simulate outside household transmission. The simulation is run until all the starting households have completed their local epidemics. This will yield a sample that can be used to obtain an unbiased estiamted of the household secondary attack rate. It is necessary to use a larger number of starting infections to obtain the household secondary attack rate." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 7, + "source": [ + "params['starting_infections'] = 500 # increasing the number of starting infections will increase the sample size\r\n", + "controller = BranchingProcessController(bpm.HouseholdLevelTracing(params))\r\n", + "controller.graph_view.set_display(False) # don't draw the network\r\n", + "controller.csv_view.set_display(False) # don't save csv info\r\n", + "controller.run_hh_sar_simulation()" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 8, + "source": [ + "controller.statistics_view.household_secondary_attack_rate_summary(use_first_generation_only=True)" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Household secondary attack rate summary:\n", + "499 local household epidemics were eligible to be included.\n", + "Only the first generation of the household epidemic was included in this calculation.\n", + "There were 1053 non-index susceptible individuals exposed, of which 809 were infected.\n", + "This yields a household secondary attack rate of 77% (95% CI: 74-79%).\n" + ] + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "The previous estimate of the household secondary attack rate was around 60%, however using this alternative simulation method we find that the secondary attack rate is estimated to be much closer to 80%. This is why it is important to use the correct method when estimating the household secondary attack rate." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Calibration" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "We have seen that it is necessary to run two different types of simulations in our to compute two metrics required for calibrating the model. To make our lives easier, we have written several configuration classes. These classes contain methods that will run the model, evaluate the fit of the model, and find the optimal combination of parameters." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Suppose we want a model with a household secondary attack rate of 25%, and a growth rate of 10% per day." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 9, + "source": [ + "calibrator = StandardCalibrationHouseholdLevelTracing(\r\n", + " desired_growth_rate=0.1, \r\n", + " desired_hh_sar=0.2, \r\n", + " asymptomatic_prob=0.2, \r\n", + " asymptomatic_relative_infectivity=0.35,\r\n", + " infection_reporting_prob=0.2,\r\n", + " reduce_contacts_by=0.5)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "We first perform a trial for two proposed parameters." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 10, + "source": [ + "calibrator.eval_metrics(\r\n", + " household_pairwise_survival_prob = 0.9,\r\n", + " outside_household_infectivity_scaling = 0.2\r\n", + ")" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Estimating growth rate using 30 time periods\n" + ] + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "{'growth_rate': 0.0034215244104689423, 'hh_sar': 0.10606060606060606}" + ] + }, + "metadata": {}, + "execution_count": 10 + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "So, we were quite far away. The calibrator class will evaluate the performance of these parameters." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 11, + "source": [ + "calibrator.evaluate_fit(\r\n", + " household_pairwise_survival_prob = 0.9,\r\n", + " outside_household_infectivity_scaling = 0.2\r\n", + ")" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Estimating growth rate using 30 time periods\n" + ] + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "0.18708774795757915" + ] + }, + "metadata": {}, + "execution_count": 11 + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "The household secondary attack rate is too low, we could decrease the household pairwise survival probability to increase the household secondary attack rate." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 12, + "source": [ + "calibrator.eval_metrics(\r\n", + " household_pairwise_survival_prob = 0.8,\r\n", + " outside_household_infectivity_scaling = 0.2\r\n", + ")" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Estimating growth rate using 30 time periods\n" + ] + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "{'growth_rate': 0.0653517347589822, 'hh_sar': 0.22484134179510426}" + ] + }, + "metadata": {}, + "execution_count": 12 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 13, + "source": [ + "calibrator.evaluate_fit(\r\n", + " household_pairwise_survival_prob = 0.8,\r\n", + " outside_household_infectivity_scaling = 0.2\r\n", + ")" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Estimating growth rate using 30 time periods\n" + ] + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "0.066385683081691" + ] + }, + "metadata": {}, + "execution_count": 13 + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "The fit has improved, as the return value is smaller. Moving on, let's use the optimisation routine to find the optimal combination of parameters. To do this, we need to supply a range of parameters. These could be set fairly wide, as the expense of requiring more simulation runs to find the optimal combination of parameters. A machine learning hyperparameter optimisation framework will handle the rest. Some of the runs could take a while, particularly if they overshoot the growth rate, as the computation time is related to the size of the epidemic." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 14, + "source": [ + "calibrator.optimise(\r\n", + " household_pairwise_survival_prob_range = [0.75, 0.95],\r\n", + " outside_household_infectivity_scaling_range = [0.2, 0.4],\r\n", + " state_criteria={\r\n", + " 'max_time': 25,\r\n", + " 'infection_threshold': 1e5\r\n", + " }\r\n", + ")" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stderr", + "text": [ + "[INFO 08-06 10:14:42] ax.modelbridge.dispatch_utils: Using GPEI (Bayesian optimization) since there are more continuous parameters than there are categories for the unordered categorical parameters.\n", + "[INFO 08-06 10:14:42] ax.modelbridge.dispatch_utils: Using Bayesian Optimization generation strategy: GenerationStrategy(name='Sobol+GPEI', steps=[Sobol for 5 trials, GPEI for subsequent trials]). Iterations after 5 will take longer to generate due to model-fitting.\n", + "[INFO 08-06 10:14:42] ax.service.managed_loop: Started full optimization with 20 steps.\n", + "[INFO 08-06 10:14:42] ax.service.managed_loop: Running optimization trial 1...\n", + "[INFO 08-06 10:14:49] ax.service.managed_loop: Running optimization trial 2...\n", + "[INFO 08-06 10:14:56] ax.service.managed_loop: Running optimization trial 3...\n", + "[INFO 08-06 10:15:05] ax.service.managed_loop: Running optimization trial 4...\n", + "[INFO 08-06 10:15:10] ax.service.managed_loop: Running optimization trial 5...\n", + "[INFO 08-06 10:15:16] ax.service.managed_loop: Running optimization trial 6...\n", + "[INFO 08-06 10:15:24] ax.service.managed_loop: Running optimization trial 7...\n", + "[INFO 08-06 10:15:29] ax.service.managed_loop: Running optimization trial 8...\n", + "[INFO 08-06 10:15:32] ax.service.managed_loop: Running optimization trial 9...\n", + "[INFO 08-06 10:15:37] ax.service.managed_loop: Running optimization trial 10...\n", + "[INFO 08-06 10:15:41] ax.service.managed_loop: Running optimization trial 11...\n", + "[INFO 08-06 10:15:47] ax.service.managed_loop: Running optimization trial 12...\n", + "[INFO 08-06 10:15:52] ax.service.managed_loop: Running optimization trial 13...\n", + "[INFO 08-06 10:16:02] ax.service.managed_loop: Running optimization trial 14...\n", + "[INFO 08-06 10:16:06] ax.service.managed_loop: Running optimization trial 15...\n", + "[INFO 08-06 10:16:10] ax.service.managed_loop: Running optimization trial 16...\n", + "[INFO 08-06 10:16:15] ax.service.managed_loop: Running optimization trial 17...\n", + "[INFO 08-06 10:16:21] ax.service.managed_loop: Running optimization trial 18...\n", + "[INFO 08-06 10:16:27] ax.service.managed_loop: Running optimization trial 19...\n", + "[INFO 08-06 10:16:31] ax.service.managed_loop: Running optimization trial 20...\n" + ] + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "({'outside_household_infectivity_scaling': 0.2369218112772867,\n", + " 'household_pairwise_survival_prob': 0.8104007705085332},\n", + " ({'objective': 0.026804953615459474},\n", + " {'objective': {'objective': 7.980087959078822e-06}}))" + ] + }, + "metadata": {}, + "execution_count": 14 + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "The optimal parameter combinations can be recovered using:" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 15, + "source": [ + "calibrator.best_parameters" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "{'outside_household_infectivity_scaling': 0.2369218112772867,\n", + " 'household_pairwise_survival_prob': 0.8104007705085332}" + ] + }, + "metadata": {}, + "execution_count": 15 + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Finally, we can easily sample the fitting metrics to check that the fit looks good." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 17, + "source": [ + "outputs = calibrator.get_fitted_model_metric_samples(n_obs = 20, state_criteria={'max_time': 25})" + ], + "outputs": [ + { + "output_type": "error", + "ename": "TypeError", + "evalue": "get_fitted_model_metric_samples() got an unexpected keyword argument 'state_criteria'", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mTypeError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0moutputs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcalibrator\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_fitted_model_metric_samples\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_obs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstate_criteria\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m{\u001b[0m\u001b[1;34m'max_time'\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;36m25\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;31mTypeError\u001b[0m: get_fitted_model_metric_samples() got an unexpected keyword argument 'state_criteria'" + ] + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, + "source": [ + "growth_rates = [\r\n", + " output['growth_rate'] for output in outputs\r\n", + "]\r\n", + "\r\n", + "hh_sars = [\r\n", + " output['hh_sar'] for output in outputs\r\n", + "]\r\n", + "\r\n", + "outputs" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "[{'growth_rate': 0.11082584992572188, 'hh_sar': 0.18533157663592445},\n", + " {'growth_rate': 0.10738751955334369, 'hh_sar': 0.21236059479553904},\n", + " {'growth_rate': 0.09084650192266462, 'hh_sar': 0.19239274657231314},\n", + " {'growth_rate': 0.10075203747153494, 'hh_sar': 0.2336119665640123},\n", + " {'growth_rate': 0.10735735237940272, 'hh_sar': 0.2076707202993452},\n", + " {'growth_rate': 0.08313646738902405, 'hh_sar': 0.23312331233123312},\n", + " {'growth_rate': 0.08981749442949848, 'hh_sar': 0.2202781516375056},\n", + " {'growth_rate': 0.09322273991873418, 'hh_sar': 0.21407121407121407},\n", + " {'growth_rate': 0.09289832274649205, 'hh_sar': 0.21107580571947346},\n", + " {'growth_rate': 0.09344663727571464, 'hh_sar': 0.20498614958448755},\n", + " {'growth_rate': 0.096859550802151, 'hh_sar': 0.21289151157512484},\n", + " {'growth_rate': 0.10469928901617155, 'hh_sar': 0.2064975522919448},\n", + " {'growth_rate': 0.10857606282994096, 'hh_sar': 0.20827710301394511},\n", + " {'growth_rate': 0.09836023406359318, 'hh_sar': 0.20940959409594095},\n", + " {'growth_rate': 0.08472873842505456, 'hh_sar': 0.20318181818181819},\n", + " {'growth_rate': 0.0941299371383913, 'hh_sar': 0.2197851387645479},\n", + " {'growth_rate': 0.09565612682725907, 'hh_sar': 0.20922795797167656},\n", + " {'growth_rate': 0.07646034952842615, 'hh_sar': 0.20435967302452315},\n", + " {'growth_rate': 0.09383057360324933, 'hh_sar': 0.17884702678166137},\n", + " {'growth_rate': 0.08110175647926467, 'hh_sar': 0.19276556776556777}]" + ] + }, + "metadata": {}, + "execution_count": 18 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, + "source": [ + "sns.displot(growth_rates)" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 26 + }, + { + "output_type": "display_data", + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPUElEQVR4nO3dfYxlB1nH8d9DW9qioKBbhLYroAQlJCoZEIshCsRUNCAEBeILJuhiDATwLRj/kT/8A6OGaIyyCgFRK6gQBaVQlRcNL7IgYmvxBXlpaW23GgVjgm15/GNudVzb3en0nnnm7nw+yWTuPfeee56cvfvNmTP33qnuDgD7717TAwAcVgIMMESAAYYIMMAQAQYYcu70ADtdfvnlfeWVV06PAbBudWcLD9QR8C233DI9AsC+OVABBjhMBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMWfTDeKrqE0k+m+T2JLd199aS2wPYJPvxaWjf3N0+ZQfgFE5BAAxZ+gi4k7y9qjrJK7v7+Kl3qKpjSY4lydGjRxceh3W6+NKjueH666bH2JNzzjs/t9/6uekx9uTBl1yaT1/3qekxWINa8s/SV9WDu/uGqrooyVVJXtjd776r+29tbfWJEycWm4f1qqo865XvmR5jT17//Ms2evYl/9+yiP3/QPbuvmH1/eYkb0ry2CW3B7BJFgtwVX1BVd33jstJviXJ1UttD2DTLHkO+IFJ3lRVd2znt7vbH3wDWFkswN39T0m+ZqnHB9h0XoYGMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYsniAq+qcqvqrqnrL0tsC2CT7cQT8oiTX7sN2ADbKogGuqkuSfFuSX19yOwCbaOkj4Fck+Ykkn7+rO1TVsao6UVUnTp48ufA4wKSLLz2aqtrIr4svPbr2/XHu2h9xpaq+PcnN3f3Bqvqmu7pfdx9PcjxJtra2eql5gHk3XH9dnvXK90yPsSevf/5la3/MJY+AH5/kqVX1iSS/k+SJVfWbC24PYKMsFuDu/snuvqS7H5Lk2Un+rLu/Z6ntAWwarwMGGLLYOeCduvudSd65H9sC2BSOgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMWC3BVXVBVf1lVf11V11TVy5baFsAmOnfBx/5ckid2939U1XlJ/qKq3trd71twmwAbY7EAd3cn+Y/V1fNWX73U9gA2zaLngKvqnKr6cJKbk1zV3e+/k/scq6oTVXXi5MmTS44DcKAsGuDuvr27vzbJJUkeW1WPupP7HO/ure7eOnLkyJLjABwo+/IqiO7+tyTvTHL5fmwPYBMs+SqII1X1xavLFyZ5cpKPLrU9gE2z5KsgHpTktVV1TrZD/4bufsuC2wPYKEu+CuIjSb5uqccH2HTeCQcwRIABhggwwBABBhgiwABDdhXgqnr8bpYBsHu7PQL+pV0uA2CXTvs64Kr6hiSXJTlSVT+y46b7JTlnycEAznZneiPGvZN84ep+992x/DNJnrnUUACHwWkD3N3vSvKuqnpNd39yn2YCOBR2+1bk86vqeJKH7Fynu5+4xFAAh8FuA/y7SX41ya8nuX25cQAOj90G+Lbu/pVFJwE4ZHb7MrQ3V9UPV9WDquoBd3wtOhnAWW63R8DPXX3/8R3LOsnD1jsOwOGxqwB390OXHgTgsNlVgKvq++5seXf/xnrHATg8dnsK4jE7Ll+Q5ElJPpREgAH2aLenIF6483pVfVGS1y0yEcAhsdePo/zPJA9f5yAAh81uzwG/Oduveki2P4Tnq5O8YamhAA6D3Z4D/rkdl29L8snuvn6BeQAOjV2dglh9KM9Hs/2JaPdP8l9LDgVwGOz2L2J8V5K/TPKdSb4ryfurysdRAtwDuz0F8VNJHtPdNydJVR1J8idJfm+pwQDOdrt9FcS97ojvyr/cjXUBuBO7PQK+sqreluSK1fVnJfnjZUYCOBzO9DfhvjLJA7v7x6vqGUm+MUkleW+S39qH+QDOWmc6jfCKJJ9Nku5+Y3f/SHe/JNtHv69YdjSAs9uZAvyQ7v7IqQu7+0S2/zwRAHt0pgBfcJrbLlznIACHzZkC/IGq+sFTF1bV85J8cJmRAA6HM70K4sVJ3lRV353/De5WknsnefqCcwGc9U4b4O6+KcllVfXNSR61WvxH3f1ni08GcJbb7ecBvyPJOxaeBeBQ8W42gCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhiwW4qi6tqndU1bVVdU1VvWipbQFsol39SaI9ui3Jj3b3h6rqvkk+WFVXdfffLrhNgI2x2BFwd9/Y3R9aXf5skmuTXLzU9gA2zZJHwP+jqh6S5OuSvP9ObjuW5FiSHD16dE+Pf/GlR3PD9dfdgwnnnHPe+bn91s9Nj8Emude5qarpKViDxQNcVV+Y5PeTvLi7P3Pq7d19PMnxJNna2uq9bOOG66/Ls175nns055TXP/+yjZ6dAZ+/zXPmLLHoqyCq6rxsx/e3uvuNS24LYNMs+SqISvKqJNd29y8stR2ATbXkEfDjk3xvkidW1YdXX09ZcHsAG2Wxc8Dd/RdJ/KYA4C54JxzAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGDIYgGuqldX1c1VdfVS2wDYZEseAb8myeULPj7ARlsswN397iT/utTjA2w654ABhowHuKqOVdWJqjpx8uTJ6XEA9s14gLv7eHdvdffWkSNHpscB2DfjAQY4rJZ8GdoVSd6b5BFVdX1VPW+pbQFsonOXeuDufs5Sjw1wNnAKAmCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYIMMAQAQYYIsAAQwQYYIgAAwwRYIAhAgwwRIABhggwwBABBhgiwABDBBhgiAADDBFggCECDDBEgAGGCDDAEAEGGCLAAEMEGGCIAAMMEWCAIQIMMESAAYYsGuCquryq/q6q/rGqXrrktgA2zWIBrqpzkvxykm9N8sgkz6mqRy61PYBNs+QR8GOT/GN3/1N3/1eS30nytAW3B7BRqruXeeCqZya5vLt/YHX9e5N8fXe/4JT7HUtybHX1EUn+LsmXJrllkcHWx4zrYcb12IQZk82Yc4kZb+nuy09deO6aN7JT3cmy/1f77j6e5Pj/WbHqRHdvLTXYOphxPcy4HpswY7IZc+7njEuegrg+yaU7rl+S5IYFtwewUZYM8AeSPLyqHlpV907y7CR/uOD2ADbKYqcguvu2qnpBkrclOSfJq7v7ml2ufvzMdxlnxvUw43pswozJZsy5bzMu9ks4AE7PO+EAhggwwJDFA3ymtyPXtl9c3f6Rqnr0avkjqurDO74+U1UvXt3201X16R23PWVixtVtL6mqa6rq6qq6oqouWC1/QFVdVVX/sPp+/wM440Hajy9azXfNHf/Oq+Vr3Y8Lzrnf+/Krquq9VfW5qvqx3aw78Jzcy4wHaT++uqpurqqrT1m+vv3Y3Yt9ZfuXbx9L8rAk907y10keecp9npLkrdl+3fDjkrz/Lh7nn5N8+er6Tyf5sekZk1yc5ONJLlxdf0OS719d/tkkL11dfmmSlx/AGQ/KfnxUkquT3Cfbvxj+kyQPX/d+XHjO/d6XFyV5TJKf2bnd06078Jzcy4wHYj+ubntCkkcnufqU5Wvbj0sfAe/m7chPS/Ibve19Sb64qh50yn2elORj3f3JAzjjuUkurKpzs/0f84Yd67x2dfm1Sb7jAM64Tvdkxq9O8r7u/s/uvi3Ju5I8fcc669qPS865Tmecsbtv7u4PJLn1bqy7r8/JPc64TvdkxnT3u5P865087tr249IBvjjJdTuuX79adnfv8+wkV5yy7AWrHw9ffQ9/lNrzjN396SQ/l+RTSW5M8u/d/fbVfR7Y3Tcmyer7RQdwxuQA7MdsH1U+oaq+pKruk+0j0DvexLPO/bjknMn+7su9rLvfz8m9rnsQ9uPprG0/Lh3g3bwd+bT3qe03cTw1ye/uuP1XknxFkq/NdlR+fmLG1ZPjaUkemuTBSb6gqr7nHsxyV5aa8UDsx+6+NsnLk1yV5Mps/6h42z2Y5XSWmnO/9+US694dS814UPbjvlg6wLt5O/KZ7vOtST7U3TfdsaC7b+ru27v780l+Lds/akzM+OQkH+/uk919a5I3JrlsdZ+b7jgFsPp+80Gb8QDtx3T3q7r70d39hGz/2PcPq/uscz8uNufAvtzLuvv9nLzb6x6g/Xg6a9uPSwd4N29H/sMk37f6zfPjsv0j8o07bn9OTjn9cMo54qdn+0fDiRk/leRxVXWfqqpsn6u+dsc6z11dfm6SPzhoMx6g/Ziqumj1/WiSZ+R//83XuR8Xm3NgX+5l3f1+Tt7tdQ/Qfjyd9e3Hvf72brdf2T5P9vfZ/m3kT62W/VCSH1pdrmx/cPvHkvxNkq0d694nyb8k+aJTHvN1q/t+ZLUzHjQ448uSfDTbT5TXJTl/tfxLkvxpto+Q/jTJAw7gjAdpP/55kr/N9o/1T9qxfK37ccE593tfflm2j/A+k+TfVpfvd1frDj0n9zLjQdqPV2T7NMitq+XPW/d+9FZkgCHeCQcwRIABhggwwBABBhgiwABDBBhgiAADDPlvTa2yWoFd9H0AAAAASUVORK5CYII=" + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, + "source": [ + "sns.displot(hh_sars)" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 27 + }, + { + "output_type": "display_data", + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAP9ElEQVR4nO3de4yldX3H8feHXRAULLZdrVl2s5gaqzVeyHhbDIlompWaGhsrGrWloV1MowGtNm1Mm/pH/2hiGkzTtGzQeqOKF0iqtihGUBt0dUA0i4utopYttjt4KRgTYeHbP87ZdFz3cnZmnvnOM/N+JSd7bs/5/X5zhjfPPueyqSokSavvlO4JSNJGZYAlqYkBlqQmBliSmhhgSWqyuXsCi+3atatuuOGG7mlI0krL0a5cU3vA9957b/cUJGnVrKkAS9JGYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJoMGOMnZST6S5M4k+5M8b8jxJGlMhv42tHcAN1TVy5OcBjxy4PEkaTQGC3CSRwMXAJcAVNUDwANDjSdJYzPkIYgnAAvAPyb5SpKrkzzqyDsl2Z1kPsn8wsLCgNPRerJ123aStJ22btve/SPQOpCh/ln6JHPAF4Hzq2pvkncA91XVnx9rm7m5uZqfnx9kPlpfknDxVbe0jX/tZTsZ6r8drUur/oXsB4ADVbV3evkjwHkDjidJozJYgKvqv4G7kzxpetULga8PNZ4kjc3Q74J4A3DN9B0QdwG/P/B4kjQagwa4qm4H5oYcQ5LGyk/CSVITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNNg/54Em+A9wPPAQcqqq5IceTpDEZNMBTL6iqe1dhHEkaFQ9BSFKToQNcwKeS3Jpk99HukGR3kvkk8wsLCwNPR1ohp2wmSctp67bt3avXChn6EMT5VXVPkscCNya5s6o+t/gOVbUH2AMwNzdXA89HWhkPH+Liq25pGfray3a2jKuVN+gecFXdM/3zIHA98Owhx5OkMRkswEkeleSsw+eB3wD2DTWeJI3NkIcgHgdcn+TwOP9UVTcMOJ4kjcpgAa6qu4CnD/X4kjR2vg1NkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWoyeICTbErylSQfH3osSRqT1dgDvhzYvwrjSNKoDBrgJOcAvwlcPeQ4kjRGQ+8BXwn8CfDwse6QZHeS+STzCwsLA09HktaOwQKc5CXAwaq69Xj3q6o9VTVXVXNbtmwZajqStOYMuQd8PvBbSb4DfBC4MMn7BxxPkkZlsABX1Z9V1TlVtQN4JfCZqnrNUONJ0tj4PmBJarJ5NQapqpuBm1djLEkaC/eAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJanJTAFOcv4s10mSZjfrHvDfznidJGlGm493Y5LnATuBLUnetOimRwObhpyYJK13xw0wcBpw5vR+Zy26/j7g5UNNSpI2guMGuKo+C3w2ybur6rurNCdJ2hBOtAd82COS7AF2LN6mqi4cYlKStBHMGuAPA/8AXA08NNx0JGnjmDXAh6rq7wediSRtMLO+De1jSf4oyeOT/OLh06Azk6R1btY94N+b/vmWRdcV8ISVnY4kbRwzBbiqzh16IpK00cwU4CS/e7Trq+q9KzsdSdo4Zj0E8axF508HXgjcBhhgSVqiWQ9BvGHx5SS/ALxvkBlJ0gax1K+j/AnwxJWciCRtNLMeA/4Yk3c9wORLeJ4MfGioSUnSRjDrMeC3Lzp/CPhuVR0YYD6StGHMdAhi+qU8dzL5RrTHAA+caJskpyf5UpKvJrkjyduWN1VJWl9m/RcxXgF8Cfgd4BXA3iQn+jrKnwIXVtXTgWcAu5I8dxlzlaR1ZdZDEG8FnlVVBwGSbAE+DXzkWBtUVQE/nl48dXqqY91fkjaaWd8Fccrh+E59f5Ztk2xKcjtwELixqvYe5T67k8wnmV9YWJhxOpI0frMG+IYkn0xySZJLgE8A/3Kijarqoap6BnAO8OwkTz3KffZU1VxVzW3ZsuUkpi5J43aifxPuV4HHVdVbkvw28HwgwBeAa2YdpKp+lORmYBewb+nTlaT140R7wFcC9wNU1XVV9aaqeiOTvd8rj7dhki1Jzp6ePwN4EZN3UkiSOPGLcDuq6mtHXllV80l2nGDbxwPvSbKJSeg/VFUfX9o0JWn9OVGATz/ObWccb8NpuJ950jOSpA3iRIcgvpzkD4+8MsmlwK3DTEmSNoYT7QFfAVyf5NX8f3DngNOAlw04L0la944b4Kr6H2BnkhcAh99C9omq+szgM5OkdW7W7wO+Cbhp4LlI0oay1O8DliQtkwGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCaDBTjJtiQ3Jdmf5I4klw81liSN0eYBH/sQ8MdVdVuSs4Bbk9xYVV8fcExJGo3B9oCr6ntVddv0/P3AfmDrUONJ0tisyjHgJDuAZwJ7j3Lb7iTzSeYXFhZWYzrrytZt20nSctq6bXv38jemUza3Pefdz/t6+30f8hAEAEnOBD4KXFFV9x15e1XtAfYAzM3N1dDzWW/uOXA3F191S8vY1162s2XcDe/hQ23POfQ+7+vt933QPeAkpzKJ7zVVdd2QY0nS2Az5LogA7wT2V9XfDDWOJI3VkHvA5wOvBS5Mcvv0dNGA40nSqAx2DLiq/g3IUI8vSWPnJ+EkqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCaDBTjJu5IcTLJvqDEkacyG3AN+N7BrwMeXpFEbLMBV9TngB0M9viSN3ebuCSTZDewG2L59+5IeY+u27dxz4O6VnNbMNp36CB568KctY7c7ZTNJumchjVZ7gKtqD7AHYG5urpbyGPccuJuLr7plRec1q2sv29k29uHx2zx8qPXnLo2d74KQpCYGWJKaDPk2tA8AXwCelORAkkuHGkuSxmiwY8BV9aqhHluS1gMPQUhSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQYNcJJdSb6R5JtJ/nTIsSRpbAYLcJJNwN8BLwaeArwqyVOGGk+SxmbIPeBnA9+sqruq6gHgg8BLBxxPkkYlVTXMAycvB3ZV1R9ML78WeE5Vvf6I++0Gdk8vPgn4xgpO45eBe1fw8Tqtp7WA61nL1tNaYG2s596q2nXklZsHHDBHue7nal9Ve4A9g0wgma+quSEee7Wtp7WA61nL1tNaYG2vZ8hDEAeAbYsunwPcM+B4kjQqQwb4y8ATk5yb5DTglcA/DzieJI3KYIcgqupQktcDnwQ2Ae+qqjuGGu8YBjm00WQ9rQVcz1q2ntYCa3g9g70IJ0k6Pj8JJ0lNDLAkNRllgE/0Eeckv5bkC0l+muTNR9z2xiR3JNmX5ANJTl+9mR/dMtdz+XQtdyS5YtUmfQwzrOXVSb42Pd2S5Omzbtthmet5V5KDSfat7qyPbanrSbItyU1J9k9/1y5f/dn/3FyXupbTk3wpyVena3nb6s9+qqpGdWLygt63gCcApwFfBZ5yxH0eCzwL+CvgzYuu3wp8GzhjevlDwCUjXs9TgX3AI5m8oPpp4IlrfC07gcdMz78Y2DvrtmNaz/TyBcB5wL7OdazQ8/N44Lzp+bOAf+98fpa5lgBnTs+fCuwFntuxjjHuAZ/wI85VdbCqvgw8eJTtNwNnJNnMJFzd701eznqeDHyxqn5SVYeAzwIvW41JH8Msa7mlqn44vfhFJu8Pn2nbBstZD1X1OeAHqzXZGSx5PVX1vaq6bXr+fmA/kx2aLstZS1XVj6fXnzo9tbwbYYwB3grcvejyAWb8Raiq/wLeDvwn8D3gf6vqUys+w5Oz5PUw2fu9IMkvJXkkcBE/++GX1Xaya7kU+NclbrsalrOetWhF1pNkB/BMJnuOXZa1liSbktwOHARurKqWtQz5UeShzPQR56NumDyGyf8lzwV+BHw4yWuq6v0rN72TtuT1VNX+JH8N3Aj8mMlfww6t4NxO1sxrSfICJv9RPP9kt11Fy1nPWrTs9SQ5E/gocEVV3bfiM5zdstZSVQ8Bz0hyNnB9kqdW1aofqx/jHvByPuL8IuDbVbVQVQ8C1zE5TtRpWR/Zrqp3VtV5VXUBk7/u/scKz+9kzLSWJE8DrgZeWlXfP5ltV9ly1rMWLWs9SU5lEt9rquq6ged6Iivy3FTVj4CbgZ/7opxV0XUQfaknJnvtdzHZiz188P3Xj3Hfv+RnX7R6DnAHk2O/Ad4DvGGs65le99jpn9uBO5m+6LBW1zKd5zeBnUv9OYxhPYtu38HaeRFuOc9PgPcCV3avYwXWsgU4e3r+DODzwEta1tH9g1ziD/8iJq/Cfgt46/S61wGvm57/FSb/h7yPyaGGA8Cjp7e9bRqqfcD7gEeMfD2fB74+/QV84QjWcjXwQ+D26Wn+eNt2n5a5ng8wea3hwelzdulY18Pkr+8FfG3RbReNdC1PA74yXcs+4C+61uBHkSWpyRiPAUvSumCAJamJAZakJgZYkpoYYElqYoAlqYkBlqQm/wcX4Uvj3h4kpwAAAABJRU5ErkJggg==" + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [], + "metadata": {} + } + ], + "metadata": { + "orig_nbformat": 4, + "language_info": { + "name": "python", + "version": "3.9.4", + "mimetype": "text/x-python", + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "pygments_lexer": "ipython3", + "nbconvert_exporter": "python", + "file_extension": ".py" + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3.9.4 64-bit (conda)" + }, + "interpreter": { + "hash": "5e088ab515b9ff01afc6092114ac7786286eca1d23e3587660ba5a8d6e06cf28" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} \ No newline at end of file diff --git a/examples/household_level_contact_tracing.ipynb b/examples/household_level_contact_tracing.ipynb index e8ebe7e..68daed7 100644 --- a/examples/household_level_contact_tracing.ipynb +++ b/examples/household_level_contact_tracing.ipynb @@ -129,7 +129,7 @@ "metadata": {}, "outputs": [], "source": [ - "[node for node in controller.model.network.all_nodes() if node.isolated]" + "[node for node in controller.model.network.all_nodes() if node.infection.isolated]" ] }, { @@ -211,7 +211,7 @@ "metadata": {}, "outputs": [], "source": [ - "controller.model.network.node(1).testing_delay" + "controller.model.network.node(1).tracing.testing_delay" ] }, { @@ -339,7 +339,7 @@ " node \n", " for node \n", " in controller.model.network.all_nodes()\n", - " if node.propensity_imperfect_isolation\n", + " if node.tracing_adherence.propensity_imperfect_isolation\n", "]" ] }, @@ -353,7 +353,7 @@ " controller.model.infection.contact_rate_reduction.get_contact_rate_reduction(node)\n", " for node \n", " in controller.model.network.all_nodes()\n", - " if node.propensity_imperfect_isolation\n", + " if node.tracing_adherence.propensity_imperfect_isolation\n", "]" ] }, diff --git a/examples/run_testing_contact_model.py b/examples/run_testing_contact_model.py index 2c986fc..d4d52c9 100644 --- a/examples/run_testing_contact_model.py +++ b/examples/run_testing_contact_model.py @@ -1,14 +1,16 @@ import os -from copy import copy +from copy import copy, deepcopy from household_contact_tracing.branching_process_controller import BranchingProcessController from household_contact_tracing.branching_process_models import IndividualTracingDailyTesting def main(): - #example_1() - #example_2() + example_1() + example_2() example_3() + example_4() + def example_1(): @@ -44,9 +46,9 @@ def prob_testing_positive_function(time_relative_to_symptom_onset): controller.timeline_view.set_display(True) controller.graph_pyvis_view.open_in_browser = True controller.graph_view.set_display(True) - controller.run_simulation(15) + controller.run_simulation({"max_time": 15}) controller.graph_pyvis_view.set_display(False) - controller.run_simulation(20) + controller.run_simulation({"max_time": 20}) @@ -92,7 +94,7 @@ def prob_pcr_positive(infectious_age): controller.shell_view.set_display(False) controller.timeline_view.set_display(True) controller.graph_view.set_display(True) - controller.run_simulation(2) + controller.run_simulation({"max_time": 2}) # Re-initialise and re-run model multiple times and save result to specified file # instead of default ('/temp/sumulation_ouput_[todays date].csv') @@ -102,7 +104,7 @@ def prob_pcr_positive(infectious_age): for idx in range(0, 10): controller.model = IndividualTracingDailyTesting(params) controller.csv_view.display_params = ["household_pairwise_survival_prob", "asymptomatic_relative_infectivity"] - controller.run_simulation(20) + controller.run_simulation({"max_time": 20}) # define some easy to look at test sensitivity functions @@ -151,7 +153,67 @@ def example_3(): controller.graph_view.set_display(True) controller.timeline_view.set_display(True) - controller.run_simulation(16) + controller.run_simulation({"max_time": 16}) + +def example_4(): + params = {"outside_household_infectivity_scaling": 0.3, + "contact_tracing_success_prob": 0.7, + "overdispersion": 0.32, + "asymptomatic_prob": 0.2, + "asymptomatic_relative_infectivity": 0.35, + "infection_reporting_prob": 0.3, + "LFA_testing_requires_confirmatory_PCR": False, + "test_delay": 1, + "contact_trace_delay": 1, + "incubation_period_delay": 5, + "symptom_reporting_delay": 1, + "household_pairwise_survival_prob": 0.2, + "propensity_risky_behaviour_lfa_testing": 0, + "global_contact_reduction_risky_behaviour": 0, + "household_positive_policy": "lfa_testing_no_quarantine" + } + # Copy of test_TestingContactModel.py - test_traced_nodes_are_lateral_flow_tested(simple_model_high_test_prob): + """Checks that a node who is traced is placed under lateral flow testing. + + To do this we: + * Initialises a model with 100% contact tracing success probability + * Create a new infection outside the initial household + * Household 1 traces household 2 with 100% success probability and delay 1 + * Simulate one day twice + + """ + def prob_testing_positive_lfa_func(infectious_age): + if infectious_age in [4, 5, 6]: + return 1 + else: + return 0 + + def prob_testing_positive_pcr_func(infectious_age): + if infectious_age in [4, 5, 6]: + return 0 + else: + return 0 + + params["contact_tracing_success_prob"] = 1 + + model = IndividualTracingDailyTesting(params) + model.intervention.increment_tracing.prob_pcr_positive = prob_testing_positive_pcr_func + model.prob_lfa_positive = prob_testing_positive_lfa_func + + model.infection.new_outside_household_infection(time=0, infecting_node=model.network.node(1)) + + model.intervention.increment_tracing.contact_trace_household( + house_to=model.network.household(2), + house_from=model.network.household(1), + days_since_contact_occurred=0, + contact_trace_delay=0, + time=0 + ) + + model.simulate_one_step() + model.simulate_one_step() + + assert model.network.node(2).lfd_testing.being_lateral_flow_tested is True if __name__ == "__main__": diff --git a/examples/simple_test_script.py b/examples/simple_test_script.py index 4551003..9e054fd 100644 --- a/examples/simple_test_script.py +++ b/examples/simple_test_script.py @@ -34,7 +34,7 @@ # Create controller and add model, then run controller = BranchingProcessController(bpm.HouseholdLevelTracing(params)) -controller.run_simulation(10) +controller.run_simulation({"max_time": 10}) # Update parameters params['infection_reporting_prob'] = 0.5 @@ -42,7 +42,7 @@ # Re initialise with new parameters and Re-run controller.model = bpm.HouseholdLevelTracing(params) -controller.run_simulation(10) +controller.run_simulation({"max_time": 10}) # Add further parameters @@ -56,7 +56,7 @@ # Switch on a view (e.g. the timeline graph views) controller.timeline_view.set_display(True) controller.graph_view.set_display(True) -controller.run_simulation(10) +controller.run_simulation({"max_time": 10}) # Repeat runs and output to a named CSV file @@ -68,4 +68,4 @@ controller.csv_view.filename = save_path controller.csv_view.display_params = ['number_of_days_to_trace_backwards', 'number_of_days_to_trace_forwards'] - controller.run_simulation(5) + controller.run_simulation({"max_time": 5}) diff --git a/household_contact_tracing/behaviours/infection/contact_rate_reduction.py b/household_contact_tracing/behaviours/infection/contact_rate_reduction.py index aced721..9bd587e 100644 --- a/household_contact_tracing/behaviours/infection/contact_rate_reduction.py +++ b/household_contact_tracing/behaviours/infection/contact_rate_reduction.py @@ -58,9 +58,9 @@ def get_contact_rate_reduction(self, node: Node) -> int: intervention parameters """ - if node.isolated and node.propensity_imperfect_isolation: + if node.infection.isolated and node.tracing_adherence.propensity_imperfect_isolation: return self.global_contact_reduction_imperfect_quarantine - elif node.isolated and not node.propensity_imperfect_isolation: + elif node.infection.isolated and not node.tracing_adherence.propensity_imperfect_isolation: # return 1 means 100% of contacts are stopped return 1 else: @@ -81,15 +81,16 @@ def get_contact_rate_reduction(self, node: Node) -> int: """ # the isolated status should never apply to an individual who will not uptake intervention - if node.isolated and not node.propensity_imperfect_isolation: + if node.infection.isolated and not node.tracing_adherence.propensity_imperfect_isolation: # perfect intervention return 1 - elif node.isolated and node.propensity_imperfect_isolation: + elif node.infection.isolated and node.tracing_adherence.propensity_imperfect_isolation: # imperfect intervention return self.global_contact_reduction_imperfect_quarantine - elif node.being_lateral_flow_tested and node.propensity_risky_behaviour_lfa_testing: + elif (node.lfd_testing.being_lateral_flow_tested and + node.lfd_testing_adherence.propensity_risky_behaviour_lfa_testing): # engaging in risky behaviour while testing negative return self.global_contact_reduction_risky_behaviour diff --git a/household_contact_tracing/behaviours/infection/new_infection.py b/household_contact_tracing/behaviours/infection/new_infection.py index 9642556..ecb3b27 100644 --- a/household_contact_tracing/behaviours/infection/new_infection.py +++ b/household_contact_tracing/behaviours/infection/new_infection.py @@ -20,7 +20,7 @@ class NewInfection(ABC, Parameterised): Attributes ---------- - Todo: fill in descriptions of each attribute + # Todo: fill in descriptions of each attribute network: ContactTracingNetwork The store of Nodes and households used in the simulation symptom_reporting_delay @@ -42,10 +42,6 @@ class NewInfection(ABC, Parameterised): ------- new_infection(self, time: int, household: Household, infecting_node: Optional[Node] = None) - Add a new infected Node to the model. - :param time: The current simulation time. - :param household: The Household to create the new infection in. - :param infecting_node: The source of the new infection. """ @@ -137,12 +133,6 @@ def new_infection(self, time: int, household: Household, infecting_node: Optiona # Symptom onset time symptom_onset_time = time + self.incubation_period(asymptomatic) - # If the node is asymptomatic, we need to generate a pseudo symptom onset time - if asymptomatic: - pseudo_symptom_onset_time = self.incubation_period(asymptomatic=False) - else: - pseudo_symptom_onset_time = symptom_onset_time - # When a node reports its infection if not asymptomatic and np.random.binomial(1, self.infection_reporting_prob) == 1: will_report_infection = True @@ -151,10 +141,6 @@ def new_infection(self, time: int, household: Household, infecting_node: Optiona will_report_infection = False time_of_reporting = float('Inf') - # We assign each node a recovery period of 21 days, after 21 days the probability of - # causing a new infections is 0, due to the generation time distribution - recovery_time = time + 14 - # If the household has the propensity to use the contact tracing app, decide # if the node uses the app. if household.propensity_trace_app: @@ -162,28 +148,58 @@ def new_infection(self, time: int, household: Household, infecting_node: Optiona else: has_trace_app = False + tracing_attributes = { + 'contact_traced': household.contact_traced, + 'has_contact_tracing_app': has_trace_app, + 'symptom_onset_time': symptom_onset_time, + 'testing_delay': self.testing_delay(), + 'time_of_reporting': time_of_reporting, + 'will_report_infection': will_report_infection, + } + isolation_uptake = self.will_uptake_isolation() + tracing_adherence_attributes = {'will_uptake_isolation': isolation_uptake, + 'propensity_imperfect_isolation': self.get_propensity_imperfect_isolation() + } + + # If the node is asymptomatic, we need to generate a pseudo symptom onset time + if asymptomatic: + pseudo_symptom_onset_time = self.incubation_period(asymptomatic=False) + else: + pseudo_symptom_onset_time = symptom_onset_time + + returning_travellers_attributes = {'pseudo_symptom_onset_time': pseudo_symptom_onset_time} + + if infecting_node: + infecting_node_id = infecting_node.id + else: + infecting_node_id = None + if household.isolated and isolation_uptake: node_is_isolated = True else: node_is_isolated = False - new_node = self.network.add_node(time_infected=time, - household_id=household.id, isolated=node_is_isolated, - will_uptake_isolation=isolation_uptake, - propensity_imperfect_isolation=self.get_propensity_imperfect_isolation(), - asymptomatic=asymptomatic, contact_traced=household.contact_traced, - symptom_onset_time=symptom_onset_time, - pseudo_symptom_onset_time=pseudo_symptom_onset_time, - recovery_time=recovery_time, - will_report_infection=will_report_infection, - time_of_reporting=time_of_reporting, - has_contact_tracing_app=has_trace_app, - testing_delay=self.testing_delay(), - infecting_node=infecting_node) + # Each node has a recovery period after which the node recovers and cannot cause a new infection + recovery_time = time + 14 - # Each house now stores the ID's of which nodes are stored inside the house, + infection_attributes = { + 'time_infected': time, + 'asymptomatic': asymptomatic, + 'infecting_node': infecting_node_id, + 'isolated': node_is_isolated, + 'recovery_time': recovery_time, + } + + new_node = self.network.add_node(household_id=household.id, + infection_attributes=infection_attributes, + tracing_adherence_attributes=tracing_adherence_attributes, + returning_travellers_attributes=returning_travellers_attributes, + tracing_attributes=tracing_attributes + ) + + # Each house stores the IDs of which nodes are stored inside the house, # so that quarantining can be done at the household level household.nodes.append(new_node) @@ -217,23 +233,6 @@ def new_infection(self, time: int, household: Household, infecting_node: Optiona node_being_lateral_flow_tested = False time_started_lfa_testing = float('Inf') - additional_attributes = { - 'being_lateral_flow_tested': node_being_lateral_flow_tested, - 'time_started_lfa_testing': time_started_lfa_testing, - 'received_positive_test_result': False, - 'received_result': None, - 'avenue_of_testing': None, - 'positive_test_time': None, - 'node_will_take_up_lfa_testing': node_will_take_up_lfa_testing, - 'confirmatory_PCR_result_was_positive': None, - 'taken_confirmatory_PCR_test': False, - 'confirmatory_PCR_test_time': None, - 'confirmatory_PCR_test_result_time': None, - 'propensity_risky_behaviour_lfa_testing': - self.will_engage_in_risky_behaviour_while_being_lfa_tested(), - 'propensity_to_miss_lfa_tests': self.propensity_to_miss_lfa_tests() - } - asymptomatic = self.is_asymptomatic_infection() # Symptom onset time @@ -271,22 +270,62 @@ def new_infection(self, time: int, household: Household, infecting_node: Optiona else: node_is_isolated = False - new_node = self.network.add_node(time_infected=time, + lfd_testing_attributes = { + 'avenue_of_testing': None, + 'being_lateral_flow_tested': node_being_lateral_flow_tested, + 'time_started_lfa_testing': time_started_lfa_testing, + 'positive_test_time': None, + 'taken_confirmatory_PCR_test': False, + 'confirmatory_PCR_test_time': None, + 'confirmatory_PCR_test_result_time': None + } + + lfd_testing_adherence_attributes = { + 'node_will_take_up_lfa_testing': node_will_take_up_lfa_testing, + 'confirmatory_PCR_result_was_positive': None, + 'propensity_risky_behaviour_lfa_testing': self.will_engage_in_risky_behaviour_while_being_lfa_tested(), + 'propensity_to_miss_lfa_tests': self.propensity_to_miss_lfa_tests() + } + + tracing_attributes = { + 'received_positive_test_result': False, + 'received_result': None, + 'contact_traced': household.contact_traced, + 'has_contact_tracing_app': has_trace_app, + 'symptom_onset_time': symptom_onset_time, + 'testing_delay': self.testing_delay(), + 'time_of_reporting': time_of_reporting, + 'will_report_infection': will_report_infection, + } + + tracing_adherence_attributes = { + 'will_uptake_isolation': isolation_uptake, + 'propensity_imperfect_isolation': self.get_propensity_imperfect_isolation() + } + + returning_travellers_attributes = {'pseudo_symptom_onset_time': pseudo_symptom_onset_time} + + infecting_node_id = None + if infecting_node: + infecting_node_id = infecting_node.id + + infection_attributes = { + 'time_infected': time, + 'asymptomatic': asymptomatic, + 'infecting_node': infecting_node_id, + 'isolated': node_is_isolated, + 'recovery_time': recovery_time, + } + + new_node = self.network.add_node( household_id=household.id, - isolated=node_is_isolated, - will_uptake_isolation=isolation_uptake, - propensity_imperfect_isolation=self.get_propensity_imperfect_isolation(), - asymptomatic=asymptomatic, - contact_traced=household.contact_traced, - symptom_onset_time=symptom_onset_time, - pseudo_symptom_onset_time=pseudo_symptom_onset_time, - recovery_time=recovery_time, - will_report_infection=will_report_infection, - time_of_reporting=time_of_reporting, - has_contact_tracing_app=has_trace_app, - testing_delay=self.testing_delay(), - additional_attributes=additional_attributes, - infecting_node=infecting_node) + infection_attributes=infection_attributes, + tracing_adherence_attributes=tracing_adherence_attributes, + returning_travellers_attributes=returning_travellers_attributes, + lfd_testing_attributes=lfd_testing_attributes, + lfd_testing_adherence_attributes=lfd_testing_adherence_attributes, + tracing_attributes=tracing_attributes + ) # Each house now stores the ID's of which nodes are stored inside the house, # so that quarantining can be done at the household level diff --git a/household_contact_tracing/behaviours/intervention/increment_tracing.py b/household_contact_tracing/behaviours/intervention/increment_tracing.py index 0826466..143ebd6 100644 --- a/household_contact_tracing/behaviours/intervention/increment_tracing.py +++ b/household_contact_tracing/behaviours/intervention/increment_tracing.py @@ -66,25 +66,25 @@ def increment_contact_tracing(self, time: int): # Isolate all households under observation that now display symptoms (excludes those # who will not take up intervention if prob <1) for node in self.network.all_nodes(): - if node.symptom_onset_time <= time: - if node.contact_traced: - if not node.isolated: - if not node.completed_isolation: + if node.tracing.symptom_onset_time <= time: + if node.tracing.contact_traced: + if not node.infection.isolated: + if not node.tracing.completed_isolation: node.household.isolate_household(time) # Propagate the contact tracing for all households that self-reported and have had their # test results come back for node in self.network.all_nodes(): - if node.time_of_reporting + node.testing_delay == time: + if node.tracing.time_of_reporting + node.tracing.testing_delay == time: if not node.household.propagated_contact_tracing: self.propagate_contact_tracing(node.household, time) # Propagate the contact tracing for all households that are isolated due to exposure, # have developed symptoms and had a test come back for node in self.network.all_nodes(): - if node.symptom_onset_time <= time: + if node.tracing.symptom_onset_time <= time: if not node.household.propagated_contact_tracing: - if node.household.isolated_time + node.testing_delay <= time: + if node.household.isolated_time + node.tracing.testing_delay <= time: self.propagate_contact_tracing(node.household, time) # Update the contact tracing index of households @@ -171,9 +171,9 @@ def update_contact_tracing_index(self, time: int): # The testing delay must be passed # The testing delay starts when the house have been isolated and symptoms have # onset - critical_time = max(node.symptom_onset_time, household.isolated_time) + critical_time = max(node.tracing.symptom_onset_time, household.isolated_time) - if critical_time + node.testing_delay <= time: + if critical_time + node.tracing.testing_delay <= time: household.contact_tracing_index = 0 for index_1_hh in household.contact_traced_households: @@ -196,30 +196,30 @@ def prob_pcr_positive(self, fn: Callable[[int], float]): self._prob_pcr_positive = fn def pcr_test_node(self, node: Node, time: int): - node.received_result = True - infectious_age_when_tested = time - node.testing_delay - node.time_infected + node.tracing.received_result = True + infectious_age_when_tested = time - node.tracing.testing_delay - node.infection.time_infected prob_positive_result = self.prob_pcr_positive(infectious_age_when_tested) - node.avenue_of_testing = TestType.pcr + node.lfd_testing.avenue_of_testing = TestType.pcr if np.random.binomial(1, prob_positive_result) == 1: - node.received_positive_test_result = True - node.positive_test_time = time + node.tracing.received_positive_test_result = True + node.lfd_testing.positive_test_time = time else: - node.received_positive_test_result = False + node.tracing.received_positive_test_result = False def receive_pcr_test_results(self, time: int): # self reporting infections for node in self.network.all_nodes(): - if node.time_of_reporting + node.testing_delay == time: - if not node.contact_traced: - if not node.received_result: + if node.tracing.time_of_reporting + node.tracing.testing_delay == time: + if not node.tracing.contact_traced: + if not node.tracing.received_result: self.pcr_test_node(node, time) # contact traced nodes for node in self.network.all_nodes(): - if node.symptom_onset_time + node.testing_delay == time: - if node.contact_traced: - if not node.received_result: + if node.tracing.symptom_onset_time + node.tracing.testing_delay == time: + if node.tracing.contact_traced: + if not node.tracing.received_result: self.pcr_test_node(node, time) def increment_contact_tracing(self, time: int): @@ -243,15 +243,15 @@ def increment_contact_tracing(self, time: int): self.receive_pcr_test_results(time) for node in self.network.all_nodes(): - if node.symptom_onset_time <= time: - if node.received_positive_test_result: - if not node.isolated: - if not node.completed_isolation: + if node.tracing.symptom_onset_time <= time: + if node.tracing.received_positive_test_result: + if not node.infection.isolated: + if not node.tracing.completed_isolation: node.household.isolate_household(time) for node in self.network.all_nodes(): - if node.received_result: - if not node.propagated_contact_tracing: + if node.tracing.received_result: + if not node.tracing.propagated_contact_tracing: self.propagate_contact_tracing(node, time) def propagate_contact_tracing(self, node: Node, time: int): @@ -260,11 +260,13 @@ def propagate_contact_tracing(self, node: Node, time: int): when a household that is under surveillance develops symptoms + gets tested. """ # update the propagation data - node.propagated_contact_tracing = True + node.tracing.propagated_contact_tracing = True # Contact tracing attempted for the household that infected the household currently # propagating the infection - infected_by_node = node.infecting_node + infected_by_node = None + if node.infection.infecting_node_id: + infected_by_node = node.household.network.node(node.infection.infecting_node_id) # If the node was globally infected, we are backwards tracing and the infecting node is # not None @@ -272,22 +274,22 @@ def propagate_contact_tracing(self, node: Node, time: int): # if the infector is not already isolated and the time the node was infected captured # by going backwards - # the node.time_infected is when they had a contact with their infector. - if not infected_by_node.isolated and node.time_infected >= node.symptom_onset_time - \ + # the node.infection.time_infected is when they had a contact with their infector. + if not infected_by_node.infection.isolated and node.infection.time_infected >= node.tracing.symptom_onset_time - \ self.number_of_days_to_trace_backwards: # Then attempt to contact trace the household of the node that infected you - self.attempt_contact_trace_of_household( + self.contact_trace_household( house_to=infected_by_node.household, house_from=node.household, time=time, - days_since_contact_occurred=time - node.time_infected + days_since_contact_occurred=time - node.infection.time_infected ) # spread_to_global_node_time_tuples stores a list of tuples, where the first element is # the node_id of a node who was globally infected by the node, and the second element is # the time of transmission - for global_infection in node.spread_to_global_node_time_tuples: + for global_infection in node.infection.spread_to_global_node_time_tuples: # Get the child node_id and the time of transmission/time of contact child_node_id, time_t = global_infection @@ -295,23 +297,23 @@ def propagate_contact_tracing(self, node: Node, time: int): child_node = self.network.node(child_node_id) # If the node was infected 2 days prior to symptom onset, or 7 days post and is not already isolated - if time_t >= node.symptom_onset_time - self.number_of_days_to_trace_backwards and \ - time_t <= node.symptom_onset_time + self.number_of_days_to_trace_forwards and \ - not child_node.isolated: + if time_t >= node.tracing.symptom_onset_time - self.number_of_days_to_trace_backwards and \ + time_t <= node.tracing.symptom_onset_time + self.number_of_days_to_trace_forwards and \ + not child_node.infection.isolated: - self.attempt_contact_trace_of_household( + self.contact_trace_household( house_to=child_node.household, house_from=node.household, days_since_contact_occurred=time - time_t, time=time ) - def attempt_contact_trace_of_household(self, - house_to: Household, - house_from: Household, - days_since_contact_occurred: int, - time: int, - contact_trace_delay: int = 0): + def contact_trace_household(self, + house_to: Household, + house_from: Household, + days_since_contact_occurred: int, + time: int, + contact_trace_delay: int = 0): # Decide if the edge was traced by the app app_traced = self.network.is_edge_app_traced(self.network.get_edge_between_household(house_from, house_to)) @@ -358,32 +360,32 @@ def receive_pcr_test_results(self, time: int): super().receive_pcr_test_results(time) else: for node in self.network.all_nodes(): - if node.time_of_reporting + node.testing_delay == time: - if not node.contact_traced: - if not node.received_result: - if not node.being_lateral_flow_tested: + if node.tracing.time_of_reporting + node.tracing.testing_delay == time: + if not node.tracing.contact_traced: + if not node.tracing.received_result: + if not node.lfd_testing.being_lateral_flow_tested: self.pcr_test_node(node, time) def increment_contact_tracing(self, time: int): for node in self.network.all_nodes(): - if node.received_positive_test_result: - if node.avenue_of_testing == TestType.pcr: - if not node.propagated_contact_tracing: + if node.tracing.received_positive_test_result: + if node.lfd_testing.avenue_of_testing == TestType.pcr: + if not node.tracing.propagated_contact_tracing: self.propagate_contact_tracing(node, time) if not self.LFA_testing_requires_confirmatory_PCR: for node in self.network.all_nodes(): - if node.received_positive_test_result: - if node.avenue_of_testing == TestType.lfa: - if not node.propagated_contact_tracing: + if node.tracing.received_positive_test_result: + if node.lfd_testing.avenue_of_testing == TestType.lfa: + if not node.tracing.propagated_contact_tracing: self.propagate_contact_tracing(node, time) elif self.LFA_testing_requires_confirmatory_PCR: for node in self.network.all_nodes(): - if node.confirmatory_PCR_test_result_time == time: - if node.confirmatory_PCR_result_was_positive: - if node.avenue_of_testing == TestType.lfa: - if not node.propagated_contact_tracing: + if node.lfd_testing.confirmatory_PCR_test_result_time == time: + if node.lfd_testing_adherence.confirmatory_PCR_result_was_positive: + if node.lfd_testing.avenue_of_testing == TestType.lfa: + if not node.tracing.propagated_contact_tracing: self.propagate_contact_tracing(node, time) def propagate_contact_tracing(self, node: Node, time: int): @@ -397,81 +399,84 @@ def propagate_contact_tracing(self, node: Node, time: int): # 2) Trace on confirmatory PCR result # update the propagation data - node.propagated_contact_tracing = True + node.tracing.propagated_contact_tracing = True # Contact tracing attempted for the household that infected the household currently # propagating the infection - infected_by_node = node.infecting_node + + infected_by_node = None + if node.infection.infecting_node_id: + infected_by_node = node.household.network.node(node.infection.infecting_node_id) # If the node was globally infected, we are backwards tracing and the infecting node # is not None if not node.locally_infected() and infected_by_node: # if the infector is not already isolated and the time the node was infected captured - # by going backwards the node.time_infected is when they had a contact with their + # by going backwards the node.infection.time_infected is when they had a contact with their # infector. - if node.avenue_of_testing == TestType.pcr: + if node.lfd_testing.avenue_of_testing == TestType.pcr: - if not infected_by_node.isolated and \ - node.time_infected >= node.symptom_onset_time - \ - self.number_of_days_to_trace_backwards: + if not infected_by_node.infection.isolated and \ + node.infection.time_infected >= node.tracing.symptom_onset_time - \ + self.number_of_days_to_trace_backwards: # Then attempt to contact trace the household of the node that infected you - self.attempt_contact_trace_of_household( + self.contact_trace_household( house_to=infected_by_node.household, house_from=node.household, - days_since_contact_occurred=time - node.time_infected, + days_since_contact_occurred=time - node.infection.time_infected, time=time) - elif node.avenue_of_testing == TestType.lfa: + elif node.lfd_testing.avenue_of_testing == TestType.lfa: if not self.LFA_testing_requires_confirmatory_PCR: - if not infected_by_node.isolated and node.time_infected >= \ - node.positive_test_time - self.number_of_days_prior_to_LFA_result_to_trace: + if not infected_by_node.tracing.isolated and node.infection.time_infected >= \ + node.lfd_testing.positive_test_time - self.number_of_days_prior_to_LFA_result_to_trace: # Then attempt to contact trace the household of the node that infected you - self.attempt_contact_trace_of_household( + self.contact_trace_household( house_to=infected_by_node.household, house_from=node.household, - days_since_contact_occurred=time - node.time_infected, - time=time ) + days_since_contact_occurred=time - node.infection.time_infected, + time=time) # spread_to_global_node_time_tuples stores a list of tuples, where the first element is # the node_id of a node who was globally infected by the node, and the second element is # the time of transmission - for global_infection in node.spread_to_global_node_time_tuples: + for global_infection in node.infection.spread_to_global_node_time_tuples: # Get the child node_id and the time of transmission/time of contact child_node_id, time_t = global_infection child_node = self.network.node(child_node_id) - if node.avenue_of_testing == TestType.pcr: + if node.lfd_testing.avenue_of_testing == TestType.pcr: # If the node was infected 2 days prior to symptom onset, or 7 days post and is # not already isolated - if time_t >= node.symptom_onset_time - self.number_of_days_to_trace_backwards: - if time_t <= node.symptom_onset_time + self.number_of_days_to_trace_forwards: - if not child_node.isolated: + if time_t >= node.tracing.symptom_onset_time - self.number_of_days_to_trace_backwards: + if time_t <= node.tracing.symptom_onset_time + self.number_of_days_to_trace_forwards: + if not child_node.infection.isolated: - self.attempt_contact_trace_of_household( - house_to=child_node.household, - house_from=node.household, - days_since_contact_occurred=time - time_t, - time=time) + self.contact_trace_household( + house_to=child_node.household, + house_from=node.household, + days_since_contact_occurred=time - time_t, + time=time) - elif node.avenue_of_testing == TestType.lfa: + elif node.lfd_testing.avenue_of_testing == TestType.lfa: if not self.LFA_testing_requires_confirmatory_PCR: # If the node was infected 2 days prior to symptom onset, or 7 days post and # is not already isolated - if time_t >= node.positive_test_time - \ + if time_t >= node.lfd_testing.positive_test_time - \ self.number_of_days_prior_to_LFA_result_to_trace: - self.attempt_contact_trace_of_household( + self.contact_trace_household( house_to=child_node.household, house_from=node.household, days_since_contact_occurred=time - time_t, - time=time) \ No newline at end of file + time=time) diff --git a/household_contact_tracing/behaviours/intervention/isolation.py b/household_contact_tracing/behaviours/intervention/isolation.py index 59fb389..24cf62e 100644 --- a/household_contact_tracing/behaviours/intervention/isolation.py +++ b/household_contact_tracing/behaviours/intervention/isolation.py @@ -89,9 +89,9 @@ def isolate_self_reporting_cases(self, time: int): who will not uptake intervention """ for node in self.network.all_nodes(): - if node.will_uptake_isolation: - if node.time_of_reporting == time: - node.isolated = True + if node.tracing_adherence.will_uptake_isolation: + if node.tracing.time_of_reporting == time: + node.infection.isolated = True def update_households_contact_traced(self, time: int): """Update the contact traced status for all households that have had the @@ -105,7 +105,7 @@ def update_households_contact_traced(self, time: int): def update_isolation(self, time: int): for node in self.network.all_nodes(): - if node.time_of_reporting + node.testing_delay == time: + if node.tracing.time_of_reporting + node.tracing.testing_delay == time: if not node.household.isolated: if not node.household.contact_traced: node.household.isolate_household(time) @@ -123,8 +123,8 @@ def update_households_contact_traced(self, time: int): def update_isolation(self, time: int): for node in self.network.all_nodes(): - if node.time_of_reporting + node.testing_delay == time: - if node.received_positive_test_result: + if node.tracing.time_of_reporting + node.tracing.testing_delay == time: + if node.tracing.received_positive_test_result: if not node.household.isolated: if not node.household.contact_traced: node.household.isolate_household(time) @@ -153,16 +153,16 @@ def update_households_contact_traced(self, time: int): household.update_network() traced_node = household.find_traced_node() # the traced node is now being lateral flow tested - if traced_node.node_will_take_up_lfa_testing: - if not traced_node.received_positive_test_result: - traced_node.being_lateral_flow_tested = True - traced_node.time_started_lfa_testing = time + if traced_node.lfd_testing_adherence.node_will_take_up_lfa_testing: + if not traced_node.tracing.received_positive_test_result: + traced_node.lfd_testing.being_lateral_flow_tested = True + traced_node.lfd_testing.time_started_lfa_testing = time def update_isolation(self, time: int): for node in self.network.all_nodes(): - if node.positive_test_time == time: - if node.avenue_of_testing == TestType.pcr: - if node.received_positive_test_result: + if node.lfd_testing.positive_test_time == time: + if node.lfd_testing.avenue_of_testing == TestType.pcr: + if node.tracing.received_positive_test_result: if not node.household.applied_household_positive_policy: node.household.apply_positive_policy(time, self.household_positive_policy) @@ -174,7 +174,7 @@ def act_on_confirmatory_pcr_results(self, time: int): * Contact tracing is propagated """ for node in self.network.all_nodes(): - if node.confirmatory_PCR_test_result_time == time: + if node.lfd_testing.confirmatory_PCR_test_result_time == time: node.household.apply_positive_policy(time, self.household_positive_policy) def isolate_positive_lateral_flow_tests(self, time: int, positive_nodes: List[Node]): @@ -185,30 +185,30 @@ def isolate_positive_lateral_flow_tests(self, time: int, positive_nodes: List[No """ for node in positive_nodes: - node.received_positive_test_result = True + node.tracing.received_positive_test_result = True - if node.will_uptake_isolation: - node.isolated = True + if node.tracing_adherence.will_uptake_isolation: + node.infection.isolated = True - node.avenue_of_testing = TestType.lfa - node.positive_test_time = time - node.being_lateral_flow_tested = False + node.lfd_testing.avenue_of_testing = TestType.lfa + node.lfd_testing.positive_test_time = time + node.lfd_testing.being_lateral_flow_tested = False if not node.household.applied_household_positive_policy and \ not self.LFA_testing_requires_confirmatory_PCR: node.household.apply_positive_policy(time, self.household_positive_policy) - def act_on_positive_LFA_tests(self, time: int, positive_nodes: List[Node]): + def act_on_positive_lfa_tests(self, time: int, positive_nodes: List[Node]): """For nodes who test positive on their LFA test, take the appropriate action depending on the policy """ self.isolate_positive_lateral_flow_tests(time, positive_nodes) if self.LFA_testing_requires_confirmatory_PCR: - self.confirmatory_pcr_test_LFA_nodes(time, positive_nodes) + self.confirmatory_pcr_test_lfa_nodes(time, positive_nodes) - def confirmatory_pcr_test_LFA_nodes(self, time: int, positive_nodes: List[Node]): + def confirmatory_pcr_test_lfa_nodes(self, time: int, positive_nodes: List[Node]): """Nodes who receive a positive LFA result will be tested using a PCR test.""" for node in positive_nodes: - if not node.taken_confirmatory_PCR_test: + if not node.lfd_testing.taken_confirmatory_PCR_test: node.take_confirmatory_pcr_test(time, self.prob_pcr_positive) diff --git a/household_contact_tracing/branching_process_controller.py b/household_contact_tracing/branching_process_controller.py index d301ee0..11b5c8a 100644 --- a/household_contact_tracing/branching_process_controller.py +++ b/household_contact_tracing/branching_process_controller.py @@ -1,9 +1,13 @@ +from typing import List, Optional +from household_contact_tracing.views.branching_process_view import BranchingProcessView +from household_contact_tracing.views.statistics_view import StatisticsView from household_contact_tracing.branching_process_model import BranchingProcessModel from household_contact_tracing.views.shell_view import ShellView from household_contact_tracing.views.csv_file_view import CSVFileView from household_contact_tracing.views.graph_view import GraphView from household_contact_tracing.views.graph_pyvis_view import GraphPyvisView from household_contact_tracing.views.timeline_graph_view import TimelineGraphView +from household_contact_tracing.views.statistics_view import StatisticsView class BranchingProcessController: @@ -26,12 +30,12 @@ class BranchingProcessController: set_graphic_displays(self, display: bool) choose whether to show the graphical outputs - run_simulation(self, max_time: int = 20, infection_threshold: int = 5000) + run_simulation(self, max_time: int = 20, max_active_infections: int = 5000) runs the simulation """ - def __init__(self, model: BranchingProcessModel): + def __init__(self, model: BranchingProcessModel, additional_views: Optional[List[BranchingProcessView]] = []): """ Constructor for BranchingProcessController @@ -47,6 +51,13 @@ def __init__(self, model: BranchingProcessModel): self.timeline_view = TimelineGraphView(model) self.shell_view = ShellView(model) self.csv_view = CSVFileView(model) + self.statistics_view = StatisticsView(model) + + # initialise any views that are required, but included as defaults + for view in additional_views: + + initialised_view = view(model) + setattr(self, initialised_view.view_name, initialised_view) self.set_graphic_displays(False) @@ -79,15 +90,32 @@ def set_graphic_displays(self, display: bool): self.graph_pyvis_view.set_display(display) self.timeline_view.set_display(display) - def run_simulation(self, max_time: int = 20, infection_threshold: int = 5000): + def run_simulation(self, state_criteria: dict): """ Run the simulation until it stops (e.g times out, too many infectious nodes or goes extinct) Parameters: - max_time (int): The maximum number of iterations (eg. days) to be run (simulation stops if reached) - infection_threshold (int): The maximum number of infectious nodes (simulation stops if reached) + state_criteria: Named variables which are evaluated each step of the model to determine + whether the state of the model will change. + + Returns: + None + """ + + self._model.run_simulation(state_criteria) + + def run_hh_sar_simulation(self, state_criteria: dict = {}): + """ + This simulation method with only simulate the infection process for households in the first + generation of the epidemic, and will continue until all nodes in the initial households of the + epidemic are recovered. This is primarily useful when we are estimating the household secondary + attack rate. If we simulated onwards transmission, and examined households where the local + epidemic was completed, we would end up with a biased sample - the longer local epidemics would + be less likely to be included in the sample. Returns: None """ - self._model.run_simulation(max_time, infection_threshold) + self._model.run_hh_sar_simulation(state_criteria = state_criteria) + + diff --git a/household_contact_tracing/branching_process_model.py b/household_contact_tracing/branching_process_model.py index e3b1f24..eceeaa5 100644 --- a/household_contact_tracing/branching_process_model.py +++ b/household_contact_tracing/branching_process_model.py @@ -20,7 +20,7 @@ class BranchingProcessModel(ABC, Parameterised): Methods ------- - run_simulation(self, max_time: int, infection_threshold: int) -> None: + run_simulation(self, max_time: int, max_active_infections: int) -> None: runs the simulation """ @@ -42,6 +42,7 @@ def __init__(self): # Set state self._state = ReadyState(self) + self.state_criteria = [] @property def state(self) -> BranchingProcessState: @@ -74,13 +75,31 @@ def root_dir(self) -> str: return self.__ROOT_DIR @abstractmethod - def run_simulation(self, max_time: int, infection_threshold: int) -> None: + def run_simulation(self, state_criteria: dict) -> None: """ Run the simulation until it stops (e.g times out, too many infectious nodes or goes extinct) + Parameters: + state_criteria: Named variables which are evaluated each step of the model to determine + whether the state of the model will change. + + Returns: + None + """ + + @abstractmethod + def run_hh_sar_simulation(self) -> None: + """ + This simulation method with only simulate the infection process for households in the first + generation of the epidemic, and will continue until all nodes in the initial households of the + epidemic are recovered. This is primarily useful when we are estimating the household secondary + attack rate. If we simulated onwards transmission, and examined households where the local + epidemic was completed, we would end up with a biased sample - the longer local epidemics would + be less likely to be included in the sample. + Parameters: max_time (int): The maximum number of iterations (eg. days) to be run (simulation stops if reached) - infection_threshold (int): The maximum number of infectious nodes (simulation stops if reached) + max_active_infections (int): The maximum number of infectious nodes (simulation stops if reached) Returns: None diff --git a/household_contact_tracing/branching_process_models.py b/household_contact_tracing/branching_process_models.py index c38db6d..8462616 100644 --- a/household_contact_tracing/branching_process_models.py +++ b/household_contact_tracing/branching_process_models.py @@ -1,3 +1,4 @@ +import math import os from typing import Callable from copy import deepcopy @@ -14,6 +15,7 @@ import household_contact_tracing.behaviours.intervention.increment_tracing as increment import household_contact_tracing.behaviours.intervention.isolation as isolation import household_contact_tracing.behaviours.infection.new_infection as new_infection +from household_contact_tracing.utilities import ParameterError class HouseholdLevelTracing(BranchingProcessModel): @@ -32,7 +34,7 @@ class HouseholdLevelTracing(BranchingProcessModel): Methods ------- - run_simulation(self, max_time: int, infection_threshold: int = 1000) -> None + run_simulation(self, max_time: int, max_active_infections: int = 1000) -> None Runs the simulation up to a maximum number of increments and max allowed number of infected nodes. @@ -78,6 +80,44 @@ def _initialise_intervention(self): increment.IncrementTracingHouseholdLevel, self.params) + def evaluate_model_state(self, ): + """Determine whether the state of the model has changed by evaluating the data from the last simulation step + against criteria which trigger a change of state.""" + + if self.time >= self.state_criteria["max_time"]: + # Simulation ends if max_time is reached + self.state.switch(TimedOutState, {"total_increments": self.time, + "non_recovered_nodes": self.network.count_non_recovered_nodes(), + "total_nodes": self.network.node_count}) + elif self.network.count_non_recovered_nodes() == self.state_criteria["min_non_recovered_nodes"]: + # Simulation ends if no more infectious nodes + self.state.switch(ExtinctState, {"total_increments": self.time, + "non_recovered_nodes": self.network.count_non_recovered_nodes(), + "total_nodes": self.network.node_count}) + elif self.network.count_non_recovered_nodes() > self.state_criteria["infection_threshold"]: + # Simulation ends if number of infectious nodes > threshold + self.state.switch(MaxNodesInfectiousState, {"total_increments": self.time, + "non_recovered_nodes": 0, + "total_nodes": self.network.node_count}) + + def set_default_state_criteria(self): + """Set default values for the state criteria if they have not yet been set.""" + valid_state_criteria = ["max_time", "min_non_recovered_nodes", "infection_threshold"] + + for criterion in self.state_criteria: + if criterion not in valid_state_criteria: + raise ParameterError(f"Criterion '{criterion}', is not a valid state criterion.\n" + f"Valid state criteria are: {valid_state_criteria}.") + + if "infection_threshold" not in self.state_criteria: + self.state_criteria["infection_threshold"] = 10000 + + if "max_time" not in self.state_criteria: + self.state_criteria["max_time"] = 40 + + if "min_non_recovered_nodes" not in self.state_criteria: + self.state_criteria["min_non_recovered_nodes"] = 0 + def simulate_one_step(self): """Simulates one day of the infection and contact tracing.""" @@ -100,23 +140,24 @@ def simulate_one_step(self): # increment time self.time += 1 - def run_simulation(self, max_time: int, infection_threshold: int = 1000) -> None: + def run_simulation(self, state_criteria: dict) -> None: """ Runs the simulation: Sets model state, Announces start/stopped and step increments to observers Arguments: - max_time -- The maximum number of step increments to perform (stops if self.time >= - max_time). Self.time is cumulative throughout multiple calls to run_simulation. - infection_threshold -- The maximum number of infectious nodes allowed, - before stopping simulation + state_criteria: Named variables which are evaluated each step of the model to determine + whether the state of the model will change. Returns: None """ + self.state_criteria = state_criteria + + self.set_default_state_criteria() # Switch model to RunningState - self._state.switch(RunningState, max_time=max_time, infection_threshold=infection_threshold) + self._state.switch(RunningState, self.state_criteria) while type(self.state) is RunningState: prev_network = deepcopy(self.network) @@ -131,29 +172,60 @@ def run_simulation(self, max_time: int, infection_threshold: int = 1000) -> None # Call parent completed step super()._completed_step_increment() - if self.time >= max_time: - # Simulation ends if max_time is reached - self.state.switch(TimedOutState, - total_increments=self.time, - non_recovered_nodes=self.network.count_non_recovered_nodes(), - total_nodes=self.network.node_count - ) - elif self.network.count_non_recovered_nodes() == 0: - # Simulation ends if no more infectious nodes - self.state.switch(ExtinctState, - total_increments=self.time, - non_recovered_nodes=0, - total_nodes=self.network.node_count) - elif self.network.count_non_recovered_nodes() > infection_threshold: - # Simulation ends if number of infectious nodes > threshold - self.state.switch(MaxNodesInfectiousState, - total_increments=self.time, - non_recovered_nodes=0, - total_nodes=self.network.node_count) + self.evaluate_model_state() # Tell parent simulation stopped super()._simulation_stopped() + def run_hh_sar_simulation(self, state_criteria: dict) -> None: + """ Runs the simulation only for the first generation of the household epidemic: + Sets model state, + Announces start/stopped and step increments to observers + + This simulation method with only simulate the infection process for households in the first + generation of the epidemic, and will continue until all nodes in the initial households of the + epidemic are recovered. This is primarily useful when we are estimating the household secondary + attack rate. If we simulated onwards transmission, and examined households where the local + epidemic was completed, we would end up with a biased sample - the longer local epidemics would + be less likely to be included in the sample. + + Returns: + None + """ + self.state_criteria = state_criteria + + self.set_default_state_criteria() + + # Switch model to RunningState + self._state.switch(RunningState, self.state_criteria) + + while type(self.state) is RunningState: + #prev_network = deepcopy(self.network) + + # This chunk of code executes a days worth of infections and recoveries, but no tracing + self.infection.increment(self.time) + self.infection.perform_recoveries(self.time) + + # if an infection is in a second generation household, set them to recovered so that + # they do not infect. This is mainly for computational ease + for node in self.network.all_nodes(): + if node.household.id not in self.infection.starting_households: + node.recovered = True + + self.time += 1 + + # If graph changed, tell parent + #if not prev_network == self.network: + # BranchingProcessModel.graph_changed(self) + + # Call parent completed step + super()._completed_step_increment() + + # the simulation ends when all nodes in the initial generation have recovered + self.evaluate_model_state() + + # Tell parent simulation stopped + super()._simulation_stopped() class IndividualLevelTracing(HouseholdLevelTracing): """ @@ -189,16 +261,16 @@ def prob_pcr_positive(self, fn: Callable[[int], float]): @staticmethod def default_prob_lfa_positive(infectious_age): """Default LFA test result probability.""" - if infectious_age in [4, 5, 6]: - return 1 + if infectious_age in [2, 3, 4, 5, 6, 7]: + return 0.75 else: return 0 @staticmethod def default_prob_pcr_positive(infectious_age): """Default PCR test result probability.""" - if infectious_age in [4, 5, 6]: - return 0 + if infectious_age in [1, 2, 3, 4, 5, 6, 7, 8]: + return 1 else: return 0 @@ -226,27 +298,6 @@ def _initialise_intervention(self): return new_intervention - def simulate_one_step(self): - """Simulates one day of the infection and contact tracing.""" - - # Perform one day of the infection - self.infection.increment(self.time) - # isolate nodes reached by tracing, isolate nodes due to self-reporting - self.intervention.isolation.isolate_self_reporting_cases(self.time) - # isolate self-reporting-nodes while they wait for tests - self.intervention.isolation.update_households_contact_traced(self.time) - self.intervention.isolation.update_isolation(self.time) - for step in range(5): - self.intervention.increment_tracing.increment_contact_tracing(self.time) - # node recoveries - self.infection.perform_recoveries(self.time) - # release nodes from quarantine or intervention if the time has arrived - self.intervention.completed_isolation(self.time) - self.intervention.completed_quarantine(self.time) - # increment time - self.time += 1 - - class IndividualTracingDailyTesting(IndividualLevelTracing): """A class used to represent a simulation of contact tracing of households along with contacting every individual and their contacts, whether they have tested positive or not, along @@ -271,10 +322,10 @@ def _initialise_infection(self): def _initialise_intervention(self): """ Initialise an Intervention class, passing in the required behaviours into its constructor """ - new_intervention = Intervention(self.network, - isolation.DailyTestingIsolation, - increment.IncrementTracingIndividualDailyTesting, - self.params) + new_intervention = Intervention(self.network, + isolation.DailyTestingIsolation, + increment.IncrementTracingIndividualDailyTesting, + self.params) # Set a new positive pcr probability function new_intervention.increment_tracing.prob_pcr_positive = self.prob_pcr_positive @@ -293,7 +344,7 @@ def simulate_one_step(self): self.intervention.isolation.update_isolation(self.time) # isolate self reporting nodes positive_nodes = self.intervention.lft_nodes(self.time, self.prob_lfa_positive) - self.intervention.isolation.act_on_positive_LFA_tests(self.time, positive_nodes) + self.intervention.isolation.act_on_positive_lfa_tests(self.time, positive_nodes) # if we require PCR tests, to confirm infection we act on those if self.intervention.increment_tracing.LFA_testing_requires_confirmatory_PCR: self.intervention.increment_tracing.act_on_confirmatory_pcr_results(self.time) diff --git a/household_contact_tracing/branching_process_state.py b/household_contact_tracing/branching_process_state.py index da0bf5b..9d06a60 100644 --- a/household_contact_tracing/branching_process_state.py +++ b/household_contact_tracing/branching_process_state.py @@ -30,13 +30,14 @@ class BranchingProcessState(ABC): def __init__(self, simulation_model): self._simulation_model = simulation_model - def switch(self, state: Type[BranchingProcessState], **state_info): + def switch(self, state: Type[BranchingProcessState], state_info: dict): """ Switch to a new state if new state is allowed (for current state) (If not, raise ValueError) Parameters: state: The new state to be switched to + state_info: A dictionary of information about the model at the time of the change of state. Returns: None diff --git a/household_contact_tracing/calibration.py b/household_contact_tracing/calibration.py new file mode 100644 index 0000000..5bd1b33 --- /dev/null +++ b/household_contact_tracing/calibration.py @@ -0,0 +1,221 @@ +""" +Code that deals with common hyperparameter optimisation routines. + +Often, we want to calibrate an epidemic to a desired growth rate and household secondary attack rate. +""" + +from abc import ABC +from household_contact_tracing.branching_process_models import HouseholdLevelTracing +from household_contact_tracing.branching_process_controller import BranchingProcessController +from ax import optimize +from copy import Error, copy +import math + +class Calibration(ABC): + """ + Base class for hyperparameter optimisation of infection dynamics. + """ + + def __init__(self) -> None: + pass + + def setup_and_run_model(self): + pass + + def compute_evaluation(self) -> float: + pass + + def optimize(self): + pass + + def plot_results(self): + pass + + def evaluate_fit(self): + pass + +class StandardCalibrationHouseholdLevelTracing(Calibration): + """Our standard calibration of HouseholdLevelTracing tunes the models growth rate and + household secondary attack rate given inputs: asymptomatic prob, asymptomatic relative infectiousness, + symptom reporting probability. + + The calibration is carried out in the absence of contact tracing to define a baseline epidemic. Contact tracing + can then be evaluated against the baseline epidemic. + """ + + def __init__( + self, + #household_pairwise_survival_prob: float, + desired_growth_rate: float, + desired_hh_sar: float, + asymptomatic_prob: float, + asymptomatic_relative_infectivity: float, + infection_reporting_prob: float, + reduce_contacts_by: float, + starting_infections: int = 100, + starting_infections_hh_sar: int = 1000 + ): + + # initialise non-infection parameters that are held constant between simulations + self.fixed_params = { + 'household_pairwise_survival_prob': 0.2, + 'contact_tracing_success_prob': 0.0, + 'overdispersion': 0.32, + 'infection_reporting_prob': 0.25, + 'contact_trace': False, + 'test_delay': 2, + 'contact_trace_delay': 1, + 'incubation_period_delay': 5, + 'symptom_reporting_delay': 1, + 'do_2_step': False, + 'reduce_contacts_by': 0.6, + 'prob_has_trace_app': 0, + 'hh_propensity_to_use_trace_app': 1, + 'test_before_propagate_tracing': True, + 'starting_infections': 100, + 'node_will_uptake_isolation_prob': 1, + 'self_isolation_duration': 0, + 'quarantine_duration': 0, + 'transmission_probability_multiplier': 1, + 'propensity_imperfect_quarantine': 0, + 'global_contact_reduction_imperfect_quarantine': 0 + } + + # set the inputted defaults + self.fixed_params['asymptomatic_prob'] = asymptomatic_prob + self.fixed_params['asymptomatic_relative_infectivity'] = asymptomatic_relative_infectivity + self.fixed_params['infection_reporting_prob'] = infection_reporting_prob + self.fixed_params['reduce_contacts_by'] = reduce_contacts_by + self.fixed_params['starting_infections'] = starting_infections + + self.desired_growth_rate = desired_growth_rate + self.desired_hh_sar = desired_hh_sar + self.starting_infections_hh_sar = starting_infections_hh_sar + + self.optimisation_complete = False + + def eval_metrics( + self, + household_pairwise_survival_prob: float, + outside_household_infectivity_scaling: float, + state_criteria: dict = {}, + verbose: bool = True) -> float: + """Sets up a model, runs it, and returns the evaluated growth rate. + + Args: + outside_household_infectivity_scaling (float): controls how infectious global contacts are + max_time (int, optional): upper limit of days to simulate. Defaults to 20. + max_active_infections (int, optional): simulation ends early if maximum number of infections is exceeded. Defaults to 1e5. + """ + + params = copy(self.fixed_params) + params['outside_household_infectivity_scaling'] = outside_household_infectivity_scaling + params['household_pairwise_survival_prob'] = household_pairwise_survival_prob + + # run a simulation to get the growth rate of the epidemic + controller = BranchingProcessController(HouseholdLevelTracing(params)) + controller.csv_view.set_display(False) + controller.run_simulation(state_criteria) + + # use a different simulation method to get the household secondary attack rate of the epidemic + params['starting_infections'] = self.starting_infections_hh_sar # use a higher number of starting infections + controller_hh_sar = BranchingProcessController(HouseholdLevelTracing(params)) + controller_hh_sar.csv_view.set_display(False) + controller_hh_sar.run_hh_sar_simulation( + state_criteria = { + 'infection_threshold': math.inf, + 'max_time': math.inf + } + ) + + return { + 'growth_rate': controller.statistics_view.get_growth_rate(verbose = verbose), + 'hh_sar': controller_hh_sar.statistics_view.get_hh_sar() + } + + + def evaluate_fit( + self, + household_pairwise_survival_prob, + outside_household_infectivity_scaling, + verbose: bool = True, + state_criteria: dict = {}) -> float: + + metrics = self.eval_metrics( + household_pairwise_survival_prob, + outside_household_infectivity_scaling, + state_criteria, + verbose + ) + + return abs(self.desired_growth_rate - metrics['growth_rate']) + abs(self.desired_hh_sar - metrics['hh_sar']) + + def optimise(self, + outside_household_infectivity_scaling_range: list[float], + household_pairwise_survival_prob_range: list[float], + total_trials: int = 20, + state_criteria: dict = {}): + """Performs the hyperparameter optimization step with proposals from the specified ranges. + + Args: + outside_household_infectivity_scaling_range (list[float]): The lower and upper values for this parameter. + total_trials (int): The total number of trials to perform. Defaults to 20 + """ + + self.best_parameters, self.values, self.experiment, self.model = optimize( + parameters=[ + { + "name": "outside_household_infectivity_scaling", + "type": "range", + "bounds": outside_household_infectivity_scaling_range, + "value_type": "float" + }, + { + "name": "household_pairwise_survival_prob", + "type": "range", + "bounds": household_pairwise_survival_prob_range, + "value_type": "float" + } + ], + evaluation_function = lambda pars: self.evaluate_fit( + household_pairwise_survival_prob = pars["household_pairwise_survival_prob"], + outside_household_infectivity_scaling = pars["outside_household_infectivity_scaling"], + verbose = False, + state_criteria = state_criteria + ), + minimize = True, + total_trials = total_trials + ) + + self.optimisation_complete = True + + return self.best_parameters, self.values + + def get_fitted_model_metric_samples( + self, + n_obs: int = 10, + state_criteria: dict = {}) -> list[dict]: + """If optimisation has been completed, this method generates sample of the growth rate using + the results from the optimisation step. + + + Args: + n_obs (int, optional): Number of fitted sample of the growth rate to get. Defaults to 10. + + Returns: + [list]: A list containing fitted samples of the growth rate. + """ + + if self.optimisation_complete: + + return [ + self.eval_metrics( + household_pairwise_survival_prob = self.best_parameters['household_pairwise_survival_prob'], + outside_household_infectivity_scaling = self.best_parameters['outside_household_infectivity_scaling'], + state_criteria = state_criteria, + verbose = False) + for _ in range(n_obs) + ] + + else: + raise Error('Optimisation has not yet been performed. Please run optimise before trying to get fitted samples.') diff --git a/household_contact_tracing/exceptions.py b/household_contact_tracing/exceptions.py new file mode 100644 index 0000000..775b840 --- /dev/null +++ b/household_contact_tracing/exceptions.py @@ -0,0 +1,22 @@ +# error handling module + +from household_contact_tracing.branching_process_state import BranchingProcessState + + +class Error(Exception): + """Base class for exceptions in this module""" + pass + + +class ModelStateError(Error): + """Exception raised when model is in an inappropriate state + when a method or function is called. + + Args: + state (BranchingProcessState): the state of the model when the error occurred. + message (str): explanation of why the state was incorrect + """ + + def __init__(self, state: BranchingProcessState, message: str): + self.state = state + self.message = message diff --git a/household_contact_tracing/infection.py b/household_contact_tracing/infection.py index 5090c92..7800dc3 100644 --- a/household_contact_tracing/infection.py +++ b/household_contact_tracing/infection.py @@ -51,8 +51,6 @@ def __init__(self, self.household_pairwise_survival_prob = 0.2 self.update_params(params) - # Perform initial set-up calculations - # Precomputing the cdf's for generating the overdispersed contact data household_size = len(self.total_contact_means) self.cdf_dict = {i + 1: compute_negbin_cdf(self.total_contact_means[i], self.overdispersion) @@ -80,6 +78,8 @@ def __init__(self, self.new_infection = new_infection(self.network, params) self.contact_rate_reduction = contact_rate_reduction(params) + self.starting_households = [] + # Initialise starting infections of simulation self.initialise() @@ -89,13 +89,16 @@ def initialise(self): new_household = self.new_household.new_household(0, None) self.new_infection.new_infection(0, new_household) + # store the set of starting households + self.starting_households = [household.id for household in self.network.all_households] + def increment(self, time): """Create a new days worth of infections.""" for node in self.network.active_infections: household = node.household # Extracting useful parameters from the node - days_since_infected = time - node.time_infected + days_since_infected = time - node.infection.time_infected outside_household_contacts = -1 local_contacts = -1 @@ -122,7 +125,7 @@ def increment(self, time): # will be thinned again local_infection_probs = self.get_infection_prob(local=True, infectious_age=days_since_infected, - asymptomatic=node.asymptomatic) + asymptomatic=node.infection.asymptomatic) local_infective_contacts = npr.binomial(local_contacts, local_infection_probs) @@ -145,12 +148,12 @@ def increment(self, time): self.new_within_household_infection(time=time, infecting_node=node) # Update how many contacts the node made - node.outside_house_contacts_made += outside_household_contacts + node.infection.outside_house_contacts_made += outside_household_contacts # How many outside household contacts cause new infections global_infection_probs = self.get_infection_prob(local=False, infectious_age=days_since_infected, - asymptomatic=node.asymptomatic) + asymptomatic=node.infection.asymptomatic) outside_household_new_infections = npr.binomial( outside_household_contacts, global_infection_probs @@ -160,7 +163,7 @@ def increment(self, time): self.new_outside_household_infection(time=time, infecting_node=node) node_time_tuple = (self.network.node_count, time) - node.spread_to_global_node_time_tuples.append(node_time_tuple) + node.infection.spread_to_global_node_time_tuples.append(node_time_tuple) def contacts_made_today(self, household_size) -> int: """Generates the number of contacts made today by a node, given the house size of the node. @@ -174,7 +177,8 @@ def contacts_made_today(self, household_size) -> int: obs = sum([int(cdf[i] < random) for i in range(100)]) return obs - def compute_hh_infection_probs(self, pairwise_survival_prob: float) -> np.ndarray: + @staticmethod + def compute_hh_infection_probs(pairwise_survival_prob: float) -> np.ndarray: # Precomputing the infection probabilities for the within household epidemics. contact_prob = 0.8 day_0_infection_prob = current_hazard_rate(0, pairwise_survival_prob) / contact_prob @@ -269,5 +273,5 @@ def perform_recoveries(self, time: int): recovered state """ for node in self.network.all_nodes(): - if node.recovery_time == time: - node.recovered = True + if node.infection.recovery_time == time: + node.infection.recovered = True diff --git a/household_contact_tracing/intervention.py b/household_contact_tracing/intervention.py index a09d907..5d56f27 100644 --- a/household_contact_tracing/intervention.py +++ b/household_contact_tracing/intervention.py @@ -54,9 +54,9 @@ def lft_nodes(self, time: int, prob_lfa_positive: Callable) -> List[Node]: positive_nodes = [] for node in self.network.all_nodes(): - if node.being_lateral_flow_tested: + if node.lfd_testing.being_lateral_flow_tested: if node.will_lfa_test_today(self.node_daily_prob_lfa_test): - if not node.received_positive_test_result: + if not node.tracing.received_positive_test_result: if node.lfa_test_node(time, prob_lfa_positive): positive_nodes.append(node) return positive_nodes @@ -77,20 +77,20 @@ def completed_quarantine(self, time: int): # (if they do not self-report they will not isolate; if contact traced, they will be # quarantining for the quarantine duration) # if node.household_id == node.infected_by_node().household_id: - if node.infecting_node: - if (node.infection_status(time) == InfectionStatus.unknown_infection) & node.isolated: + if node.infection.infecting_node_id: + if (node.infection_status(time) == InfectionStatus.unknown_infection) & node.infection.isolated: if node.locally_infected(): if time >= (node.household.earliest_recognised_symptom_onset(model_time=time) + self.quarantine_duration): node.isolated = False - node.completed_isolation = True + node.tracing.completed_isolation = True # For nodes who do not self-report, and are not in the same household as # their infector (if they do not self-report they will not isolate; if contact # traced, they will be quarantining for the quarantine duration) - elif node.contact_traced & (time >= node.time_infected + self.quarantine_duration): - node.isolated = False - node.completed_isolation = True + elif node.tracing.contact_traced & (time >= node.infection.time_infected + self.quarantine_duration): + node.infection.isolated = False + node.tracing.completed_isolation = True def completed_isolation(self, time: int): """ @@ -106,20 +106,20 @@ def completed_isolation(self, time: int): """ for node in self.network.all_nodes(): - if node.isolated: + if node.infection.isolated: infection_status = node.infection_status(time) if infection_status in [InfectionStatus.known_infection, InfectionStatus.self_recognised_infection]: - if node.avenue_of_testing == TestType.lfa: - if time >= node.positive_test_time + self.self_isolation_duration: - node.isolated = False - node.completed_isolation = True + if node.lfd_testing.avenue_of_testing == TestType.lfa: + if time >= node.lfd_testing.positive_test_time + self.self_isolation_duration: + node.infection.isolated = False + node.tracing.completed_isolation = True else: - if time >= node.symptom_onset_time + self.self_isolation_duration: + if time >= node.tracing.symptom_onset_time + self.self_isolation_duration: # this won't include nodes who tested positive due to LF tests who do not # have symptoms node.isolated = False - node.completed_isolation = True + node.tracing.completed_isolation = True def completed_lateral_flow_testing(self, time: int): """If a node is currently in lateral flow testing, and has completed this period then we @@ -137,7 +137,9 @@ def completed_lateral_flow_testing(self, time: int): """ for node in self.network.all_nodes(): - if time >= node.time_started_lfa_testing + self.lateral_flow_testing_duration \ - and node.being_lateral_flow_tested: - node.being_lateral_flow_tested = False - node.completed_lateral_flow_testing_time = time + # Todo: Check Ann's addition of node.lfd_testing.time_started_lfa_testing (not null) condition + if node.lfd_testing.time_started_lfa_testing and\ + time >= node.lfd_testing.time_started_lfa_testing + self.lateral_flow_testing_duration \ + and node.lfd_testing.being_lateral_flow_tested: + node.lfd_testing.being_lateral_flow_tested = False + node.lfd_testing.completed_lateral_flow_testing_time = time diff --git a/household_contact_tracing/network.py b/household_contact_tracing/network.py index cff569c..3a487b7 100644 --- a/household_contact_tracing/network.py +++ b/household_contact_tracing/network.py @@ -5,7 +5,8 @@ import numpy from dataclasses import dataclass -from household_contact_tracing.parameterised import Parameterised +from household_contact_tracing.node_attributes import LFDTestingAdherenceAttributes, TracingAdherenceAttributes, \ + ReturningTravellerAttributes, LFDTestingAttributes, TracingAttributes, InfectionAttributes class EdgeType(Enum): @@ -132,12 +133,12 @@ def all_households(self) -> Iterator[Household]: def active_infections(self): """Returns a list of nodes who have not yet recovered. - Nodes can still infect unless they have been isolated. + Nodes can still infect unless they have been isolated. Returns: list: list of nodes able to infect """ - return [node for node in self.all_nodes() if not node.recovered] + return [node for node in self.all_nodes() if not node.infection.recovered] def is_isomorphic(self, network: Network) -> bool: """ Determine whether graphs have identical network structures.""" @@ -177,7 +178,7 @@ def household(self, house_id: int) -> Household: def count_non_recovered_nodes(self) -> int: """Returns the number of nodes not in the recovered state.""" - return len([node for node in self.all_nodes() if not node.recovered]) + return len([node for node in self.all_nodes() if not node.infection.recovered]) def get_edge_between_household(self, house1: Household, house2: Household) -> Tuple[int, int]: """Get the id's of the two nodes that connect households.""" @@ -188,34 +189,29 @@ def get_edge_between_household(self, house1: Household, house2: Household) -> Tu def is_edge_app_traced(self, edge: Tuple[int, int]) -> bool: """Returns whether two nodes have the contract tracing app.""" - node_1_app = self.node(edge[0]).has_contact_tracing_app - node_2_app = self.node(edge[1]).has_contact_tracing_app + node_1_app = self.node(edge[0]).tracing.has_contact_tracing_app + node_2_app = self.node(edge[1]).tracing.has_contact_tracing_app return node_1_app and node_2_app - def add_node(self, time_infected, household_id, isolated, will_uptake_isolation, - propensity_imperfect_isolation, asymptomatic, symptom_onset_time, - pseudo_symptom_onset_time, recovery_time, will_report_infection, - time_of_reporting, has_contact_tracing_app, contact_traced, testing_delay=0, - additional_attributes: Optional[dict] = None, - infecting_node: Optional[Node] = None, completed_isolation=False) -> Node: + def add_node(self, household_id, + infection_attributes: Optional[dict] = None, + lfd_testing_adherence_attributes: Optional[dict] = None, + tracing_attributes: Optional[dict] = None, + tracing_adherence_attributes: Optional[dict] = None, + returning_travellers_attributes: Optional[dict] = None, + lfd_testing_attributes: Optional[dict] = None + ) -> Node: new_node_id = self.node_count + 1 self.graph.add_node(new_node_id) new_node_household = self.household(household_id) - node = Node(node_id=new_node_id, time_infected=time_infected, - household=new_node_household, isolated=isolated, - will_uptake_isolation=will_uptake_isolation, - propensity_imperfect_isolation=propensity_imperfect_isolation, - asymptomatic=asymptomatic, symptom_onset_time=symptom_onset_time, - pseudo_symptom_onset_time=pseudo_symptom_onset_time, - recovery_time=recovery_time, - will_report_infection=will_report_infection, - time_of_reporting=time_of_reporting, - has_contact_tracing_app=has_contact_tracing_app, - contact_traced=contact_traced, - testing_delay=testing_delay, - additional_attributes=additional_attributes, - infecting_node=infecting_node, - completed_isolation=completed_isolation) + node = Node(node_id=new_node_id, + household=new_node_household, + infection_attributes=infection_attributes, + lfd_testing_adherence_attributes=lfd_testing_adherence_attributes, + tracing_attributes=tracing_attributes, + tracing_adherence_attributes=tracing_adherence_attributes, + returning_travellers_attributes=returning_travellers_attributes, + lfd_testing_attributes=lfd_testing_attributes) self.graph.nodes[new_node_id]['node_obj'] = node return node @@ -238,7 +234,7 @@ def label_edges_inside_household(self, household: Household, new_edge_type: Edge self.graph.edges[edge[0], edge[1]].update({"edge_type": new_edge_type}) -class Node(Parameterised): +class Node: """ A class used to store contact tracing node data. Uses networkx as storage tool. @@ -256,75 +252,48 @@ class Node(Parameterised): """ - def __init__(self, node_id: int, time_infected: int, household: Household, isolated: bool, - will_uptake_isolation: bool, propensity_imperfect_isolation: bool, - asymptomatic: bool, symptom_onset_time: float, pseudo_symptom_onset_time: int, - recovery_time: int, will_report_infection: bool, time_of_reporting: int, - has_contact_tracing_app: bool, contact_traced: bool, testing_delay: int = 0, - completed_isolation=False, outside_house_contacts_made=0, recovered=False, - infecting_node: Optional[Node] = None, additional_attributes: dict = None): + def __init__(self, node_id: int, household: Household, + infection_attributes: dict = None, + lfd_testing_adherence_attributes: dict = None, + tracing_attributes: dict = None, + tracing_adherence_attributes: dict = None, + returning_travellers_attributes: dict = None, + lfd_testing_attributes: dict = None): self.id = node_id - self.time_infected = time_infected self.household = household - self.isolated = isolated - self.will_uptake_isolation = will_uptake_isolation - self.propensity_imperfect_isolation = propensity_imperfect_isolation - self.asymptomatic = asymptomatic - self.symptom_onset_time = symptom_onset_time - self.pseudo_symptom_onset_time = pseudo_symptom_onset_time - self.recovery_time = recovery_time - self.will_report_infection = will_report_infection - self.time_of_reporting = time_of_reporting - self.has_contact_tracing_app = has_contact_tracing_app - self.testing_delay = testing_delay - self.contact_traced = contact_traced - self.outside_house_contacts_made = outside_house_contacts_made - self.spread_to_global_node_time_tuples = [] - self.recovered = recovered - self.propagated_contact_tracing = False - self.infecting_node = infecting_node if infecting_node else None - self.completed_isolation = completed_isolation - self.received_result = False - self.received_positive_test_result = False - - self.being_lateral_flow_tested = None - self.time_started_lfa_testing = None - self.avenue_of_testing: Optional[TestType] = None - self.positive_test_time = None - self.node_will_take_up_lfa_testing = None - self.confirmatory_PCR_result_was_positive: Optional[bool] = None - self.taken_confirmatory_PCR_test: Optional[bool] = None - self.confirmatory_PCR_test_result_time = None - self.propensity_risky_behaviour_lfa_testing = None - self.propensity_to_miss_lfa_tests = None - - # Update instance variables with anything in `additional_attributes` - self.update_params(additional_attributes) + + # Instantiate attribute classes + self.infection = InfectionAttributes(infection_attributes) + self.lfd_testing_adherence = LFDTestingAdherenceAttributes(lfd_testing_adherence_attributes) + self.tracing = TracingAttributes(tracing_attributes) + self.tracing_adherence = TracingAdherenceAttributes(tracing_adherence_attributes) + self.returning_travellers = ReturningTravellerAttributes(returning_travellers_attributes) + self.lfd_testing = LFDTestingAttributes(lfd_testing_attributes) def time_relative_to_symptom_onset(self, time: int) -> int: # asymptomatics do not have a symptom onset time # pseudo_symptom_onset time is a fake onset we give them # so we can work out when they test positive - return time - self.pseudo_symptom_onset_time + return time - self.returning_travellers.pseudo_symptom_onset_time def locally_infected(self) -> bool: - if self.infecting_node: - return self.infecting_node.household == self.household + if self.infection.infecting_node_id: + return self.household.network.node(self.infection.infecting_node_id).household == self.household else: return False def infection_status(self, time_now: int) -> InfectionStatus: - if self.contact_traced: - if self.symptom_onset_time + self.testing_delay <= time_now: + if self.tracing.contact_traced: + if self.tracing.symptom_onset_time + self.tracing.testing_delay <= time_now: return InfectionStatus.known_infection - if self.symptom_onset_time <= time_now: + if self.tracing.symptom_onset_time <= time_now: return InfectionStatus.self_recognised_infection else: - if self.will_report_infection: - if self.time_of_reporting + self.testing_delay <= time_now: + if self.tracing.will_report_infection: + if self.tracing.time_of_reporting + self.tracing.testing_delay <= time_now: return InfectionStatus.known_infection - if self.time_of_reporting <= time_now: + if self.tracing.time_of_reporting <= time_now: return InfectionStatus.self_recognised_infection return InfectionStatus.unknown_infection @@ -334,32 +303,32 @@ def node_type(self, time=None) -> NodeType: params time (int): The current increment / step number (e.g. day number) of the simulation """ - if self.being_lateral_flow_tested: - if self.isolated: + if self.lfd_testing.being_lateral_flow_tested: + if self.infection.isolated: return NodeType.being_lateral_flow_tested_isolated else: return NodeType.being_lateral_flow_tested_not_isolated - elif self.isolated: + elif self.infection.isolated: return NodeType.isolated - elif not self.asymptomatic: - if self.will_report_infection: + elif not self.infection.asymptomatic: + if self.tracing.will_report_infection: return NodeType.symptomatic_will_report_infection else: return NodeType.symptomatic_will_not_report_infection - elif self.received_positive_test_result: - if self.avenue_of_testing == TestType.pcr: + elif self.tracing.received_positive_test_result: + if self.lfd_testing.avenue_of_testing == TestType.pcr: return NodeType.received_pos_test_pcr else: return NodeType.received_pos_test_lfa - elif self.received_result and self.avenue_of_testing == TestType.pcr: + elif self.tracing.received_result and self.lfd_testing.avenue_of_testing == TestType.pcr: return NodeType.received_neg_test_pcr - elif self.taken_confirmatory_PCR_test: - if time and time >= self.confirmatory_PCR_test_result_time: - if self.confirmatory_PCR_result_was_positive: + elif self.lfd_testing.taken_confirmatory_PCR_test: + if time and time >= self.lfd_testing.confirmatory_PCR_test_result_time: + if self.lfd_testing_adherence.confirmatory_PCR_result_was_positive: return NodeType.confirmatory_pos_pcr_test else: return NodeType.confirmatory_neg_pcr_test - elif self.asymptomatic: + elif self.infection.asymptomatic: return NodeType.asymptomatic else: return NodeType.default @@ -367,20 +336,20 @@ def node_type(self, time=None) -> NodeType: def take_confirmatory_pcr_test(self, time: int, prob_pcr_positive: Callable): """Given a the time relative to a nodes symptom onset, will that node test positive.""" - infectious_age_when_tested = time - self.time_infected + infectious_age_when_tested = time - self.infection.time_infected - self.confirmatory_PCR_test_result_time = time + self.testing_delay - self.taken_confirmatory_PCR_test = True + self.lfd_testing.confirmatory_PCR_test_result_time = time + self.tracing.testing_delay + self.lfd_testing.taken_confirmatory_PCR_test = True if numpy.random.binomial(1, prob_pcr_positive(infectious_age_when_tested)) == 1: - self.confirmatory_PCR_result_was_positive = True + self.lfd_testing_adherence.confirmatory_PCR_result_was_positive = True else: - self.confirmatory_PCR_result_was_positive = False + self.lfd_testing_adherence.confirmatory_PCR_result_was_positive = False def will_lfa_test_today(self, daily_prob_lfa_test: float) -> bool: """Determine whether a node will do an LFT test today.""" - if not self.propensity_to_miss_lfa_tests: + if not self.lfd_testing_adherence.propensity_to_miss_lfa_tests: return True if numpy.random.binomial(1, daily_prob_lfa_test) == 1: @@ -390,7 +359,7 @@ def will_lfa_test_today(self, daily_prob_lfa_test: float) -> bool: def lfa_test_node(self, time: int, prob_lfa_positive: Callable): """Given a the time relative to a nodes symptom onset, will that node test positive""" - infectious_age = time - self.time_infected + infectious_age = time - self.infection.time_infected prob_positive_result = prob_lfa_positive(infectious_age) @@ -468,7 +437,7 @@ def get_recognised_symptom_onsets(self, model_time: int): infection_status = household_node.infection_status(model_time) if infection_status in [InfectionStatus.known_infection, InfectionStatus.self_recognised_infection]: - recognised_symptom_onsets.append(household_node.symptom_onset_time) + recognised_symptom_onsets.append(household_node.tracing.symptom_onset_time) return recognised_symptom_onsets def get_positive_test_times(self, model_time: int) -> List[int]: @@ -476,8 +445,8 @@ def get_positive_test_times(self, model_time: int) -> List[int]: for node in self.nodes: if node.infection_status(model_time) == InfectionStatus.known_infection: - if node.received_positive_test_result: - positive_test_times.append(node.positive_test_time) + if node.tracing.received_positive_test_result: + positive_test_times.append(node.lfd_testing.positive_test_time) return positive_test_times def earliest_recognised_symptom_onset(self, model_time: int): @@ -515,9 +484,9 @@ def isolate_household(self, time: int): # Update isolated and contact traced status for Nodes in Household for node in self.nodes: - node.contact_traced = True - if node.will_uptake_isolation: - node.isolated = True + node.tracing.contact_traced = True + if node.tracing_adherence.will_uptake_isolation: + node.infection.isolated = True self._update_edges_on_isolation() @@ -548,11 +517,11 @@ def start_lateral_flow_testing_household(self, time: int): self.being_lateral_flow_tested_start_time = time for node in self.nodes: - if node.node_will_take_up_lfa_testing: - if not node.received_positive_test_result: - if not node.being_lateral_flow_tested: - node.being_lateral_flow_tested = True - node.time_started_lfa_testing = time + if node.lfd_testing_adherence.node_will_take_up_lfa_testing: + if not node.tracing.received_positive_test_result: + if not node.lfd_testing.being_lateral_flow_tested: + node.lfd_testing.being_lateral_flow_tested = True + node.lfd_testing.time_started_lfa_testing = time def start_lateral_flow_testing_household_and_quarantine(self, time): """Sets the household to the lateral flow testing status so that new within household @@ -566,14 +535,14 @@ def start_lateral_flow_testing_household_and_quarantine(self, time): self.contact_traced = True for node in self.nodes: - if node.node_will_take_up_lfa_testing: - if not node.received_positive_test_result: - if not node.being_lateral_flow_tested: - node.being_lateral_flow_tested = True - node.time_started_lfa_testing = time + if node.lfd_testing_adherence.node_will_take_up_lfa_testing: + if not node.tracing.received_positive_test_result: + if not node.lfd_testing.being_lateral_flow_tested: + node.lfd_testing.being_lateral_flow_tested = True + node.lfd_testing.time_started_lfa_testing = time - if node.will_uptake_isolation: - node.isolated = True + if node.tracing_adherence.will_uptake_isolation: + node.infection.isolated = True def apply_positive_policy(self, time: int, household_positive_policy: str): """Depending on the positive policy, different interventions are made to the household @@ -614,7 +583,7 @@ def update_network(self): # Update the nodes to the contact traced status for node in self.nodes: - node.contact_traced = True + node.tracing.contact_traced = True # Colour the edges within household self.network.label_edges_inside_household(self, EdgeType.within_house) @@ -622,7 +591,7 @@ def update_network(self): def isolate_if_symptomatic_nodes(self, time: int): """If there are any symptomatic nodes in the household then isolate the household.""" for node in self.nodes: - if node.symptom_onset_time <= time and not node.completed_isolation: + if node.tracing.symptom_onset_time <= time and not node.tracing.completed_isolation: self.isolate_household(time) break @@ -630,5 +599,20 @@ def quarantine_traced_node(self): traced_node = self.find_traced_node() # the traced node should go into quarantine - if not traced_node.isolated and traced_node.will_uptake_isolation: - traced_node.isolated = True + if not traced_node.infection.isolated and traced_node.tracing_adherence.will_uptake_isolation: + traced_node.infection.isolated = True + + @property + def local_epidemic_completed(self): + """ + Returns true if all infections in the household have recovered, + which is defined as being 10 + """ + return all([node.infection.recovered for node in self.nodes]) + + @property + def household_epidemic_size(self): + """Returns the current size of the household epidemic, i.e: the number of household members + that are, or were, infected. + """ + return self.size - self.susceptibles diff --git a/household_contact_tracing/node_attributes.py b/household_contact_tracing/node_attributes.py new file mode 100644 index 0000000..5d85f24 --- /dev/null +++ b/household_contact_tracing/node_attributes.py @@ -0,0 +1,158 @@ +from typing import Optional, List, Tuple +from household_contact_tracing.parameterised import Parameterised + + +class InfectionAttributes(Parameterised): + """ + A class used to store Node attributes relating to infection + Inherits from Parameterised to handle validation and updating of large number of parameters + + Attributes + ---------- + asymptomatic + infecting_node_id + isolated + outside_house_contacts_made + recovered + recovery_time + spread_to_global_node_time_tuples + time_infected + """ + + def __init__(self, attributes: dict): + self.asymptomatic: Optional[bool] = None + self.infecting_node_id: Optional[int] = None + self.isolated: Optional[bool] = None + self.outside_house_contacts_made: int = 0 + self.recovered: bool = False + self.recovery_time: Optional[int] = None + self.spread_to_global_node_time_tuples: List[Tuple[int, int]] = [] + self.time_infected: Optional[int] = None + + # Update instance variables with anything in attributes + self.update_params(attributes) + + +class LFDTestingAttributes(Parameterised): + """ + A class used to store Node attributes relating to LFD Testing + Inherits from Parameterised to handle validation and updating of large number of parameters + + Attributes + ---------- + avenue_of_testing + being_lateral_flow_tested + positive_test_time + taken_confirmatory_PCR_test + time_started_lfa_testing + confirmatory_PCR_test_result_time + completed_lateral_flow_testing_time + lateral_flow_testing_duration + """ + + def __init__(self, attributes: dict): + self.avenue_of_testing: Optional[int] = None + self.being_lateral_flow_tested: Optional[bool] = None + self.positive_test_time: Optional[int] = None + self.taken_confirmatory_PCR_test: Optional[bool] = None + self.time_started_lfa_testing: Optional[int] = None + self.confirmatory_PCR_test_result_time: Optional[int] = None + self.completed_lateral_flow_testing_time: Optional[bool] = None + self.lateral_flow_testing_duration: Optional[int] = 0 + + # Update instance variables with anything in attributes + self.update_params(attributes) + + +class LFDTestingAdherenceAttributes(Parameterised): + """ + A class used to store Node attributes relating to LFD Testing Adherence + Inherits from Parameterised to handle validation and updating of large number of parameters + + Attributes + ---------- + confirmatory_PCR_result_was_positive + node_will_take_up_lfa_testing + propensity_risky_behaviour_lfa_testing + propensity_to_miss_lfa_tests + """ + + def __init__(self, attributes: dict): + self.confirmatory_PCR_result_was_positive: Optional[bool] = None + self.node_will_take_up_lfa_testing: Optional[bool] = None + self.propensity_risky_behaviour_lfa_testing: Optional[float] = None + self.propensity_to_miss_lfa_tests: Optional[float] = None + + # Update instance variables with anything in attributes + self.update_params(attributes) + + +class ReturningTravellerAttributes(Parameterised): + """ + A class used to store Node attributes relating to Returning Travellers + Inherits from Parameterised to handle validation and updating of large number of parameters + + Attributes + ---------- + pseudo_symptom_onset_time + """ + + def __init__(self, attributes: dict): + self.pseudo_symptom_onset_time: Optional[int] = None + + # Update instance variables with anything in attributes + self.update_params(attributes) + + +class TracingAttributes(Parameterised): + """ + A class used to store Node attributes relating to Contact Tracing + Inherits from Parameterised to handle validation and updating of large number of parameters + + Attributes + ---------- + contact_traced + has_contact_tracing_app + propagated_contact_tracing + received_positive_test_result + received_result + symptom_onset_time + testing_delay + time_of_reporting + will_report_infection + completed_isolation + """ + + def __init__(self, attributes: dict): + self.contact_traced: Optional[bool] = None + self.has_contact_tracing_app: Optional[bool] = None + self.propagated_contact_tracing: bool = False + self.received_positive_test_result: bool = False + self.received_result: bool = False + self.symptom_onset_time: Optional[int] = None + self.testing_delay: Optional[int] = None + self.time_of_reporting: Optional[int] = None + self.will_report_infection: Optional[bool] = None + self.completed_isolation: Optional[bool] = None + + # Update instance variables with anything in attributes + self.update_params(attributes) + + +class TracingAdherenceAttributes(Parameterised): + """ + A class used to store Node attributes relating to Contact Tracing Adherence + Inherits from Parameterised to handle validation and updating of large number of parameters + + Attributes + ---------- + propensity_imperfect_isolation + will_uptake_isolation + """ + + def __init__(self, attributes: dict): + self.propensity_imperfect_isolation: Optional[float] = None + self.will_uptake_isolation: Optional[bool] = None + + # Update instance variables with anything in attributes + self.update_params(attributes) diff --git a/household_contact_tracing/queueing_processes.py b/household_contact_tracing/queueing_processes.py new file mode 100644 index 0000000..e6b19fb --- /dev/null +++ b/household_contact_tracing/queueing_processes.py @@ -0,0 +1,654 @@ +''' +Contains queueing process objects that are used to model testing delays and probability of +not being able to get testing when there are constrained processing resource. +Processing resources can refer to either swabbing capacity, or genetic sequencing capacity. +''' + +from datetime import time +import pandas as pd +import numpy as np +import numpy.random as npr +from typing import Callable, List + +class Queue: + + def __init__( + self, + days_to_simulate: int, + capacity: List[int] + ): + """Creates a queueing process object that stores the current and previous states of the queue. + Also contains methods for accessing the data in the queue, and changing the data in the queue. + + Args: + days_to_simulate (int): The total number of days that will be simulated + capacity (list): The processing capacity (integer values) at each timepoint + max_time_in_queue (int): Maximum days from symptom onset to ineligibility for processing + verbose (bool, optional): If true prints some outputs. Defaults to False. + """ + + self.days_to_simulate = days_to_simulate + self.capacity = capacity + + # default values + self.time = 0 + + # initialise a dataframe that stores summaries of the queue at each timepoint + self.create_queue_df() + self.create_applicants_df() + + + def create_queue_df(self): + """ + Queue df stores summaries of the overall status of the queue at each timepoint. It is updated as the calculations progresses + """ + + # create a dataframe to store information about the overall queueing process + self.queue_df = pd.DataFrame({ + 'time': list(range(self.days_to_simulate)), + 'capacity': self.capacity + }) + + # create some empty columns for storing results + self.queue_df['new_applicants'] = '' + self.queue_df['spillover_to_next_day'] = '' + self.queue_df['total_applications_today'] = '' + self.queue_df['capacity_exceeded'] = '' + self.queue_df['capacity_exceeded_by'] = '' + self.queue_df['number_processed_today'] = '' + self.queue_df['number_left_queue_not_tested'] = '' + + + def create_applicants_df(self): + """Applicants df store information about everyone who has applied for a test at each timepoint. It is updated as the calculation progresses. + """ + + self.applicant_df = pd.DataFrame() + + # create empty columns for applicants + self.applicant_df['id'] = '' + self.applicant_df['processed'] = '' + self.applicant_df['waiting_to_be_processed'] = '' + self.applicant_df['left_queue_not_processed'] = '' + self.applicant_df['time_symptom_onset'] = '' + self.applicant_df['time_joined_queue'] = '' + self.applicant_df['time_processed'] = '' + self.applicant_df['time_received_result'] = '' + self.applicant_df['time_will_leave_queue'] = '' + + def add_new_applicants( + self, + ids: list, + time: int, + symptom_onset_times: list, + queue_leaving_times: list + ): + """ + Adds new applicants to the queue. + """ + + new_applicant_df = pd.DataFrame( + { + 'id': ids, + 'time_symptom_onset': symptom_onset_times, + 'time_will_leave_queue': queue_leaving_times + } + ) + + # initialise other columns with default values + new_applicant_df['processed'] = False + new_applicant_df['waiting_to_be_processed'] = True # default value, initially the queue is empty + new_applicant_df['left_queue_not_processed'] = '' + new_applicant_df['time_joined_queue'] = time + new_applicant_df['time_processed'] = '' + new_applicant_df['time_received_result'] = '' + + self.applicant_df = self.applicant_df.append(new_applicant_df, ignore_index = True) + + def swab_applicants(self, + to_be_processed: list, + processing_delays: list): + """For a list of applicants who were successful in getting thorugh the queue, update their variables associated with processing + Args: + to_be_processed (list): A list of integers, referring the rows of the applicant_dataframe that will get processed + """ + + # The columns that will be updated + columns_to_update = [ + 'waiting_to_be_processed', + 'time_processed', + 'left_queue_not_processed', + 'processed' + ] + + # record an attribtue of which individuals were processed today for use later + self.todays_processed_index = to_be_processed + + # update the above status to show they have been processed + self.applicant_df.loc[to_be_processed, columns_to_update] = [False, self.time, False, True] + + # work out when they receive their result, and update the data + self.applicant_df.loc[to_be_processed, 'time_received_result'] = self.time + np.array(processing_delays) + + # update the queue_df table with the number of individuals processed today + self.queue_df.loc[self.queue_df.time == self.time, ['number_processed_today']] = len(to_be_processed) + + + def update_queue_leaver_status(self): + """These individuals have been in the queue too long. They are no longer trying/able to get a swab. + """ + + # These people will leave the queue today + self.leavers = (self.applicant_df.time_will_leave_queue <= self.time) & (self.applicant_df.waiting_to_be_processed == True) + + + # record the number of people who carry over to the next day + if self.todays_capacity > len(self.current_applicants): + spillover_to_next_day = 0 + else: + spillover_to_next_day = len(self.current_applicants) - sum(self.leavers) - self.todays_capacity + + # + self.queue_df.loc[self.time, ['spillover_to_next_day', 'number_left_queue_not_tested']] = [spillover_to_next_day, sum(self.leavers)] + + # Set their waiting to be processed status to False + self.applicant_df.loc[self.leavers, ['waiting_to_be_processed', 'left_queue_not_processed']] = [False, True] + + + @property + def current_applicants(self) -> list: + """Gets the indexes of individuals waiting to be processed. + + Returns: + list: The indexes of individuals waiting to be processed + """ + return list(self.applicant_df[self.applicant_df.waiting_to_be_processed].index) + + + @property + def todays_capacity(self) -> int: + """Gets the number of processes that can be performed today. + + Returns: + int: The number of processes that can be performed today + """ + return int(self.queue_df[self.queue_df.time == self.time].capacity) + + + @property + def number_processes_performed_today(self) -> int: + """ + Gets the number of processes that have been completed. This will be the number of new applicants + or the processing capacity. + + Returns: + int: The number of processes that have been performed today + """ + + return sum(self.applicant_df.time_processed == self.time) + + +# controller layout +class QueueController: + + def __init__(self): + """ + Queue controllers control how the queue is run. In some + cases you may want to run the queue on it's own, or you may want + the queue to be interacting with a branching process model in various ways. + """ + + self.completed = False + + def process_todays_new_demand(self): + pass + + def process_queue(self): + pass + + def simulate_one_day(self): + pass + + +class DeterministicQueue(QueueController): + # TODO: Rename, it's not deterministic, but the inputs are + + def __init__( + self, + days_to_simulate: int, + demand: List[int], + capacity: List[int], + max_time_in_queue: int, + processing_delay_dist: Callable, + symptom_onset_delay_dist: Callable, + selection_method: str): + """A simple queueing process object that does not interact with a branching process model. + + The test demand and capacity are pre-determined, and the model works out what happens to the queue. + + Args: + days_to_simulate (int): Number of simulation steps to be performed + demand (List[int]): The number of new test seekers at each time step. + capacity (List[int]): The processing capacity of the queue at each time step. + max_time_in_queue (int): How long since symptom onset that an individual can remain in the queue + they become ineligible for testing + processing_delay_dist (Callable): A callable that returns integer test processing delays + symptom_onset_delay_dist (Callable): A callable the returns integer delays of the time from symptom onset to booking a test. + selection_method ('uniform', 'newest'): Method for selecting which applicants to process when demand exceeds capacity. + """ + + # initialise the queue + self.queue = Queue( + days_to_simulate = days_to_simulate, + capacity = capacity + ) + + # set parameters + self.demand = demand + self.processing_delay_dist = processing_delay_dist + self.symptom_onset_delay_dist = symptom_onset_delay_dist + self.max_time_in_queue = max_time_in_queue + self.days_to_simulate = days_to_simulate + self.selection_method = selection_method + + # ease of acccess stuff + self.time = self.queue.time + + + def add_new_test_seekers(self): + """ + Adds new test seekers to the queue. + + For this model, the new test seekers at each time point are defined a priori. + """ + symptom_onset_times = [ + self.time - self.symptom_onset_delay_dist() for _ in range(self.demand[self.time]) + ] + + queue_leaving_times = [ + onset_time + self.max_time_in_queue for onset_time in symptom_onset_times + ] + + self.queue.add_new_applicants( + ids = [''] * self.demand[self.time], + time = self.time, + symptom_onset_times = symptom_onset_times, + queue_leaving_times = queue_leaving_times + ) + + def select_applicants_for_processing(self, remaining_processing_capacity: int) -> list: + """Given the current demand and remaining testing capacity, compute which individuals get selected for testing. + + Args: + current_queue_applicants (list): A list of id's of individuals who are waiting to get processed. + remaining_processing_capacity (int): The remaining capacity for individuals to get processed. + + Returns: + list: The list of processed individuals. + """ + + if self.selection_method == 'uniform': + return( + npr.choice( + a = self.queue.current_applicants, + size = remaining_processing_capacity, + replace = False + ) + ) + elif self.selection_method == 'newest': + return( + self.queue.applicant_df.sort_values('time_joined_queue')[0:remaining_processing_capacity] + ) + + def process_queue(self): + """ + Performs processing of individuals up to capacity, and updates the dataframes that store the calculations. + """ + + # Note: this method is set up so that it can be called multiple times in one day + # in case new applicants are added multiple times in a day. This is sometimes useful + + number_applicants = len(self.queue.current_applicants) + + # update queue_df with the number of applicants today + self.queue.queue_df.loc[self.queue.queue_df.time == self.time, ['total_applications_today']] = [number_applicants] + + # how much processing capacity do we have remaining? The method + remaining_processing_capacity = self.queue.todays_capacity - self.queue.number_processes_performed_today + + # is todays remaining capacity exceeded? + if number_applicants <= remaining_processing_capacity: + # if capacity not exceeded, then everyone gets processed + + processing_delays = [ + self.processing_delay_dist() for _ in range(number_applicants) + ] + + self.queue.swab_applicants( + to_be_processed = self.queue.current_applicants, + processing_delays = processing_delays) + + else: + # Then processing capacity is being exceeded. We process up to capacity. + # We must select who gets processed, at the moment there is only one method + # implemented that does this, that picks a subset without replacement + + self.select_applicants_for_processing(remaining_processing_capacity) + + successful_applicants = npr.choice( + a = self.queue.current_applicants, + size = remaining_processing_capacity, + replace = False + ) + + processing_delays = [ + self.processing_delay_dist() for _ in range(remaining_processing_capacity) + ] + + self.queue.swab_applicants( + to_be_processed = successful_applicants, + processing_delays = processing_delays) + + def update_queue_leaver_status(self): + """These individuals have been in the queue too long. They are no longer trying/able to get a swab. + """ + + # These people will leave the queue today + self.leavers = (self.queue.applicant_df.time_will_leave_queue <= self.time) & (self.queue.applicant_df.waiting_to_be_processed == True) + + self.queue.queue_df.loc[self.time, 'number_left_queue_not_tested'] = [sum(self.leavers)] + + # Set their waiting to be processed status to False + self.queue.applicant_df.loc[self.leavers, ['waiting_to_be_processed', 'left_queue_not_processed']] = [False, True] + + # work out who will come back the next day + # not left and not processed + returners_index = self.queue.applicant_df.waiting_to_be_processed == True + + self.queue.queue_df.loc[self.time, 'spillover_to_next_day'] = [sum(returners_index)] + + def simulate_one_day(self): + """ + Simulates one day of the queue. + """ + + # steps required to simulate one day + self.add_new_test_seekers() + self.update_queue_leaver_status() + self.process_queue() + + self.queue.time += 1 + + def run_simulation(self): + """Runs the queueing process model. + """ + + while self.time < self.days_to_simulate: + + self.simulate_one_day() + + self.time += 1 + + +class DeterministicQueueVariantSequencing(QueueController): + + def __init__( + self, + days_to_simulate: int, + demand: List[int], + demand_variant: List[int], + capacity: List[int], + max_time_in_queue: int, + processing_delay_dist: Callable, + symptom_onset_delay_dist: Callable, + selection_method: str): + """A simple queueing process object that does not interact with a branching process model. + + The test demand and capacity are pre-determined, and the model works out what happens to the queue. + + Args: + days_to_simulate (int): Number of simulation steps to be performed + demand (List[int]): The number of new test seekers at each time step. + capacity (List[int]): The processing capacity of the queue at each time step. + max_time_in_queue (int): How long since symptom onset that an individual can remain in the queue + they become ineligible for testing + processing_delay_dist (Callable): A callable that returns integer test processing delays + symptom_onset_delay_dist (Callable): A callable the returns integer delays of the time from symptom onset to booking a test. + selection_method ('uniform', 'newest'): Method for selecting which applicants to process when demand exceeds capacity. + """ + + # initialise the queue + self.queue = Queue( + days_to_simulate = days_to_simulate, + capacity = capacity + ) + + # add an empty column to store variant status + self.queue.applicant_df['variant'] = '' + + # set parameters + self.demand = demand + self.demand_variant = demand_variant + self.processing_delay_dist = processing_delay_dist + self.symptom_onset_delay_dist = symptom_onset_delay_dist + self.max_time_in_queue = max_time_in_queue + self.days_to_simulate = days_to_simulate + self.selection_method = selection_method + + # ease of acccess stuff + self.time = self.queue.time + + + def add_new_queue_joiners(self): + """ + Adds new test seekers to the queue. + + For this model, the new test seekers at each time point are defined a priori. + """ + total_new_joiners = self.demand[self.time] + self.demand_variant[self.time] + + symptom_onset_times = [ + self.time - self.symptom_onset_delay_dist() for _ in range(total_new_joiners) + ] + + queue_leaving_times = [ + onset_time + self.max_time_in_queue for onset_time in symptom_onset_times + ] + + self.queue.add_new_applicants( + ids = [''] * total_new_joiners, + time = self.time, + symptom_onset_times = symptom_onset_times, + queue_leaving_times = queue_leaving_times + ) + + # work out which of the new joiners are variants + variant_ids = npr.choice( + a = list(range(total_new_joiners)), + size = self.demand_variant[self.time], + replace = False) + + # by default cases are not variants + variant_status = [False]*total_new_joiners + for _ in variant_ids: + variant_status[_] = True + + # set the variant status on the applicant dataframe column + todays_joiner_index = self.queue.applicant_df.time_joined_queue == self.time + self.queue.applicant_df.loc[todays_joiner_index, 'variant'] = variant_status + + + def select_applicants_for_processing(self, remaining_processing_capacity: int) -> list: + """Given the current demand and remaining testing capacity, compute which individuals get selected for testing. + + Args: + current_queue_applicants (list): A list of id's of individuals who are waiting to get processed. + remaining_processing_capacity (int): The remaining capacity for individuals to get processed. + + Returns: + list: The list of processed individuals. + """ + + if self.selection_method == 'uniform': + return( + npr.choice( + a = self.queue.current_applicants, + size = remaining_processing_capacity, + replace = False + ) + ) + elif self.selection_method == 'newest': + return( + self.queue.applicant_df.sort_values('time_joined_queue')[0:remaining_processing_capacity] + ) + + def process_queue(self): + """ + Performs processing of individuals up to capacity, and updates the dataframes that store the calculations. + """ + + # Note: this method is set up so that it can be called multiple times in one day + # in case new applicants are added multiple times in a day. This is sometimes useful + + number_applicants = len(self.queue.current_applicants) + + # update queue_df with the number of applicants today + self.queue.queue_df.loc[self.queue.queue_df.time == self.time, ['total_applications_today']] = [number_applicants] + + # how much processing capacity do we have remaining? The method + remaining_processing_capacity = self.queue.todays_capacity - self.queue.number_processes_performed_today + + # is todays remaining capacity exceeded? + if number_applicants <= remaining_processing_capacity: + # if capacity not exceeded, then everyone gets processed + + processing_delays = [ + self.processing_delay_dist() for _ in range(number_applicants) + ] + + self.queue.swab_applicants( + to_be_processed = self.queue.current_applicants, + processing_delays = processing_delays) + + else: + # Then processing capacity is being exceeded. We process up to capacity. + # We must select who gets processed, at the moment there is only one method + # implemented that does this, that picks a subset without replacement + + self.select_applicants_for_processing(remaining_processing_capacity) + + successful_applicants = npr.choice( + a = self.queue.current_applicants, + size = remaining_processing_capacity, + replace = False + ) + + processing_delays = [ + self.processing_delay_dist() for _ in range(remaining_processing_capacity) + ] + + self.queue.swab_applicants( + to_be_processed = successful_applicants, + processing_delays = processing_delays) + + def update_queue_leaver_status(self): + """These individuals have been in the queue too long. They are no longer trying/able to get a swab. + """ + + # These people will leave the queue today + self.leavers = (self.queue.applicant_df.time_will_leave_queue <= self.time) & (self.queue.applicant_df.waiting_to_be_processed == True) + + self.queue.queue_df.loc[self.time, 'number_left_queue_not_tested'] = [sum(self.leavers)] + + # Set their waiting to be processed status to False + self.queue.applicant_df.loc[self.leavers, ['waiting_to_be_processed', 'left_queue_not_processed']] = [False, True] + + # work out who will come back the next day + # not left and not processed + returners_index = self.queue.applicant_df.waiting_to_be_processed == True + + self.queue.queue_df.loc[self.time, 'spillover_to_next_day'] = [sum(returners_index)] + + def simulate_one_day(self): + """ + Simulates one day of the queue. + """ + + # steps required to simulate one day + self.add_new_queue_joiners() + self.update_queue_leaver_status() + self.process_queue() + + self.queue.time += 1 + + def run_simulation(self): + """Runs the queueing process model. + """ + + while self.time < self.days_to_simulate: + + self.simulate_one_day() + + self.time += 1 + + +class QueueBranchingProcessController(): + + def __init__( + self, + queue: Queue): + + self.queue = queue + + def get_todays_queue_output(self): + """ + Provides outputs from the queueing process that can be passed to a + branching process model. + + Returns: + dict: output dict, with the ids and number of processed individuals + """ + + processed_individuals = self.queue.applicant_df.loc[self.queue.todays_processed_index] + + output = { + 'leaving_the_queue_node_ids': self.queue.todays_leavers, + 'processed_individuals': processed_individuals + } + + return output + +class QueueAnalyzer(): + + def __init__( + self, + queue: Queue): + + self.queue = queue + self.applicant_df = queue.applicant_df + self.queue_df = queue.queue_df + + def get_prob_getting_processed(self, time_joined_queue: int): + """ + Returns the probability of getting processed if you join the queue on a specified day + + Args: + time_joined_queue (int): The day of interest + """ + valid_individuals = (self.queue.applicant_df.time_joined_queue == time_joined_queue) & (self.applicant_df.waiting_to_be_processed == False) + left_queue_not_processed = self.applicant_df[valid_individuals].left_queue_not_processed + return 1 - left_queue_not_processed.mean() + + def get_delays_for(self, time_joined_queue: int, delay_from_column: str, delay_to_column: str): + """ + Return a list of the delays between two timepoints who joined on a specified day + + Args: + time (int): The day on which the applicants joined the queue + delay_from_column (str): The earliest timepoint + delay_to_column (str): The latest timepoint + """ + day_index = (self.applicant_df.time_joined_queue == time_joined_queue) & (self.applicant_df.processed == True) + delay_from_column = self.applicant_df.loc[day_index, delay_from_column] + delay_to_column = self.applicant_df.loc[day_index, delay_to_column] + return delay_to_column - delay_from_column diff --git a/household_contact_tracing/schemas/node_attributes/infection.json b/household_contact_tracing/schemas/node_attributes/infection.json new file mode 100644 index 0000000..1b2259f --- /dev/null +++ b/household_contact_tracing/schemas/node_attributes/infection.json @@ -0,0 +1,58 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "title": "Node infection attribute schema", + "description": "A schema for validating the infection attributes of nodes in a household branching process model", + "type": "object", + "anyOf": [ + {"required": ["asymptomatic"]}, + {"required": ["infecting_node"]}, + {"required": ["isolated"]}, + {"required": ["outside_house_contacts_made"]}, + {"required": ["recovered"]}, + {"required": ["recovery_time"]}, + {"required": ["spread_to_global_node_time_tuples"]}, + {"required": ["time_infected"]} + ], + "properties": { + "asymptomatic": { + "description": "Is the node asymptomatic?", + "type": "boolean" + }, + "infecting_node": { + "description": "The node that infected this node [CHECK]", + "type": ["integer", "null"] + }, + "isolated": { + "description": "Has the node isolated?", + "type": "boolean" + }, + "outside_house_contacts_made": { + "description": "Has the node made outside contacts?", + "type": "boolean" + }, + "recovered": { + "description": "Has the node recovered?", + "type": "boolean" + }, + "recovery_time": { + "description": "The time taken for the node to recover", + "type": ["number", "null"] + }, + "spread_to_global_node_time_tuples": { + "description": "List of node/times where spread to global.", + "type": "array", + "items": { + "type": "array", + "items": { + "type": "integer", + "type": "number" + } + } + }, + "time_infected": { + "description": "The time that the node was infected.", + "type": ["number", "null"] + } + } +} + diff --git a/household_contact_tracing/schemas/node_attributes/lfd_testing.json b/household_contact_tracing/schemas/node_attributes/lfd_testing.json new file mode 100644 index 0000000..980fa90 --- /dev/null +++ b/household_contact_tracing/schemas/node_attributes/lfd_testing.json @@ -0,0 +1,36 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "title": "Node infection attribute schema", + "description": "A schema for validating the attributes of nodes in a household branching process model", + "type": "object", + "anyOf": [ + {"required": ["avenue_of_testing"]}, + {"required": ["being_lateral_flow_tested"]}, + {"required": ["positive_test_time"]}, + {"required": ["taken_confirmatory_pcr_test"]}, + {"required": ["time_started_LFA_testing"]} + ], + "properties": { + "avenue_of_testing": { + "description": "What was the avenue of testing? 0=PCR, 1=LFA", + "type": ["integer", "null"] + }, + "being_lateral_flow_tested": { + "description": "Is the node being lateral flow tested?", + "type": "boolean" + }, + "positive_test_time": { + "description": "The time taken to receive a positive test result", + "type": ["number", "null"] + }, + "taken_confirmatory_pcr_test": { + "description": "Has the node taken a confirmatory PCR test?", + "type": "boolean" + }, + "time_started_LFA_testing": { + "description": "The time that the node started LFA testing", + "type": ["number", "null"] + } + } +} + diff --git a/household_contact_tracing/schemas/node_attributes/lfd_testing_adherence.json b/household_contact_tracing/schemas/node_attributes/lfd_testing_adherence.json new file mode 100644 index 0000000..7600cb2 --- /dev/null +++ b/household_contact_tracing/schemas/node_attributes/lfd_testing_adherence.json @@ -0,0 +1,31 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "title": "Node infection attribute schema", + "description": "A schema for validating the attributes of nodes in a household branching process model", + "type": "object", + "anyOf": [ + {"required": ["confirmatory_PCR_result_was_positive"]}, + {"required": ["node_will_take_up_LFA_testing"]}, + {"required": ["propensity_to_miss_lfa_tests"]}, + {"required": ["propensity_risky_behaviour_lfa_testing"]} + ], + "properties": { + "confirmatory_PCR_result_was_positive": { + "description": "Was the confirmatory PCR test positive?", + "type": "boolean" + }, + "node_will_take_up_LFA_testing": { + "description": "Will the node take up LFA testing?", + "type": "boolean" + }, + "propensity_to_miss_lfa_tests": { + "description": "Propensity of node to miss LFA tests", + "type": ["number", "null"] + }, + "propensity_risky_behaviour_lfa_testing": { + "description": "Propensity of node to engage in risky behaviour during LFA testing", + "type": ["number", "null"] + } + } +} + diff --git a/household_contact_tracing/schemas/node_attributes/returning_travellers.json b/household_contact_tracing/schemas/node_attributes/returning_travellers.json new file mode 100644 index 0000000..693bcc3 --- /dev/null +++ b/household_contact_tracing/schemas/node_attributes/returning_travellers.json @@ -0,0 +1,16 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "title": "Node returning travellers attribute schema", + "description": "A schema for validating the attributes of returning traveller nodes in a household branching process model", + "type": "object", + "anyOf": [ + {"required": ["pseudo_symptom_onset_time"]} + ], + "properties": { + "pseudo_symptom_onset_time": { + "description": "Pseudo time of onset of symptoms", + "type": ["number", "null"] + } + } +} + diff --git a/household_contact_tracing/schemas/node_attributes/tracing.json b/household_contact_tracing/schemas/node_attributes/tracing.json new file mode 100644 index 0000000..546906a --- /dev/null +++ b/household_contact_tracing/schemas/node_attributes/tracing.json @@ -0,0 +1,61 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "title": "Node infection attribute schema", + "description": "A schema for validating the contact tracing attributes of nodes in a household branching process model", + "type": "object", + "anyOf": [ + {"required": ["contact_traced"]}, + {"required": ["has_contact_tracing_app"]}, + {"required": ["propagated_contact_tracing"]}, + {"required": ["received_positive_test_result"]}, + {"required": ["received_result"]}, + {"required": ["symptom_onset_time"]}, + {"required": ["testing_delay"]}, + {"required": ["time_of_reporting"]}, + {"required": ["will_report_infection"]}, + {"required": ["completed_isolation"]} + ], + "properties": { + "contact_traced": { + "description": "Has the node been contact traced?", + "type": "boolean" + }, + "has_contact_tracing_app": { + "description": "Does this node have the contact tracing app?", + "type": "boolean" + }, + "propagated_contact_tracing": { + "description": "Has the node propagated contact tracing?", + "type": "boolean" + }, + "received_positive_test_result": { + "description": "Has the node received a positive test result?", + "type": "boolean" + }, + "received_result": { + "description": "Has the node received any test result?", + "type": "boolean" + }, + "symptom_onset_time": { + "description": "The time taken for symptom onset for this node", + "type": ["number", "null"] + }, + "testing_delay": { + "description": "The delay time before node was tested.", + "type": ["number", "null"] + }, + "time_of_reporting": { + "description": "The time that the node reported symptoms [CHECK].", + "type": ["number", "null"] + }, + "will_report_infection": { + "description": "Will the node report an infection?", + "type": "boolean" + }, + "completed_isolation": { + "description": "Has the node completed isolation?", + "type": "boolean" + } + } +} + diff --git a/household_contact_tracing/schemas/node_attributes/tracing_adherence.json b/household_contact_tracing/schemas/node_attributes/tracing_adherence.json new file mode 100644 index 0000000..2d8677a --- /dev/null +++ b/household_contact_tracing/schemas/node_attributes/tracing_adherence.json @@ -0,0 +1,20 @@ +{ + "$schema": "http://json-schema.org/draft-07/schema#", + "title": "Node tracing adherence attribute schema", + "description": "A schema for validating the tracing adherence attributes of nodes in a household branching process model", + "type": "object", + "anyOf": [ + {"required": ["will_uptake_isolation"]}, + {"required": ["propensity_imperfect_isolation"]} + ], + "properties": { + "will_uptake_isolation": { + "description": "Will this node uptake isolation?", + "type": "boolean" + }, + "propensity_imperfect_isolation": { + "description": "The propensity of this node to perform imperfect isolation?", + "type": ["boolean", "null"] + } + } +} \ No newline at end of file diff --git a/household_contact_tracing/utilities.py b/household_contact_tracing/utilities.py index 560b553..bc19fc3 100644 --- a/household_contact_tracing/utilities.py +++ b/household_contact_tracing/utilities.py @@ -1,23 +1,24 @@ -from typing import Type, List +from __future__ import annotations +from typing import Type, List, TYPE_CHECKING from copy import deepcopy import itertools -from household_contact_tracing.branching_process_models import HouseholdLevelTracing +if TYPE_CHECKING: + from household_contact_tracing.branching_process_models import HouseholdLevelTracing class ParameterError(Exception): """Raised if simulation parameters cannot be parsed.""" -def run_parameterised_simulation(model_type: Type[HouseholdLevelTracing], num_steps: int, - params: dict): +def run_parameterised_simulation(model_type: Type[HouseholdLevelTracing], num_steps: int, params: dict): """Assume sequence nesting is combinatorial at first.""" processed_params = process_sequences(params) model_results = [] for param_set in processed_params: model = model_type(param_set) - model.run_simulation(num_steps, 1000) + model.run_simulation({"max_time": num_steps, "infection_threshold": 1000}) model_results.append(model) print(len(model_results)) diff --git a/household_contact_tracing/views/statistics_view.py b/household_contact_tracing/views/statistics_view.py new file mode 100644 index 0000000..5e580e2 --- /dev/null +++ b/household_contact_tracing/views/statistics_view.py @@ -0,0 +1,250 @@ +# Code to estimate the growth rate of a simulated branching process +from household_contact_tracing.views.branching_process_view import BranchingProcessView +from household_contact_tracing.branching_process_model import BranchingProcessModel +from household_contact_tracing.branching_process_state import MaxNodesInfectiousState, ReadyState, RunningState, ExtinctState +from household_contact_tracing.exceptions import Error, ModelStateError +import scipy.stats as ss +import statsmodels.api as sm +import numpy as np + +class StatisticsView(BranchingProcessView): + + """ + View that estimates the growth rate of a completed simulation + """ + + def __init__(self, model: BranchingProcessModel): + + self._model = model + self.show = False + + def model_state_change(self, subject: BranchingProcessModel): + """ + Respond to changes in model state (e.g. running, extinct, timed-out) + + Parameters: + subject (BranchingProcessModel): The branching process model being displayed by this simulation view. + + Returns: + None + """ + # nothing to do here, usually it only makes sense to estimate the growth rate after the simulation is complete. + pass + + def model_step_increment(self, subject: BranchingProcessModel): + """ + Respond to single step increment in simulation + + Parameters: + subject (BranchingProcessModel): The branching process model being displayed by this simulation view. + + Returns: + None + """ + # nothing to do here, usually it only makes sense to estimate the growth rate after the simulation is complete. + pass + + def model_simulation_stopped(self, subject: BranchingProcessModel): + """ + Respond to end of simulation run + + Parameters: + subject (BranchingProcessModel): The branching process model being displayed by this simulation view. + + Returns: + None + """ + if self.show: + self._estimate_growth_rate() + + def graph_change(self, subject: BranchingProcessModel): + """ + Respond to changes in graph (nodes/households network) + + Parameters: + subject (SimulationModel): The branching process model being displayed by this simulation view. + + Returns: + None + """ + # nothing to do here, usually it only makes sense to estimate the growth rate after the simulation is complete. + pass + + def set_display(self, show: bool): + """ + Sets whether this view is displayed or not. + + Parameters: + show (bool): To display this view, set to True + + Returns: + None + """ + self.show = show + + def get_infection_times(self): + """ + Returns a list containing the times at which each node was infected + """ + return [node.time_infected for node in self._model.network.all_nodes()] + + def get_daily_incidence(self): + """Returns a list of the new infections at each time point. + + The list contains [time, incidence] pairs + """ + + infection_times = self.get_infection_times() + return([ + infection_times.count(t) + for t in range(self._model.time) + ]) + + def _estimate_growth_rate(self, discard_first_n_days: int = 10, verbose = True): + """Uses Poisson regression to estimate the growth rate of the epidemic. + The first few days of a simulation are typically discarded while the process becomes mixed + after it's artificial initial conditions + + Args: + discard_first_n_days (int, optional): estimate growth rate from data after the first n days. Defaults to 10. + """ + + if isinstance(self._model.state, ReadyState): + raise ModelStateError(self._model.state, 'Simulation has not started yet. Cannot estimate growth rate.') + + # we work out how many + time = self._model.time + num_eligible_dates = time - discard_first_n_days + + if num_eligible_dates < 2: + # there is not enough data to estimate the growth rate + + if isinstance(self._model.state, RunningState): + raise Error("""Cannot estimate growth rate due to insufficient eligible dates. + This simulation is still running, consider continuing the simulation before estimating the growth rate.""") + + elif isinstance(self._model.state, ExtinctState): + raise Error("""Cannot estimate growth rate due to insufficient eligible dates. + This simulation went extinct, possibly before discard_first_n_days. Consider starting the simulation with more infections""") + + elif isinstance(self._model.state, MaxNodesInfectiousState): + raise Error("""Cannot estimate growth rate due to insufficient eligible dates. + This simulation exceeded the maximum number of infectious nodes. Consider raising max_active_infections when simulating.""") + + else: + raise Error("""Cannot estimate growth rate due to insufficient eligible dates.""") + + else: + # there is some data that can be used to estimate the growth rate. Perform analysis + + if verbose: + print(f'Estimating growth rate using {num_eligible_dates} time periods') + + # the incidence after the first n days + y = self.get_daily_incidence()[discard_first_n_days:] + + # create a simple design matrix + X = [[t] for t in range(num_eligible_dates)] + X = sm.add_constant(X, prepend=False) + + glm_poisson = sm.GLM(y, X, family=sm.families.Poisson()) + self.glm_poisson = glm_poisson.fit() + + def get_growth_rate(self, discard_first_n_days: int = 10, verbose: bool = True): + """Returns the growth rate of the simulated epidemic, estimated using poisson regression. + + The first few days of a simulation are typically discarded while the process becomes mixed + after it's artificial initial conditions. + + Args: + discard_first_n_days (int, optional): estimate growth rate from data after the first n days. Defaults to 10. + """ + self._estimate_growth_rate(discard_first_n_days, verbose) + + return self.glm_poisson.params[0] + + def growth_rate_summary(self, discard_first_n_days: int = 10, alpha: float = 0.05, glm_summary: bool = False): + """Returns the growth rate of the simulated epidemic, estimated using poisson regression. + + The first few days of a simulation are typically discarded while the process becomes mixed + after it's artificial initial conditions. + + Args: + discard_first_n_days (int, optional): estimate growth rate from data after the first n days. Defaults to 10. + """ + + self._estimate_growth_rate(discard_first_n_days, verbose = False) + + num_eligible_dates = self._model.time - discard_first_n_days + growth_rate = self.glm_poisson.params[0] + growth_rate_ci = self.glm_poisson.conf_int(alpha = alpha, cols = [0])[0] + doubling_time = np.log(2) / np.log(1 + growth_rate) + doubling_time_ci = np.log(2) / np.log(1 + np.array(growth_rate_ci)) + + if glm_summary: + print(self.glm_poisson.summary()) + print(f'{num_eligible_dates} time periods were used to estimate the growth rate.') + print(f'The estimated growth rate was {round(growth_rate*100, 2)}% ({100*(1-alpha)}% CI: {round(growth_rate_ci[0]*100,2)}-{round(growth_rate_ci[1]*100,2)}%) per day.') + print(f'The estimated doubling time is {round(doubling_time, 2)} ({100*(1-alpha)}% CI: {round(doubling_time_ci[1],2)}-{round(doubling_time_ci[0],2)}) days.') + + def _estimate_household_secondary_attack_rate(self, use_first_generation_only: bool = False) -> None: + if isinstance(self._model.state, ReadyState): + raise ModelStateError(self._model.state, 'Simulation has not started yet. Cannot estimate growth rate.') + + if use_first_generation_only: + households_with_completed_local_epidemics = [ + household + for household + in self._model.network.all_households + if household.local_epidemic_completed + and household.id in self._model.infection.starting_households + ] + else: + households_with_completed_local_epidemics = [ + household + for household + in self._model.network.all_households + if household.local_epidemic_completed + ] + + # size of household - number of remaining susceptibles = final size. + # we subtract 1, to work out the number of non-index secondary infections + self.total_infected = sum([ + household.size - household.susceptibles - 1 + for household + in households_with_completed_local_epidemics + ]) + + # we subtract 1, to work out the number of non-index exposed individuals + self.total_exposed = sum([ + household.size - 1 for household in households_with_completed_local_epidemics + ]) + + self.n_households_with_completed_local_epidemics = len(households_with_completed_local_epidemics) + + self.household_sar = self.total_infected / self.total_exposed + + # calculating some confidence intervals using the good ol' Jefferys interval + self.household_sar_ci = ss.beta.interval(alpha = 0.95, a = self.total_infected + 0.5, b = self.total_exposed - self.total_infected + 0.5) + + def get_hh_sar(self): + self._estimate_household_secondary_attack_rate(use_first_generation_only=True) + return self.household_sar + + def household_secondary_attack_rate_summary(self, use_first_generation_only: bool = False, alpha: float = 0.95) -> None: + """Estimates the household secondary attack rate, and prints and interpretable output. + + Args: + use_first_generation_only (bool, optional): Use the first generation of the household epidemic only to estimate the household secondary attack. Defaults to False. + """ + self._estimate_household_secondary_attack_rate(use_first_generation_only) + + + print('Household secondary attack rate summary:') + print(f'{self.n_households_with_completed_local_epidemics} local household epidemics were eligible to be included.') + if use_first_generation_only: + print('Only the first generation of the household epidemic was included in this calculation.') + else: + print('All households with completed local epidemics were included. This may lead to a biased sample, as it is possible that local epidemics with a long duration were not included.') + print(f'There were {self.total_exposed} non-index susceptible individuals exposed, of which {self.total_infected} were infected.') + print(f'This yields a household secondary attack rate of {round(self.household_sar*100)}% ({int(alpha * 100)}% CI: {round(self.household_sar_ci[0]*100)}-{round(self.household_sar_ci[1]*100)}%).') diff --git a/requirements.txt b/requirements.txt index f01fa6b..ade7b51 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,4 +10,7 @@ loguru pyyaml pandas beautifulsoup4==4.9.3 -pyvis==0.1.9 \ No newline at end of file +pyvis==0.1.9 +statsmodels +torch +ax diff --git a/test/fixtures/queueing_processes/DeterministicQueue_add_new_test_seekers.pickle b/test/fixtures/queueing_processes/DeterministicQueue_add_new_test_seekers.pickle new file mode 100644 index 0000000..77307c3 Binary files /dev/null and b/test/fixtures/queueing_processes/DeterministicQueue_add_new_test_seekers.pickle differ diff --git a/test/fixtures/queueing_processes/Queue_new_applicants_fixture.pickle b/test/fixtures/queueing_processes/Queue_new_applicants_fixture.pickle new file mode 100644 index 0000000..83dff41 Binary files /dev/null and b/test/fixtures/queueing_processes/Queue_new_applicants_fixture.pickle differ diff --git a/test/fixtures/queueing_processes/Queue_swab_applicants.pickle b/test/fixtures/queueing_processes/Queue_swab_applicants.pickle new file mode 100644 index 0000000..c58c3fe Binary files /dev/null and b/test/fixtures/queueing_processes/Queue_swab_applicants.pickle differ diff --git a/test/fixtures/queueing_processes/empty_applicant_df.pickle b/test/fixtures/queueing_processes/empty_applicant_df.pickle new file mode 100644 index 0000000..a25a3f0 Binary files /dev/null and b/test/fixtures/queueing_processes/empty_applicant_df.pickle differ diff --git a/test/fixtures/queueing_processes/empty_queue_df.pickle b/test/fixtures/queueing_processes/empty_queue_df.pickle new file mode 100644 index 0000000..0dfae63 Binary files /dev/null and b/test/fixtures/queueing_processes/empty_queue_df.pickle differ diff --git a/test/fixtures/queueing_processes/update_fixtures.ipynb b/test/fixtures/queueing_processes/update_fixtures.ipynb new file mode 100644 index 0000000..049c738 --- /dev/null +++ b/test/fixtures/queueing_processes/update_fixtures.ipynb @@ -0,0 +1,1696 @@ +{ + "metadata": { + "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.9.4" + }, + "orig_nbformat": 4, + "kernelspec": { + "name": "python3", + "display_name": "Python 3.9.4 64-bit (conda)" + }, + "interpreter": { + "hash": "5e088ab515b9ff01afc6092114ac7786286eca1d23e3587660ba5a8d6e06cf28" + } + }, + "nbformat": 4, + "nbformat_minor": 2, + "cells": [ + { + "cell_type": "code", + "execution_count": 8, + "source": [ + "from household_contact_tracing.queueing_processes import Queue, DeterministicQueue\r\n", + "import pandas as pd\r\n", + "import pickle" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "This notebook is used to create the fixtures for test_queueing_processes.py. We save fixtures using pickle protocol 4 for backwards compatibility." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Queue" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 9, + "source": [ + "my_queue = Queue(days_to_simulate=10, capacity=[10]*10)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### Initialised applicant df" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 10, + "source": [ + "my_queue.applicant_df" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "Empty DataFrame\n", + "Columns: [id, processed, waiting_to_be_processed, left_queue_not_processed, time_symptom_onset, time_joined_queue, time_processed, time_received_result, time_will_leave_queue]\n", + "Index: []" + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queue
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 10 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 11, + "source": [ + "pd.read_pickle('./test/fixtures/queueing_processes/empty_applicant_df.pickle')" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "Empty DataFrame\n", + "Columns: [id, processed, waiting_to_be_processed, left_queue_not_processed, time_symptom_onset, time_joined_queue, time_processed, time_received_result, time_will_leave_queue]\n", + "Index: []" + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queue
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 11 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 12, + "source": [ + "my_queue.applicant_df.to_pickle('./test/fixtures/queueing_processes/empty_applicant_df.pickle', protocol=4)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### Initialised queue_df" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 13, + "source": [ + "my_queue.queue_df" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " time capacity new_applicants spillover_to_next_day \\\n", + "0 0 10 \n", + "1 1 10 \n", + "2 2 10 \n", + "3 3 10 \n", + "4 4 10 \n", + "5 5 10 \n", + "6 6 10 \n", + "7 7 10 \n", + "8 8 10 \n", + "9 9 10 \n", + "\n", + " total_applications_today capacity_exceeded capacity_exceeded_by \\\n", + "0 \n", + "1 \n", + "2 \n", + "3 \n", + "4 \n", + "5 \n", + "6 \n", + "7 \n", + "8 \n", + "9 \n", + "\n", + " number_processed_today number_left_queue_not_tested \n", + "0 \n", + "1 \n", + "2 \n", + "3 \n", + "4 \n", + "5 \n", + "6 \n", + "7 \n", + "8 \n", + "9 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
timecapacitynew_applicantsspillover_to_next_daytotal_applications_todaycapacity_exceededcapacity_exceeded_bynumber_processed_todaynumber_left_queue_not_tested
0010
1110
2210
3310
4410
5510
6610
7710
8810
9910
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 13 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 14, + "source": [ + "pd.read_pickle('./test/fixtures/queueing_processes/empty_queue_df.pickle')" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " time capacity new_applicants spillover_to_next_day \\\n", + "0 0 10 \n", + "1 1 10 \n", + "2 2 10 \n", + "3 3 10 \n", + "4 4 10 \n", + "5 5 10 \n", + "6 6 10 \n", + "7 7 10 \n", + "8 8 10 \n", + "9 9 10 \n", + "\n", + " total_applications_today capacity_exceeded capacity_exceeded_by \\\n", + "0 \n", + "1 \n", + "2 \n", + "3 \n", + "4 \n", + "5 \n", + "6 \n", + "7 \n", + "8 \n", + "9 \n", + "\n", + " number_processed_today number_left_queue_not_tested \n", + "0 \n", + "1 \n", + "2 \n", + "3 \n", + "4 \n", + "5 \n", + "6 \n", + "7 \n", + "8 \n", + "9 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
timecapacitynew_applicantsspillover_to_next_daytotal_applications_todaycapacity_exceededcapacity_exceeded_bynumber_processed_todaynumber_left_queue_not_tested
0010
1110
2210
3310
4410
5510
6610
7710
8810
9910
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 14 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 15, + "source": [ + "my_queue.queue_df.to_pickle('./test/fixtures/queueing_processes/empty_queue_df.pickle', protocol=4)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### add_new_applicants" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 16, + "source": [ + "my_queue.add_new_applicants(\r\n", + " ids = ['A', 'B', 'C'],\r\n", + " time = [1,2,3],\r\n", + " symptom_onset_times = [6, 6, 6],\r\n", + " queue_leaving_times=[11,12,13])" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 17, + "source": [ + "my_queue.applicant_df" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id processed waiting_to_be_processed left_queue_not_processed \\\n", + "0 A False True \n", + "1 B False True \n", + "2 C False True \n", + "\n", + " time_symptom_onset time_joined_queue time_processed time_received_result \\\n", + "0 6 1 \n", + "1 6 2 \n", + "2 6 3 \n", + "\n", + " time_will_leave_queue \n", + "0 11 \n", + "1 12 \n", + "2 13 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queue
0AFalseTrue6111
1BFalseTrue6212
2CFalseTrue6313
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 17 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 18, + "source": [ + "pd.read_pickle('./test/fixtures/queueing_processes/Queue_new_applicants_fixture.pickle')" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id processed waiting_to_be_processed left_queue_not_processed \\\n", + "0 A False True \n", + "1 B False True \n", + "2 C False True \n", + "\n", + " time_symptom_onset time_joined_queue time_processed time_received_result \\\n", + "0 6 1 \n", + "1 6 2 \n", + "2 6 3 \n", + "\n", + " time_will_leave_queue \n", + "0 11 \n", + "1 12 \n", + "2 13 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queue
0AFalseTrue6111
1BFalseTrue6212
2CFalseTrue6313
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 18 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 20, + "source": [ + "# my_queue.applicant_df.to_pickle('./test/fixtures/queueing_processes/Queue_new_applicants_fixture.pickle', protocol=4)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### swab_applicants" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 21, + "source": [ + "my_queue = Queue(days_to_simulate=10, capacity=[10]*10)\r\n", + "my_queue.add_new_applicants(\r\n", + " ids = ['A', 'B', 'C'],\r\n", + " time = [1,2,3],\r\n", + " symptom_onset_times = [6, 6, 6],\r\n", + " queue_leaving_times=[11,12,13])\r\n", + "my_queue.swab_applicants(\r\n", + " to_be_processed=[1, 2],\r\n", + " processing_delays=[1, 2])" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 22, + "source": [ + "my_queue.applicant_df" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id processed waiting_to_be_processed left_queue_not_processed \\\n", + "0 A False True \n", + "1 B True False False \n", + "2 C True False False \n", + "\n", + " time_symptom_onset time_joined_queue time_processed time_received_result \\\n", + "0 6 1 \n", + "1 6 2 0 1 \n", + "2 6 3 0 2 \n", + "\n", + " time_will_leave_queue \n", + "0 11 \n", + "1 12 \n", + "2 13 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queue
0AFalseTrue6111
1BTrueFalseFalse620112
2CTrueFalseFalse630213
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 22 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 23, + "source": [ + "pd.read_pickle('./test/fixtures/queueing_processes/Queue_swab_applicants.pickle')" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id processed waiting_to_be_processed left_queue_not_processed \\\n", + "0 A False True \n", + "1 B True False False \n", + "2 C True False False \n", + "\n", + " time_symptom_onset time_joined_queue time_processed time_received_result \\\n", + "0 6 1 \n", + "1 6 2 0 1 \n", + "2 6 3 0 2 \n", + "\n", + " time_will_leave_queue \n", + "0 11 \n", + "1 12 \n", + "2 13 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queue
0AFalseTrue6111
1BTrueFalseFalse620112
2CTrueFalseFalse630213
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 23 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 24, + "source": [ + "my_queue.applicant_df.to_pickle('./test/fixtures/queueing_processes/Queue_swab_applicants.pickle', protocol=4)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Deterministic Queue" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### Add new test seekers" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 25, + "source": [ + "def processing_delay_dist():\r\n", + " return 1\r\n", + "\r\n", + "def symptom_onset_delay_dist():\r\n", + " return 2 \r\n", + "\r\n", + "my_det_queue = DeterministicQueue(\r\n", + " days_to_simulate = 10,\r\n", + " demand = [10]*10,\r\n", + " capacity = [10]*10,\r\n", + " max_time_in_queue = 10,\r\n", + " processing_delay_dist = processing_delay_dist,\r\n", + " symptom_onset_delay_dist = symptom_onset_delay_dist,\r\n", + " selection_method = 'uniform'\r\n", + ")" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 26, + "source": [ + "my_det_queue.add_new_test_seekers()" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 27, + "source": [ + "my_det_queue.queue.applicant_df" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id processed waiting_to_be_processed left_queue_not_processed \\\n", + "0 False True \n", + "1 False True \n", + "2 False True \n", + "3 False True \n", + "4 False True \n", + "5 False True \n", + "6 False True \n", + "7 False True \n", + "8 False True \n", + "9 False True \n", + "\n", + " time_symptom_onset time_joined_queue time_processed time_received_result \\\n", + "0 -2 0 \n", + "1 -2 0 \n", + "2 -2 0 \n", + "3 -2 0 \n", + "4 -2 0 \n", + "5 -2 0 \n", + "6 -2 0 \n", + "7 -2 0 \n", + "8 -2 0 \n", + "9 -2 0 \n", + "\n", + " time_will_leave_queue \n", + "0 8 \n", + "1 8 \n", + "2 8 \n", + "3 8 \n", + "4 8 \n", + "5 8 \n", + "6 8 \n", + "7 8 \n", + "8 8 \n", + "9 8 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queue
0FalseTrue-208
1FalseTrue-208
2FalseTrue-208
3FalseTrue-208
4FalseTrue-208
5FalseTrue-208
6FalseTrue-208
7FalseTrue-208
8FalseTrue-208
9FalseTrue-208
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 27 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 28, + "source": [ + "my_det_queue.queue.applicant_df.to_pickle('./test/fixtures/queueing_processes/DeterministicQueue_add_new_test_seekers.pickle', protocol=4)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### test_process_queue_excess_capacity" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 10, + "source": [ + "def processing_delay_dist():\r\n", + " return 1\r\n", + "\r\n", + "def symptom_onset_delay_dist():\r\n", + " return 2 \r\n", + "\r\n", + "my_det_queue = DeterministicQueue(\r\n", + " days_to_simulate = 10,\r\n", + " demand = [5]*10,\r\n", + " capacity = [10]*10,\r\n", + " max_time_in_queue = 10,\r\n", + " processing_delay_dist = processing_delay_dist,\r\n", + " symptom_onset_delay_dist = symptom_onset_delay_dist\r\n", + ")\r\n", + "\r\n", + "my_det_queue.add_new_test_seekers()\r\n", + "\r\n" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 11, + "source": [ + "my_det_queue.queue.applicant_df" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id processed waiting_to_be_processed left_queue_not_processed time_symptom_onset \\\n", + "0 5 False True \n", + "1 5 False True \n", + "2 5 False True \n", + "3 5 False True \n", + "4 5 False True \n", + "\n", + " time_joined_queue time_processed time_received_result time_will_leave_queue \\\n", + "0 12 \n", + "1 12 \n", + "2 12 \n", + "3 12 \n", + "4 12 \n", + "\n", + " symptom_onset \n", + "0 2.0 \n", + "1 2.0 \n", + "2 2.0 \n", + "3 2.0 \n", + "4 2.0 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queuesymptom_onset
05FalseTrue122.0
15FalseTrue122.0
25FalseTrue122.0
35FalseTrue122.0
45FalseTrue122.0
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 11 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 12, + "source": [ + "my_det_queue.process_queue()" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 13, + "source": [ + "my_det_queue.queue.applicant_df" + ], + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + " id processed waiting_to_be_processed left_queue_not_processed time_symptom_onset \\\n", + "0 5 True False False \n", + "1 5 True False False \n", + "2 5 True False False \n", + "3 5 True False False \n", + "4 5 True False False \n", + "\n", + " time_joined_queue time_processed time_received_result time_will_leave_queue \\\n", + "0 0 1 12 \n", + "1 0 1 12 \n", + "2 0 1 12 \n", + "3 0 1 12 \n", + "4 0 1 12 \n", + "\n", + " symptom_onset \n", + "0 2.0 \n", + "1 2.0 \n", + "2 2.0 \n", + "3 2.0 \n", + "4 2.0 " + ], + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
idprocessedwaiting_to_be_processedleft_queue_not_processedtime_symptom_onsettime_joined_queuetime_processedtime_received_resulttime_will_leave_queuesymptom_onset
05TrueFalseFalse01122.0
15TrueFalseFalse01122.0
25TrueFalseFalse01122.0
35TrueFalseFalse01122.0
45TrueFalseFalse01122.0
\n", + "
" + ] + }, + "metadata": {}, + "execution_count": 13 + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, + "source": [], + "outputs": [], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/test/test_BranchingProcessSimulation.py b/test/test_BranchingProcessSimulation.py index c1bcf6f..4efedda 100644 --- a/test/test_BranchingProcessSimulation.py +++ b/test/test_BranchingProcessSimulation.py @@ -25,10 +25,10 @@ def test_asymptomatic_nodes_attributes(): lfa_test_node = test_model.network.node(1) - assert lfa_test_node.asymptomatic is True + assert lfa_test_node.infection.asymptomatic is True # Symptom onset time is infinite - assert lfa_test_node.symptom_onset_time > 10000 - assert lfa_test_node.will_report_infection is False + assert lfa_test_node.tracing.symptom_onset_time > 10000 + assert lfa_test_node.tracing.will_report_infection is False def test_symptomatic_nodes_attributes(): @@ -42,9 +42,9 @@ def test_symptomatic_nodes_attributes(): lfa_test_node = test_model.network.node(1) - assert lfa_test_node.asymptomatic is False - assert lfa_test_node.symptom_onset_time == 5 - assert lfa_test_node.will_report_infection is True + assert lfa_test_node.infection.asymptomatic is False + assert lfa_test_node.tracing.symptom_onset_time == 5 + assert lfa_test_node.tracing.will_report_infection is True @pytest.fixture diff --git a/test/test_TestingContactModel.py b/test/test_TestingContactModel.py index f3cca6b..793a8a6 100644 --- a/test/test_TestingContactModel.py +++ b/test/test_TestingContactModel.py @@ -127,13 +127,13 @@ def prob_testing_positive_pcr_func(infectious_age): model.prob_pcr_positive = prob_testing_positive_pcr_func model.prob_lfa_positive = prob_testing_positive_lfa_func - assert model.network.node(1).pseudo_symptom_onset_time == 5 + assert model.network.node(1).returning_travellers.pseudo_symptom_onset_time == 5 def test_pseudo_symptom_onset(simple_model): """Checks that it is also working for symptomatics """ - assert simple_model.network.node(1).pseudo_symptom_onset_time == 5 + assert simple_model.network.node(1).returning_travellers.pseudo_symptom_onset_time == 5 def test_time_relative_to_symptom_onset(simple_model): @@ -165,7 +165,7 @@ def test_being_lateral_flow_tested_attribute(simple_model): """Check nodes are generated with the lateral flow testing attribute """ - assert not simple_model.network.node(1).being_lateral_flow_tested + assert not simple_model.network.node(1).lfd_testing.being_lateral_flow_tested def test_get_positive_lateral_flow_nodes_default_exclusion(simple_model_high_test_prob): @@ -185,7 +185,7 @@ def test_get_positive_lateral_flow_nodes_timings(simple_model_high_test_prob): node_of_interest = model.network.node(1) - node_of_interest.being_lateral_flow_tested = True + node_of_interest.lfd_testing.being_lateral_flow_tested = True assert model.intervention.lft_nodes(model.time, model.prob_lfa_positive) == [] @@ -198,7 +198,7 @@ def test_get_positive_lateral_flow_nodes(simple_model_high_test_prob): node_of_interest = model.network.node(1) - node_of_interest.being_lateral_flow_tested = True + node_of_interest.lfd_testing.being_lateral_flow_tested = True model.time = 5 @@ -236,7 +236,7 @@ def prob_testing_positive_pcr_func(infectious_age): model.infection.new_outside_household_infection(time=0, infecting_node=model.network.node(1)) - model.intervention.increment_tracing.attempt_contact_trace_of_household( + model.intervention.increment_tracing.contact_trace_household( house_to=model.network.household(2), house_from=model.network.household(1), days_since_contact_occurred=0, @@ -247,7 +247,7 @@ def prob_testing_positive_pcr_func(infectious_age): model.simulate_one_step() model.simulate_one_step() - assert model.network.node(2).being_lateral_flow_tested is True + assert model.network.node(2).lfd_testing.being_lateral_flow_tested is True def test_isolate_positive_lateral_flow_tests(simple_model_high_test_prob: IndividualTracingDailyTesting): @@ -257,7 +257,7 @@ def test_isolate_positive_lateral_flow_tests(simple_model_high_test_prob: Indivi model.time = 5 - model.network.node(1).being_lateral_flow_tested = True + model.network.node(1).lfd_testing.being_lateral_flow_tested = True positive_nodes = model.intervention.lft_nodes(model.time, model.prob_lfa_positive) new_isolation = DailyTestingIsolation(model.network, model.params) @@ -267,11 +267,11 @@ def test_isolate_positive_lateral_flow_tests(simple_model_high_test_prob: Indivi # but they are lfa testing model.infection.new_within_household_infection(time=model.time, infecting_node=model.network.node(1)) - assert model.network.node(1).isolated + assert model.network.node(1).infection.isolated assert model.network.household(1).applied_household_positive_policy - assert model.network.node(1).received_positive_test_result - assert not model.network.node(2).isolated - assert model.network.node(2).being_lateral_flow_tested + assert model.network.node(1).tracing.received_positive_test_result + assert not model.network.node(2).infection.isolated + assert model.network.node(2).lfd_testing.being_lateral_flow_tested @pytest.fixture @@ -315,7 +315,7 @@ def test_start_lateral_flow_testing_household_and_quarantine( model.time = 5 - model.network.node(1).being_lateral_flow_tested = True + model.network.node(1).lfd_testing.being_lateral_flow_tested = True positive_nodes = model.intervention.lft_nodes(model.time, model.prob_lfa_positive) @@ -325,11 +325,11 @@ def test_start_lateral_flow_testing_household_and_quarantine( model.infection.new_within_household_infection(time=model.time, infecting_node=model.network.node(1)) - assert model.network.node(1).isolated + assert model.network.node(1).infection.isolated assert model.network.household(1).applied_household_positive_policy - assert model.network.node(1).received_positive_test_result - assert model.network.node(2).isolated - assert model.network.node(2).being_lateral_flow_tested + assert model.network.node(1).tracing.received_positive_test_result + assert model.network.node(2).infection.isolated + assert model.network.node(2).lfd_testing.being_lateral_flow_tested @pytest.fixture @@ -375,7 +375,7 @@ def test_household_contacts_quarantine_only( model.time = 5 - model.network.node(1).being_lateral_flow_tested = True + model.network.node(1).lfd_testing.being_lateral_flow_tested = True # this line is required before the isolate_positive_lateral_flow_tests func can work positive_nodes = model.intervention.lft_nodes(model.time, model.prob_lfa_positive) @@ -385,18 +385,18 @@ def test_household_contacts_quarantine_only( model.infection.new_within_household_infection(time=model.time, infecting_node=model.network.node(1)) - assert model.network.node(1).isolated + assert model.network.node(1).infection.isolated assert model.network.household(1).applied_household_positive_policy - assert model.network.node(1).received_positive_test_result - assert model.network.node(2).isolated - assert model.network.node(2).being_lateral_flow_tested + assert model.network.node(1).tracing.received_positive_test_result + assert model.network.node(2).infection.isolated + assert model.network.node(2).lfd_testing.being_lateral_flow_tested def test_risky_behaviour_attributes_default(simple_model: simple_model): """Tests that the default behaviour is no more risky behaviour """ - assert not simple_model.network.node(1).propensity_risky_behaviour_lfa_testing + assert not simple_model.network.node(1).lfd_testing_adherence.propensity_risky_behaviour_lfa_testing def test_risky_behaviour_attributes(simple_model_risky_behaviour: simple_model_risky_behaviour): @@ -404,7 +404,7 @@ def test_risky_behaviour_attributes(simple_model_risky_behaviour: simple_model_r being tested. """ - assert simple_model_risky_behaviour.network.node(1).propensity_risky_behaviour_lfa_testing + assert simple_model_risky_behaviour.network.node(1).lfd_testing_adherence.propensity_risky_behaviour_lfa_testing @pytest.fixture @@ -458,7 +458,7 @@ def test_lfa_tested_nodes_make_more_contacts_if_risky( model = simple_model_risky_behaviour_2_infections - model.network.node(1).propensity_risky_behaviour_lfa_testing = False + model.network.node(1).lfd_testing_adherence.propensity_risky_behaviour_lfa_testing = False # stop there being any within household infections # not sure if this is strictly necessary @@ -468,12 +468,12 @@ def test_lfa_tested_nodes_make_more_contacts_if_risky( model.network.household(2).susceptibles = 0 # set the nodes to being lfa tested - model.network.node(1).being_lateral_flow_tested = True - model.network.node(2).being_lateral_flow_tested = True + model.network.node(1).lfd_testing.being_lateral_flow_tested = True + model.network.node(2).lfd_testing.being_lateral_flow_tested = True for _ in range(5): model.simulate_one_step() # node 1 does not engage in risky behaviour and should not make any global contacts - assert model.network.node(1).outside_house_contacts_made == 0 - assert model.network.node(2).outside_house_contacts_made != 0 + assert model.network.node(1).infection.outside_house_contacts_made == 0 + assert model.network.node(2).infection.outside_house_contacts_made != 0 diff --git a/test/test_integration.py b/test/test_integration.py index f4cdc6f..e8dd6f9 100644 --- a/test/test_integration.py +++ b/test/test_integration.py @@ -51,7 +51,7 @@ def run_simulation(params: dict, days=10) -> BranchingProcessModel: model.""" controller = BranchingProcessController(bpm.HouseholdLevelTracing(params)) controller.set_graphic_displays(False) - controller.run_simulation(days) + controller.run_simulation({"max_time": days}) return controller.model @@ -61,7 +61,7 @@ def nodes_isolating_correctly(network: Network) -> List[bool]: isolating_correctly = [] for node in network.all_nodes(): if node.household.isolated: - if node.isolated: + if node.infection.isolated: isolating_correctly.append(True) else: isolating_correctly.append(False) @@ -171,7 +171,7 @@ def test_simple_testing(self, household_params): numpy.random.seed(42) network = self.run_simulation(household_params).network - assert network.node(1).testing_delay != 0 + assert network.node(1).tracing.testing_delay != 0 # No intervention should expire by day 10 so all whose household is isolated should # be isolating. @@ -242,7 +242,7 @@ def test_imperfect_isolation(self, household_params): numpy.random.seed(42) model = self.run_simulation(household_params) network = model.network - node_imperfect = [node.propensity_imperfect_isolation for node in network.all_nodes()] + node_imperfect = [node.tracing_adherence.propensity_imperfect_isolation for node in network.all_nodes()] assert any(node_imperfect) node_contact_rate_reduction = \ [model.infection.contact_rate_reduction.get_contact_rate_reduction(node) for node in network.all_nodes()] @@ -300,7 +300,7 @@ def run_simulation(params: dict, days=10) -> BranchingProcessModel: model.""" controller = BranchingProcessController(bpm.IndividualLevelTracing(params)) controller.set_graphic_displays(False) - controller.run_simulation(days) + controller.run_simulation({"max_time": days}) return controller.model @@ -350,7 +350,7 @@ def run_simulation(params: dict, days=10) -> BranchingProcessModel: return the model.""" controller = BranchingProcessController(bpm.IndividualTracingDailyTesting(params)) controller.set_graphic_displays(False) - controller.run_simulation(days) + controller.run_simulation({"max_time": days}) return controller.model diff --git a/test/test_queueing_processes.py b/test/test_queueing_processes.py new file mode 100644 index 0000000..44d8f4d --- /dev/null +++ b/test/test_queueing_processes.py @@ -0,0 +1,142 @@ +import pandas as pd +import pytest +from household_contact_tracing.queueing_processes import Queue, DeterministicQueue + +# Testing the Queue object + +@pytest.fixture +def simple_Queue(): + """ + Creates a Queue object with some easy parameters for testing. + """ + + my_queue = Queue(days_to_simulate=10, capacity=[10]*10) + + return my_queue + + +@pytest.fixture +def empty_applicant_df_fixture(): + """ + Loads a fixture of an empty, correctly initialised applicant_df + """ + return pd.read_pickle('./test/fixtures/queueing_processes/empty_applicant_df.pickle') + + +def test_Queue_init_applicant_df(simple_Queue, empty_applicant_df_fixture): + return pd.testing.assert_frame_equal(empty_applicant_df_fixture, simple_Queue.applicant_df) + + +@pytest.fixture +def empty_queue_df_fixture(): + """ + Loads a fixture of an empty, correctly initialised applicant_df + """ + return pd.read_pickle('./test/fixtures/queueing_processes/empty_queue_df.pickle') + + +def test_Queue_init_queue_df(simple_Queue, empty_queue_df_fixture): + return pd.testing.assert_frame_equal(empty_queue_df_fixture, simple_Queue.queue_df) + + +@pytest.fixture +def Queue_new_applicants_fixture(): + """ + Load a fixture of an applicant df where several applicants have been added + """ + return pd.read_pickle('./test/fixtures/queueing_processes/Queue_new_applicants_fixture.pickle') + + +def test_Queue_new_applicants(simple_Queue, Queue_new_applicants_fixture): + simple_Queue.add_new_applicants( + ids = ['A', 'B', 'C'], + time = [1,2,3], + symptom_onset_times = [6, 6, 6], + queue_leaving_times=[11,12,13]) + + return pd.testing.assert_frame_equal(Queue_new_applicants_fixture, simple_Queue.applicant_df) + + +@pytest.fixture +def Queue_swab_applicants_fixture(): + """Loads a fixture where some applicants have been processed. + """ + return pd.read_pickle('./test/fixtures/queueing_processes/Queue_swab_applicants.pickle') + +def test_Queue_swab_applicants(simple_Queue, Queue_swab_applicants_fixture): + """Adds some applicants, processes some of the and checks the applicant df + """ + simple_Queue.add_new_applicants( + ids = ['A', 'B', 'C'], + time = [1,2,3], + symptom_onset_times = [6, 6, 6], + queue_leaving_times=[11,12,13]) + simple_Queue.swab_applicants( + to_be_processed=[1, 2], + processing_delays=[1, 2]) + + return pd.testing.assert_frame_equal(Queue_swab_applicants_fixture, simple_Queue.applicant_df) + + +def test_Queue_current_applicants(simple_Queue): + """Checks that the waiting to be processed indexes are returned. + + Add 3 people to the queue, process 2 + """ + simple_Queue.add_new_applicants( + ids = ['A', 'B', 'C'], + time = [1,2,3], + symptom_onset_times = [6, 6, 6], + queue_leaving_times=[11,12,13]) + simple_Queue.swab_applicants( + to_be_processed=[1, 2], + processing_delays=[1, 2]) + assert simple_Queue.current_applicants == [0] + + +def test_Queue_todays_capacity(): + """Checks that the queue returns the right value for todays capacity. + """ + queue = Queue(days_to_simulate=10, capacity=list(range(10))) + + queue.time = 4 + + assert queue.todays_capacity == 4 + + +@pytest.fixture +def DeterministicQueue_add_new_test_seekers_fixture(): + return pd.read_pickle('./test/fixtures/queueing_processes/DeterministicQueue_add_new_test_seekers.pickle') + + +def test_DeterministicQueue_add_new_test_seekers(DeterministicQueue_add_new_test_seekers_fixture): + """ + Checks that the add_new_test_seekers method correctly modifies the dataframe + by adding test seekers based upon the demand + """ + def processing_delay_dist(): + return 1 + + def symptom_onset_delay_dist(): + return 2 + + my_det_queue = DeterministicQueue( + days_to_simulate = 10, + demand = [10]*10, + capacity = [10]*10, + max_time_in_queue = 10, + processing_delay_dist = processing_delay_dist, + symptom_onset_delay_dist = symptom_onset_delay_dist, + selection_method = 'uniform' + ) + + my_det_queue.add_new_test_seekers() + + return pd.testing.assert_frame_equal(DeterministicQueue_add_new_test_seekers_fixture, my_det_queue.queue.applicant_df) + + +def test_process_queue_excess_capacity(): + pass + +def test_process_queue_excess_demand(): + pass