From 044f64ccf8cae85cb0f861a00357650d07eec885 Mon Sep 17 00:00:00 2001 From: Andrejh Date: Sat, 27 Feb 2021 23:37:43 +0100 Subject: [PATCH 1/3] dotted axis and increase fontsize --- cara/montecarlo.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 65642c6d..3d02a7ac 100644 --- a/cara/montecarlo.py +++ b/cara/montecarlo.py @@ -1064,6 +1064,7 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label ax.legend(loc='upper left') ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2) + ax.spines["right"].set_visible(False) factors = [0.6 * model.exposed.activity.inhalation_rate * (1 - model.exposed.mask.η_inhale) for model in exp_models] present_indexes = np.array([exp_models[0].exposed.person_present(t) for t in times]) @@ -1075,23 +1076,27 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label for c, factor in zip(modified_concentrations, factors)] plt.suptitle(title) - plt.xlabel("Time ($h$)", fontsize=12) - plt.ylabel("Quantum concentration ($q\;m^{-3}$)\nmean values of $C(t)$", fontsize=12) + plt.xlabel("Time ($h$)", fontsize=14) + plt.ylabel("Quantum concentration ($q\;m^{-3}$)\nmean values", fontsize=14) if show_qd: ax1 = ax.twinx() for qd, label, color in zip(qds, labels, colors): ax1.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted') - ax1.set_ylabel('Dose ($q$)', fontsize=12) + ax1.spines["right"].set_linestyle("--") + ax1.spines["right"].set_linestyle((0,(1,5))) + ax1.set_ylabel('Dose ($q$)', fontsize=14) ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.2) ax2 = ax.twinx() - ax2.spines["right"].set_position(("axes", 1.2)) - ax2.set_ylabel('Dose (RNA copies)\nmean values $qD$', fontsize=12) + ax2.spines["right"].set_position(("axes", 1.15)) + ax2.spines["right"].set_linestyle((0,(1,5))) + ax2.set_ylabel('Dose (RNA copies)\nmean values', fontsize=14) ax2.set_ylim(tuple(y * exp_models[0].concentration_model.virus.qID for y in ax1.get_ylim())) #ax1.legend(loc='upper right') plt.tight_layout() plt.hlines([60], colors=['lightgrey'], linestyles=['dashed'], xmin=start, xmax=stop) - plt.text(7, 65, "$qID = 60$", color='lightgrey') - plt.show() \ No newline at end of file + plt.text(7, 65, "$qID = 60$", color='grey') + plt.show() + print() \ No newline at end of file From dadabcbec4cf4f2dc1885621e4e3af2183b32773 Mon Sep 17 00:00:00 2001 From: Andrejh Date: Mon, 1 Mar 2021 15:25:35 +0100 Subject: [PATCH 2/3] lable changes --- cara/montecarlo.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 3d02a7ac..9f755e60 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 @@ -687,9 +687,9 @@ def composite_plot_pi_vs_viral_load(baselines: typing.List[MCExposureModel], lab axs[1, 0].set_xticks([i for i in range(2, 13, 2)]) 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}$)\n$vl$') - axs[0, 0].set_ylabel('Probability of infection (%)\n$P(I|qID=60)$') - plt.suptitle(title) + 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) + plt.suptitle(title, fontsize=12) axs[0, 0].text(11, -0.01, '$(i)$') axs[1, 0].text(11, axs[1, 0].get_ylim()[1] * 0.8, '$(ii)$') @@ -704,12 +704,12 @@ def composite_plot_pi_vs_viral_load(baselines: typing.List[MCExposureModel], lab break for i, (crit, color) in enumerate(zip(crits, colors)): - axs[0, 0].text(2.5, 0.4 - i * 0.1, f'x $vl_{"{0.95}"}=' + '10^{' + str(np.round(crits[i], 1)) + '}$', fontsize=10, color=color) + axs[0, 0].text(2.5, 0.45 - i * 0.1, f'x $vl_{"{0.95}"}=' + '10^{' + str(np.round(crits[i], 1)) + '}$', fontsize=10, color=color) axs[0, 0].plot(crits[i], 0.95, 'x', color=color) if show_lines: axs[0, 0].hlines([0.5], colors=['lightgrey'], linestyles=['dashed'], xmin=2, xmax=12) - axs[0, 0].text(9.7, 0.52, "$P(I) = 0.5$", color='lightgrey') + axs[0, 0].text(9.7, 0.52, "$P(I) = 0.5$", color='grey') middle_positions = [] for line in lines: for i, point in enumerate(line): From 85564c282995590a78283bda20cdddce3f019d3f Mon Sep 17 00:00:00 2001 From: Andrejh Date: Mon, 1 Mar 2021 15:25:52 +0100 Subject: [PATCH 3/3] add ski expsoure models --- cara/mc-output.py | 32 +++++---- cara/model_scenarios.py | 151 ++++++++++++++++++++++++++++++++++++++-- 2 files changed, 163 insertions(+), 20 deletions(-) diff --git a/cara/mc-output.py b/cara/mc-output.py index e31b7d19..2689bc72 100644 --- a/cara/mc-output.py +++ b/cara/mc-output.py @@ -17,26 +17,32 @@ from cara.model_scenarios import * #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', 'Baseline +\nHEPA filter'], # colors=['tomato', '#1f77b4', 'limegreen'], -# title='$P(I|qID)$ vs $vl$ - Shared office scenario', +# title='Shared office scenario', # 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_hepa[1]], -# labels=['Windows closed', 'Baseline:\n(windows 10min/2h)', 'Baseline:\n(windows 10min/2h)+\nHEPA', 'Multiple windows open +\nHEPA filter'], -# colors=['tomato','#1f77b4', 'seagreen', 'limegreen'], -# title='', +#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 10min/2h)', 'Baseline:(windows 10min/2h)\n+ HEPA', 'Multiple windows open', 'Multiple windows open\n+masks'], +# colors=['tomato','#1f77b4', 'dodgerblue', 'seagreen', 'limegreen'], +# title='Classroom scenario', # vl_points=200) -#compare_concentration_curves([classroom_model_no_vent[1], classroom_model[1], classroom_model_with_hepa[1], classroom_model_full_open_hepa[1]], -# labels=['Windows closed', 'Baseline:(windows 10min/2h)', 'Baseline:(windows 10min/2h) + HEPA', 'Multiple windows open + HEPA filter'], +composite_plot_pi_vs_viral_load([ski_cabin_model_60[1], ski_cabin_model_30[1], ski_cabin_model_20[1], ski_cabin_model_10[1]], + labels=['60 min', '30 min', 'Baseline: 20 min', '10 min'], + colors=['tomato', 'lightsalmon', '#1f77b4', 'limegreen'], + title='Ski cabin scenario', + 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'], # colors=['tomato','#1f77b4', 'seagreen', 'limegreen'], # title='Classroom scenario' # ) -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)'], - colors=['#1f77b4', 'darkorange', 'deepskyblue', 'lightskyblue'], - title='Waiting room scenario' - ) +#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)'], +# colors=['#1f77b4', 'darkorange', 'deepskyblue', 'lightskyblue'], +# title='Waiting room scenario' +# ) #plot_pi_vs_viral_load([shared_office_model[1]], labels=['Baseline'],title='') diff --git a/cara/model_scenarios.py b/cara/model_scenarios.py index 0c39393d..fcc624da 100644 --- a/cara/model_scenarios.py +++ b/cara/model_scenarios.py @@ -311,7 +311,7 @@ classroom_model_with_hepa = [MCExposureModel( ) )for qid in (100, 60)] -classroom_model_full_open_hepa = [MCExposureModel( +classroom_model_full_open_multi = [MCExposureModel( concentration_model=MCConcentrationModel( room=models.Room(volume=160), ventilation=models.MultipleVentilation( @@ -323,7 +323,7 @@ classroom_model_full_open_hepa = [MCExposureModel( window_height=1.6, opening_length=4*0.6, ), models.HEPAFilter(active=models.PeriodicInterval(period=120, duration=120), - q_air_mech=800) + q_air_mech=0) ) ), infected=MCInfectedPopulation( @@ -344,6 +344,39 @@ classroom_model_full_open_hepa = [MCExposureModel( ) )for qid in (100, 60)] +classroom_model_full_open_multi_masks = [MCExposureModel( + concentration_model=MCConcentrationModel( + room=models.Room(volume=160), + ventilation=models.MultipleVentilation( + ventilations=( + 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=4*0.6, + ), + models.HEPAFilter(active=models.PeriodicInterval(period=120, duration=120), + q_air_mech=0) + ) + ), + infected=MCInfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), + masked=True, + virus=MCVirus(halflife=1.1, qID=qid), + expiratory_activity=2, + samples=200000, + breathing_category=3, + ) + ), + exposed=models.Population( + number=19, + presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), + activity=models.Activity.types['Seated'], + mask=models.Mask.types['Type I'] + ) +)for qid in (100, 60)] + ######### Shared office exposure models ########### shared_office_model = [MCExposureModel( concentration_model=MCConcentrationModel( @@ -434,7 +467,7 @@ shared_office_better_model = [MCExposureModel( ) for qid in (100, 60)] ######### Ski cabine exposure models ########### -ski_cabin_model = [MCExposureModel( +ski_cabin_model_20 = [MCExposureModel( concentration_model=MCConcentrationModel( room=models.Room(volume=10), ventilation=models.HVACMechanical( @@ -443,7 +476,7 @@ ski_cabin_model = [MCExposureModel( ), infected=MCInfectedPopulation( number=1, - presence=models.SpecificInterval(((0, 0.33),)), + presence=models.SpecificInterval(((0, 20/60),)), masked=True, virus=MCVirus(halflife=1.1, qID=qid), expiratory_activity=2, @@ -454,7 +487,111 @@ ski_cabin_model = [MCExposureModel( ), exposed=models.Population( number=3, - presence=models.SpecificInterval(((0, 0.33),)), + presence=models.SpecificInterval(((0, 20/60),)), + activity=models.Activity.types['Moderate activity'], + mask=models.Mask.types['Type I'] + ) +)for qid in (100, 60)] + +ski_cabin_model_10 = [MCExposureModel( + concentration_model=MCConcentrationModel( + room=models.Room(volume=10), + ventilation=models.HVACMechanical( + active=models.PeriodicInterval(period=120, duration=120), + q_air_mech=0., + ), + infected=MCInfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 10/60),)), + masked=True, + virus=MCVirus(halflife=1.1, qID=qid), + expiratory_activity=2, + samples=200000, + breathing_category=4, + expiratory_activity_weights=(0.7, 0.3, 0) + ) + ), + exposed=models.Population( + number=3, + presence=models.SpecificInterval(((0, 10/60),)), + activity=models.Activity.types['Moderate activity'], + mask=models.Mask.types['Type I'] + ) +)for qid in (100, 60)] + +ski_cabin_model_25 = [MCExposureModel( + concentration_model=MCConcentrationModel( + room=models.Room(volume=10), + ventilation=models.HVACMechanical( + active=models.PeriodicInterval(period=120, duration=120), + q_air_mech=0., + ), + infected=MCInfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 25/60),)), + masked=True, + virus=MCVirus(halflife=1.1, qID=qid), + expiratory_activity=2, + samples=200000, + breathing_category=4, + expiratory_activity_weights=(0.7, 0.3, 0) + ) + ), + exposed=models.Population( + number=3, + presence=models.SpecificInterval(((0, 25/60),)), + activity=models.Activity.types['Moderate activity'], + mask=models.Mask.types['Type I'] + ) +)for qid in (100, 60)] + +ski_cabin_model_30 = [MCExposureModel( + concentration_model=MCConcentrationModel( + room=models.Room(volume=10), + ventilation=models.HVACMechanical( + active=models.PeriodicInterval(period=120, duration=120), + q_air_mech=0., + ), + infected=MCInfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 30/60),)), + masked=True, + virus=MCVirus(halflife=1.1, qID=qid), + expiratory_activity=2, + samples=200000, + breathing_category=4, + expiratory_activity_weights=(0.7, 0.3, 0) + ) + ), + exposed=models.Population( + number=3, + presence=models.SpecificInterval(((0, 30/60),)), + activity=models.Activity.types['Moderate activity'], + mask=models.Mask.types['Type I'] + ) +)for qid in (100, 60)] + +ski_cabin_model_60 = [MCExposureModel( + concentration_model=MCConcentrationModel( + room=models.Room(volume=10), + ventilation=models.HVACMechanical( + active=models.PeriodicInterval(period=120, duration=120), + q_air_mech=0., + ), + infected=MCInfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 1),)), + masked=True, + virus=MCVirus(halflife=1.1, qID=qid), + expiratory_activity=2, + samples=200000, + breathing_category=4, + expiratory_activity_weights=(0.7, 0.3, 0) + ) + ), + exposed=models.Population( + number=3, + presence=models.SpecificInterval(((0, 1),)), activity=models.Activity.types['Moderate activity'], mask=models.Mask.types['Type I'] ) @@ -555,7 +692,7 @@ waiting_room_model_full_summer = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2),)), masked=False, - virus=MCVirus(halflife=3.8, qID=qid), + virus=MCVirus(halflife=1.1, qID=qid), expiratory_activity=4, samples=200000, breathing_category=1, @@ -611,7 +748,7 @@ waiting_room_model_periodic_summer = [MCExposureModel( number=1, presence=models.SpecificInterval(((0, 2),)), masked=False, - virus=MCVirus(halflife=3.8, qID=qid), + virus=MCVirus(halflife=1.1, qID=qid), expiratory_activity=4, samples=200000, breathing_category=1,