From df40d709782c4ee6fb2e709c53a7cee1cf0dd51c Mon Sep 17 00:00:00 2001 From: Andrejh Date: Mon, 15 Mar 2021 23:00:32 +0100 Subject: [PATCH 1/7] improve labels of plots --- cara/montecarlo.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 6ee1af99..912f3d4c 100644 --- a/cara/montecarlo.py +++ b/cara/montecarlo.py @@ -11,7 +11,7 @@ import matplotlib.pyplot as plt import matplotlib.patches as patches import matplotlib.lines as mlines from sklearn.neighbors import KernelDensity -TIME_STEP = 0.01 +TIME_STEP = 0.001 USE_SCOEH = False @@ -688,7 +688,7 @@ def composite_plot_pi_vs_viral_load(baselines: typing.List[MCExposureModel], lab axs[1, 0].set_xticklabels(['$10^{' + str(i) + '}$' for i in range(2, 13, 2)]) axs[1, 0].set_xlim(2, 12) axs[1, 0].set_xlabel('Viral load (RNA copies mL$^{-1}$)', fontsize=12) - axs[0, 0].set_ylabel('Probability of infection\n$P(I|qID=60)$', fontsize=12) + axs[0, 0].set_ylabel('Probability of infection\n$P(\,\mathtt{I}\,|\,\mathrm{vl}\,)$', fontsize=12) plt.suptitle(title, fontsize=12) axs[0, 0].text(11, -0.01, '$(i)$') @@ -827,8 +827,8 @@ def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid: """ fig, axs = plt.subplots(3, 1, sharex='all') - plt.suptitle("$F(qR|qID=$" + str(qid) + "$)$",fontsize=14, y=0.93) - fig.text(0.02, 0.5, 'Cumulative Distribution Function', va='center', rotation='vertical',fontsize=14) + plt.suptitle("$F(\mathrm{qR}|\mathrm{qID}=$" + str(qid) + "$)$",fontsize=16, y=0.93) + fig.text(0.02, 0.5, '', va='center', rotation='vertical',fontsize=16) scenarios = [MCInfectedPopulation( number=1, @@ -853,17 +853,17 @@ def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid: for i in range(3): axs[i].hist(qr_values[3 * i:3 * (i + 1)], bins=2000, histtype='step', - color=colors, cumulative=True, range=(left, right)) - axs[i].set_xlim(left, right) + color=colors, cumulative=True, range=(-7, 6)) + axs[i].set_xlim(-6, 6) axs[i].set_yticks([0, samples / 2, samples]) axs[i].set_yticklabels(['0.0', '0.5', '1.0']) axs[i].yaxis.set_label_position("right") - axs[i].set_ylabel(activities[i],fontsize=12) + axs[i].set_ylabel(activities[i], fontsize=14) axs[i].grid(linestyle='--') axs[0].legend(handles=lines, loc='upper left') - plt.xlabel('qR (q h$^{-1}$)', fontsize=12) - tick_positions = np.arange(int(np.ceil(left)), int(np.ceil(right)), 2) + plt.xlabel('$\mathrm{qR}$ (q h$^{-1}$)', fontsize=16) + tick_positions = np.arange(-6, 6, 2) plt.xticks(ticks=tick_positions, labels=['$\;10^{' + str(i) + '}$' for i in tick_positions]) fig.set_figheight(8) @@ -1258,6 +1258,6 @@ def plot_pi_vs_exposure_time(exp_models: typing.List[MCExposureModel], labels: t plt.title('') plt.xlabel(f'Travel time ({"min" if time_in_minutes else "h"})', fontsize=12) - plt.ylabel('Probability of infection\n$P(I|qID=60)$', fontsize=12) + plt.ylabel('Probability of infection\n$P(\,\mathtt{I}\,|\,\mathrm{vl}\,)$', fontsize=12) plt.legend() plt.show() \ No newline at end of file From ac11392bd330cf4569982c98d40851332934be2e Mon Sep 17 00:00:00 2001 From: Andrejh Date: Mon, 15 Mar 2021 23:01:54 +0100 Subject: [PATCH 2/7] update scenarios half-life --- cara/model_scenarios.py | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/cara/model_scenarios.py b/cara/model_scenarios.py index 105e53e0..5e50d49d 100644 --- a/cara/model_scenarios.py +++ b/cara/model_scenarios.py @@ -129,7 +129,7 @@ classroom_model = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=False, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -156,7 +156,7 @@ classroom_model_full_open = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=False, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -183,7 +183,7 @@ classroom_model_summer = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=False, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -210,7 +210,7 @@ classroom_model_full_open_summer = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=False, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -227,17 +227,15 @@ classroom_model_full_open_summer = [MCExposureModel( classroom_model_no_vent = [MCExposureModel( concentration_model=MCConcentrationModel( room=models.Room(volume=160), - ventilation=models.SlidingWindow( + ventilation=models.AirChange( active=models.PeriodicInterval(period=120, duration=120), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), - outside_temp=models.PiecewiseConstant((0, 24), (291,)), - window_height=1.6, opening_length=0., + air_exch=0.25, ), infected=MCInfectedPopulation( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=False, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -264,7 +262,7 @@ classroom_model_teacher_mask_full_open = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=True, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -297,7 +295,7 @@ classroom_model_with_hepa = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=False, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -330,7 +328,7 @@ classroom_model_full_open_multi = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=False, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -363,7 +361,7 @@ classroom_model_full_open_multi_masks = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), masked=True, - virus=MCVirus(halflife=1.1, qID=qid), + virus=MCVirus(halflife=3.8, qID=qid), expiratory_activity=2, samples=200000, breathing_category=3, @@ -439,7 +437,7 @@ shared_office_worst_model = [MCExposureModel( room=models.Room(volume=50), ventilation=models.AirChange( active=models.PeriodicInterval(period=120, duration=120), - air_exch=0., + air_exch=0.25, ), infected=MCInfectedPopulation( number=1, @@ -788,7 +786,7 @@ waiting_room_model = [MCExposureModel( room=models.Room(volume=100), ventilation=models.AirChange( active=models.PeriodicInterval(period=120, duration=120), - air_exch=0., + air_exch=0.25, ), infected=MCInfectedPopulation( number=1, From e5cbb781ee98d97af739593c484340750d43bc6d Mon Sep 17 00:00:00 2001 From: Andrejh Date: Mon, 15 Mar 2021 23:02:15 +0100 Subject: [PATCH 3/7] min air_exchange 0.25 ACH --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 64962b4c..61ae3d04 100644 --- a/cara/models.py +++ b/cara/models.py @@ -236,7 +236,7 @@ class WindowOpening(Ventilation): def air_exchange(self, room: Room, time: float) -> float: # If the window is shut, no air is being exchanged. if not self.active.triggered(time): - return 0. + return 0.25 # Reminder, no dependence on time in the resulting calculation. inside_temp = self.inside_temp.value(time) From b58840a49bbcab3087f82bc163354103723bd67a Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Thu, 18 Mar 2021 16:22:01 +0000 Subject: [PATCH 4/7] Update mc-output.py --- cara/mc-output.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/cara/mc-output.py b/cara/mc-output.py index e74d9940..77a14046 100644 --- a/cara/mc-output.py +++ b/cara/mc-output.py @@ -1,3 +1,10 @@ +""" Title: +Author: <author(s) names> +Date: <date> +Code version: <code version> +Availability: <where it's located> """ + + from cara.montecarlo import * from cara.model_scenarios import * From 381e03f6d452460e17ca1ea5806cb275ff7ce891 Mon Sep 17 00:00:00 2001 From: Andrejh <andre.henriques@cern.ch> Date: Fri, 19 Mar 2021 10:37:28 +0100 Subject: [PATCH 5/7] add new .py files for publication --- cara/mc-output-publication.py | 3 +++ cara/model_scenarios_publication.py | 3 +++ 2 files changed, 6 insertions(+) create mode 100644 cara/mc-output-publication.py create mode 100644 cara/model_scenarios_publication.py diff --git a/cara/mc-output-publication.py b/cara/mc-output-publication.py new file mode 100644 index 00000000..a4805144 --- /dev/null +++ b/cara/mc-output-publication.py @@ -0,0 +1,3 @@ +from cara.montecarlo import * +from cara.model_scenarios import * + diff --git a/cara/model_scenarios_publication.py b/cara/model_scenarios_publication.py new file mode 100644 index 00000000..689da9f3 --- /dev/null +++ b/cara/model_scenarios_publication.py @@ -0,0 +1,3 @@ +from cara import models +from cara.montecarlo import * + From 418305d50aa32828ea31cd7a48c0904cbb61cd0e Mon Sep 17 00:00:00 2001 From: Andrejh <andre.henriques@cern.ch> Date: Fri, 19 Mar 2021 10:45:21 +0100 Subject: [PATCH 6/7] updates --- cara/mc-output-publication.py | 6 ++++++ cara/mc-output.py | 33 ++++++++++++++++----------------- 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/cara/mc-output-publication.py b/cara/mc-output-publication.py index a4805144..a12bc6aa 100644 --- a/cara/mc-output-publication.py +++ b/cara/mc-output-publication.py @@ -1,3 +1,9 @@ +""" Title: COVID Airborne Risk Assessment +Author: <author(s) names> +Date: <date> +Code version: <code version> +Availability: <where it's located> """ + from cara.montecarlo import * from cara.model_scenarios import * diff --git a/cara/mc-output.py b/cara/mc-output.py index 77a14046..2c186587 100644 --- a/cara/mc-output.py +++ b/cara/mc-output.py @@ -1,10 +1,3 @@ -""" Title: <title of program/source code> -Author: <author(s) names> -Date: <date> -Code version: <code version> -Availability: <where it's located> """ - - from cara.montecarlo import * from cara.model_scenarios import * @@ -14,7 +7,7 @@ from cara.model_scenarios import * #print(np.mean(classroom_model[0].infection_probability())) #print(np.mean(classroom_model[1].infection_probability())) #print(np.mean(chorale_model.infection_probability())+np.std(chorale_model.infection_probability())) -#print(np.quantile(chorale_model.infection_probability(),0.8)) +#print(np.quantile(chorale_model[0].infection_probability(),0.8)) #print(np.quantile(chorale_model.infection_probability(),0.90)) #print(np.quantile(chorale_model.infection_probability(),0.1)) @@ -27,20 +20,20 @@ from cara.model_scenarios import * # time_in_minutes=True, # normalize_y_axis=True) -compare_viruses_qr(violins=True) +#compare_viruses_qr(violins=True) -# print_qd_info(large_population_baselines[0]) +#print_qd_info(waiting_room_model_full_winter[1]) #print(np.mean(shared_office_model[1].infection_probability())) #composite_plot_pi_vs_viral_load([shared_office_worst_model[1], shared_office_model[1], shared_office_better_model[1]], # labels=['No mask &\nwindows closed', 'Baseline:Windows\nopen (10min/2h)', 'Windows open\n(10min/2h) +\nHEPA filter'], # colors=['tomato', '#1f77b4', 'limegreen'], -# title='Shared office scenario', +# title='', # vl_points=200) #composite_plot_pi_vs_viral_load([classroom_model_no_vent[1], classroom_model[1], classroom_model_with_hepa[1], classroom_model_full_open_multi[1], classroom_model_full_open_multi_masks[1]], # labels=['Windows closed', 'Baseline:Windows\nopen (10min/2h)', 'Windows open\n(10min/2h) + HEPA', 'Multiple windows open', 'Multiple windows open\n+masks'], # colors=['tomato','#1f77b4', 'dodgerblue', 'seagreen', 'limegreen'], -# title='Classroom scenario', +# title='', # vl_points=200) #composite_plot_pi_vs_viral_load([ski_cabin_model_60[1], ski_cabin_model_30[1], ski_cabin_model_baseline_20[1], ski_cabin_model_10[1]], @@ -56,16 +49,16 @@ compare_viruses_qr(violins=True) # vl_points=200) #compare_concentration_curves([classroom_model_no_vent[1], classroom_model[1], classroom_model_with_hepa[1], classroom_model_full_open_multi[1]], -# labels=['Windows closed', 'Baseline:(windows 10min/2h)', 'Baseline:(windows 10min/2h) + HEPA', 'Multiple windows open'], +# labels=['Windows closed', 'Baseline: Windows open (10min/2h)', 'Windows open (10min/2h) + HEPA', 'Multiple windows open'], # colors=['tomato','#1f77b4', 'seagreen', 'limegreen'], -# title='Classroom scenario' +# title='' # ) - +# #compare_concentration_curves([waiting_room_model[1], waiting_room_model_full_summer[1], # waiting_room_model_full_winter[1], waiting_room_model_periodic_winter[1]], -# labels=['Baseline:(windows closed)', 'Windows open (summer)', 'Windows open (winter)', 'Windows open 5min/20min (winter)'], +# labels=['Baseline: Windows closed', 'Windows open (summer)', 'Windows open (winter)', 'Windows open 5min/20min (winter)'], # colors=['#1f77b4', 'darkorange', 'deepskyblue', 'lightskyblue'], -# title='Waiting room scenario' +# title='' # ) #plot_pi_vs_viral_load([shared_office_model[1]], labels=['Baseline'],title='') @@ -98,3 +91,9 @@ compare_viruses_qr(violins=True) # plt.title(f'Probability of infection in baseline case - {"English" if model.concentration_model.infected.qid == 60 else "Original"} variant') # plt.yticks([], []) # plt.show() + +#composite_plot_pi_vs_viral_load([shared_office_worst_model[1], shared_office_worst_model_infiltration[1]], +# labels=['No mask &\nwindows closed', 'Baseline:Windows\nopen (10min/2h)', 'Windows open\n(10min/2h) +\nHEPA filter'], +# colors=['tomato', '#1f77b4', 'limegreen'], +# title='Shared office scenario', +# vl_points=200) From 9f6e500618be39281a04dd9426ca0a2fa9962bd8 Mon Sep 17 00:00:00 2001 From: Andrejh <andre.henriques@cern.ch> Date: Fri, 19 Mar 2021 11:31:54 +0100 Subject: [PATCH 7/7] create present_qR_values --- cara/model_scenarios_publication.py | 28 ++++++++++++++++++++++++++++ cara/montecarlo.py | 13 ++++++++++++- 2 files changed, 40 insertions(+), 1 deletion(-) diff --git a/cara/model_scenarios_publication.py b/cara/model_scenarios_publication.py index 689da9f3..e78f2ac2 100644 --- a/cara/model_scenarios_publication.py +++ b/cara/model_scenarios_publication.py @@ -1,3 +1,31 @@ from cara import models from cara.montecarlo import * +######### Standard exposure models ########### +exposure_models = [MCExposureModel( + concentration_model=MCConcentrationModel( + room=models.Room(volume=45), + ventilation=models.SlidingWindow( + active=models.PeriodicInterval(period=120, duration=120), + inside_temp=models.PiecewiseConstant((0, 24), (293,)), + outside_temp=models.PiecewiseConstant((0, 24), (283,)), + window_height=1.6, opening_length=0.2, + ), + infected=MCInfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 4), (5, 9))), + masked=False, + virus=MCVirus(halflife=1.1, qID=100), + expiratory_activity=a, + samples=2000000, + breathing_category=3, + expiratory_activity_weights=(0.7, 0.3, 0) + ) + ), + exposed=models.Population( + number=2, + presence=models.SpecificInterval(((0, 4), (5, 9))), + activity=models.Activity.types['Seated'], + mask=models.Mask.types['Type I'] + ) +) for a in (1, 2, 3)] \ No newline at end of file diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 912f3d4c..fca16db8 100644 --- a/cara/montecarlo.py +++ b/cara/montecarlo.py @@ -426,11 +426,22 @@ def print_qr_info(log_qr: np.ndarray) -> None: print(f"MEAN of log_10(qR) = {np.mean(log_qr)}\n" f"SD of log_10(qR) = {np.std(log_qr)}\n" f"MEAN of qR = {np.mean(qr_values)}\n" - f"SD of qR = {np.std(qr_values)}\n") + f"SD of qR = {np.std(qr_values)}\n" + f"Median of qR = {np.quantile(qr_values, 0.5)}\n") for quantile in (0.01, 0.05, 0.25, 0.50, 0.75, 0.95, 0.99): print(f"qR_{quantile} = {np.quantile(qr_values, quantile)}") +def present_qR_values(model: MCConcentrationModel) -> None: + """ + Displays a handful of key parameters and results of qR from a given MCConcentrationModel + :param model: The MCConcentrationModel representing the scenario to be presented + :return: Nothing, parameters are printed + """ + + qRs = (np.log10(model.infected.emission_rate_when_present())) + + print_qr_info(np.asarray(qRs)) def present_model(model: MCConcentrationModel, bins: int = 200, title: str = 'Summary of $qR$ model parameters') -> None: