From 0c992e459f56eac5ab0b63b004d4a77ed86e2fa3 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 16 Sep 2021 12:12:52 +0200 Subject: [PATCH] added corrected models for the concentrations plot --- cara/apps/calculator/model_generator.py | 1 - cara/model_scenarios_paper.py | 129 ++++++++++++++---------- cara/models.py | 19 ++++ cara/results_paper.py | 27 ++--- 4 files changed, 105 insertions(+), 71 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 39edb387..d38c983d 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -247,7 +247,6 @@ class FormData: infected=self.infected_population(), ), exposed=self.exposed_population(), - fraction_deposited=self. ) def build_model(self, sample_size=_DEFAULT_MC_SAMPLE_SIZE) -> models.ExposureModel: diff --git a/cara/model_scenarios_paper.py b/cara/model_scenarios_paper.py index 2ff79e9e..0e7453f6 100644 --- a/cara/model_scenarios_paper.py +++ b/cara/model_scenarios_paper.py @@ -489,97 +489,118 @@ def shouting_heavy_exercise_exposure(): return exposure_mc ########## Concentration curves ########### -def classroom_no_mask_windows_closed_exposure(): - exposure_mc = mc.ExposureModel( +def office_model_no_mask_windows_closed(): + office_model_no_vent = mc.ExposureModel( concentration_model=mc.ConcentrationModel( - room=models.Room(volume=100, humidity=0.5), - ventilation=models.AirChange( - active=models.SpecificInterval(((0, 24),)), - air_exch=0.25, - ), + room=models.Room(volume=160, humidity=0.3), + ventilation=models.MultipleVentilation( + (models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.0), + models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25))), infected=mc.InfectedPopulation( number=1, + presence=models.SpecificInterval(present_times = ((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), virus=mc.SARSCoV2( viral_load_in_sputum=symptomatic_vl_frequencies, infectious_dose=infectious_dose_distribution, viable_to_RNA=infectious_virus_distribution, ), - presence=models.SpecificInterval(((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), mask=models.Mask.types["No mask"], - activity=activity_distributions['Light activity'], - expiration=models.Expiration.types['Talking'], - ), + activity=activity_distributions['Seated'], + expiration=models.MultipleExpiration( + expirations = (models.Expiration.types['Talking'], + models.Expiration.types['Breathing']), + weights=(1, 2) + ) + ) ), - exposed=mc.Population( - number=14, - presence=models.SpecificInterval(((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), + exposed=models.Population( + number=18, + presence=models.SpecificInterval(present_times = ((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), activity=models.Activity.types['Seated'], - mask=models.Mask.types["No mask"], - ), + mask=models.Mask.types['No mask'] + ) ) - return exposure_mc + return office_model_no_vent -def classrom_no_mask_windows_open_breaks(): - exposure_mc = mc.ExposureModel( +def office_model_no_mask_windows_open_breaks(): + office_model_no_vent = mc.ExposureModel( concentration_model=mc.ConcentrationModel( - room=models.Room(volume=100, humidity=0.5), - ventilation=models.SlidingWindow( - active=models.SpecificInterval(((1.5, 2), (3.5, 4.5), (6, 6.5))), - inside_temp=models.PiecewiseConstant((0, 24), (295,)), - outside_temp=models.PiecewiseConstant((0, 24), (291,)), - window_height=1.6, - opening_length=0.6, + room=models.Room(volume=160, humidity=0.3), + ventilation = models.MultipleVentilation( + ventilations=( + models.SlidingWindow( + active=models.SpecificInterval(present_times=((1.5, 2), (3.5, 4.5), (6, 6.5))), + inside_temp=models.PiecewiseConstant((0, 24), (295,)), + outside_temp=models.PiecewiseConstant((0, 24), (291,)), + window_height=1.6, + opening_length=0.6, + ), + models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25), + ) ), infected=mc.InfectedPopulation( number=1, + presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), virus=mc.SARSCoV2( viral_load_in_sputum=symptomatic_vl_frequencies, infectious_dose=infectious_dose_distribution, viable_to_RNA=infectious_virus_distribution, ), - presence=models.SpecificInterval(((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), mask=models.Mask.types["No mask"], - activity=activity_distributions['Light activity'], - expiration=models.Expiration.types['Talking'], - ), + activity=activity_distributions['Seated'], + expiration=models.MultipleExpiration( + expirations = (models.Expiration.types['Talking'], + models.Expiration.types['Breathing']), + weights=(1, 2) + ) + ) ), - exposed=mc.Population( - number=14, - presence=models.SpecificInterval(((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), + exposed=models.Population( + number=18, + presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), activity=models.Activity.types['Seated'], - mask=models.Mask.types["No mask"], - ), + mask=models.Mask.types['No mask'] + ) ) - return exposure_mc + return office_model_no_vent -def classrom_no_mask_windows_open_alltimes(): - exposure_mc = mc.ExposureModel( +def office_model_no_mask_windows_open_alltimes(): + office_model_no_vent = mc.ExposureModel( concentration_model=mc.ConcentrationModel( - room=models.Room(volume=100, humidity=0.5), - ventilation=models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=120), - inside_temp=models.PiecewiseConstant((0, 24), (295,)), - outside_temp=models.PiecewiseConstant((0, 24), (291,)), - window_height=1.6, opening_length=0.6, + room=models.Room(volume=160, humidity=0.3), + ventilation=models.MultipleVentilation( + ventilations=( + models.SlidingWindow( + active=models.PeriodicInterval(period=120, duration=120), + inside_temp=models.PiecewiseConstant((0, 24), (295,)), + outside_temp=models.PiecewiseConstant((0, 24), (291,)), + window_height=1.6, opening_length=0.6, + ), + models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25), + ) ), infected=mc.InfectedPopulation( number=1, + presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), virus=mc.SARSCoV2( viral_load_in_sputum=symptomatic_vl_frequencies, infectious_dose=infectious_dose_distribution, viable_to_RNA=infectious_virus_distribution, ), - presence=models.SpecificInterval(((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), mask=models.Mask.types["No mask"], - activity=activity_distributions['Light activity'], - expiration=models.Expiration.types['Talking'], - ), + activity=activity_distributions['Seated'], + expiration=models.MultipleExpiration( + expirations = (models.Expiration.types['Talking'], + models.Expiration.types['Breathing']), + weights=(1, 2) + ) + ) ), - exposed=mc.Population( - number=14, - presence=models.SpecificInterval(((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), + exposed=models.Population( + number=18, + presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))), activity=models.Activity.types['Seated'], - mask=models.Mask.types["No mask"], - ), + mask=models.Mask.types['No mask'] + ) ) - return exposure_mc \ No newline at end of file + return office_model_no_vent \ No newline at end of file diff --git a/cara/models.py b/cara/models.py index 8ebde1ae..f3b72c27 100644 --- a/cara/models.py +++ b/cara/models.py @@ -924,6 +924,25 @@ class ExposureModel: DF = IF * (0.0587 + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d))))) fraction_deposited: _VectorisedFloat = DF + def _normed_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: + """The number of virions per meter^3 between any two times, normalized + by the emission rate of the infected population""" + for start, stop in self.exposed.presence.boundaries(): + if start > time2: + normed_exposure = 0. + break + elif time2 <= stop: + normed_exposure = self.concentration_model.normed_integrated_concentration(time1, time2) + break + else: + normed_exposure = self.concentration_model.normed_integrated_concentration(time1, time2) + return normed_exposure + + def exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: + """The number of virions per meter^3 between any two times.""" + return (self._normed_exposure_between_bounds(time1, time2) * + self.concentration_model.infected.emission_rate_when_present()) + def _normed_exposure(self) -> _VectorisedFloat: """ The number of virus per meter^3, normalized by the emission rate diff --git a/cara/results_paper.py b/cara/results_paper.py index e38e3c9e..e8476c44 100644 --- a/cara/results_paper.py +++ b/cara/results_paper.py @@ -15,7 +15,7 @@ import math ######### Plot material ######### -SAMPLE_SIZE = 50000 +SAMPLE_SIZE = 250000 viral_loads = np.linspace(2, 12, 600) ############# Markers (for legend) ############# @@ -528,9 +528,9 @@ def calculate_deposition_factor(): ############ Compare concentration curves ############ def compare_concentration_curves(): - exp_models=[classroom_no_mask_windows_closed_exposure().build_model(size=SAMPLE_SIZE), - classrom_no_mask_windows_open_breaks().build_model(size=SAMPLE_SIZE), - classrom_no_mask_windows_open_alltimes().build_model(size=SAMPLE_SIZE)] + exp_models=[office_model_no_mask_windows_closed().build_model(size=SAMPLE_SIZE), + office_model_no_mask_windows_open_breaks().build_model(size=SAMPLE_SIZE), + office_model_no_mask_windows_open_alltimes().build_model(size=SAMPLE_SIZE)] labels=['Windows closed', 'Window open during breaks', 'Window open at all times'] colors=['tomato', 'lightskyblue', 'limegreen', '#1f77b4', 'seagreen', 'lightskyblue', 'deepskyblue'] @@ -550,25 +550,20 @@ def compare_concentration_curves(): ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2) ax.spines["right"].set_visible(False) - # When merged, use the integrated_concentration function. - 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]) - modified_concentrations = [np.array(c) for c in concentrations] - for mc in modified_concentrations: - mc[~present_indexes] = 0 - - cumulative_doses = [[np.trapz(c[:i + 1], times[:i + 1]) * factor for i in range(len(times))] - for c, factor in zip(modified_concentrations, factors)] + cumulative_doses = [np.cumsum([ + np.array(model.exposure_between_bounds(float(time1), float(time2))).mean() + for time1, time2 in zip(times[:-1], times[1:]) + ]) for model in exp_models] plt.xlabel("Exposure time ($h$)", fontsize=14) - plt.ylabel("Mean viral concentration\n(virion m$^{-3}$)", fontsize=14) + plt.ylabel("Mean concentration (virions m$^{-3}$)", fontsize=14) ax1 = ax.twinx() for qd, label, color in zip(cumulative_doses, labels, colors): - ax1.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted') + ax1.plot(times[:-1], qd, label='qD - ' + label, color=color, linestyle='dotted') ax1.spines["right"].set_linestyle("--") ax1.spines["right"].set_linestyle((0,(1,5))) - ax1.set_ylabel('Mean cumulative dose\n(virion)', fontsize=14) + ax1.set_ylabel('Mean cumulative dose (virions)', fontsize=14) ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.2) plt.show()