From f99dcc5e0cefa7261ec454b9c0a84aa160c28639 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 16:50:17 +0200 Subject: [PATCH 01/10] Fix wrong emission rate in test_monte_carlo_full_models.py::test_small_shared_office_Geneva; decrease tolerance and number of samples subsequently --- cara/tests/test_monte_carlo_full_models.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index b97956b9..caed9c52 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -9,8 +9,8 @@ from cara.apps.calculator.model_generator import build_expiration # TODO: seed better the random number generators np.random.seed(2000) -SAMPLE_SIZE = 600_000 -TOLERANCE = 0.06 +SAMPLE_SIZE = 500_000 +TOLERANCE = 0.05 # Load the weather data (temperature in kelvin) for Toronto. toronto_coordinates = (43.667, 79.400) @@ -342,7 +342,7 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, ["No mask", "Jul", 9.52, 9.920, 809], ["Type I", "Jul", 1.7, 0.913, 149], ["FFP2", "Jul", 0.51, 0.239, 149], - ["Type I", "Feb", 0.57, 0.272, 162], + ["Type I", "Feb", 0.57, 0.272, 149], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, From ef84f28ab6d298b52d6dc8f0940f89f716ac8031 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 17:13:15 +0200 Subject: [PATCH 02/10] Adding a method to extract the right interaction time boundaries between two given times in the short-range model (cosmetics) --- cara/models.py | 52 ++++++++++++++++++++++++++++++++------------------ 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/cara/models.py b/cara/models.py index bfa50139..8097fd71 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1168,19 +1168,44 @@ class ShortRangeModel: # calculations for the same time (e.g. at state change times). return self._normed_concentration(concentration_model, time) - def normed_exposure_between_bounds(self, concentration_model: ConcentrationModel, time1: float, time2: float): + @method_cache + def extract_between_bounds(self, time1: float, time2: float) -> typing.Tuple[float,float]: + """ + Extract the bounds of the interval resulting from the + intersection of [time1, time2] and the presence interval. + If [time1, time2] has nothing common to the presence interval, + we return (0, 0). + Raise an error if time1 and time2 are not in ascending order. + """ + if time1>time2: + raise ValueError("time1 must be less or equal to time2") + + start, stop = self.presence.boundaries()[0] + if (stop < time1) or (start > time2): + return (0, 0) + elif start <= time1 and time2<= stop: + return time1, time2 + elif start <= time1 and stop < time2: + return time1, stop + elif time1 < start and time2 <= stop: + return start, time2 + elif time1 <= start and stop < time2: + return start, stop + + def normed_exposure_between_bounds(self, concentration_model: ConcentrationModel, + time1: float, time2: float): """ Get the integrated short-range concentration of viruses in the air between the times start and stop, normalized by the virus viral load. """ - start_bound, stop_bound = self.presence.boundaries()[0] + start, stop = self.extract_between_bounds(time1, time2) jet_origin = self.expiration.jet_origin_concentration() dilution = self.dilution_factor() total_normed_concentration_diluted = ( - concentration_model.integrated_concentration(start_bound, - stop_bound)/dilution/ + concentration_model.integrated_concentration(start, stop + )/dilution/ concentration_model.virus.viral_load_in_sputum ) total_normed_concentration_interpolated = np.interp( @@ -1188,7 +1213,7 @@ class ShortRangeModel: concentration_model.infected.particle.diameter, total_normed_concentration_diluted ) - return (jet_origin/dilution * (stop_bound - start_bound) + return (jet_origin/dilution * (stop - start) ) - total_normed_concentration_interpolated @@ -1297,20 +1322,9 @@ class ExposureModel: """ deposited_exposure = 0. for interaction in self.short_range: - start, stop = interaction.presence.boundaries()[0] - if stop < time1: - continue - elif start > time2: - break - elif start <= time1 and time2<= stop: - start_bound, stop_bound = time1, time2 - elif start <= time1 and stop < time2: - start_bound, stop_bound = time1, stop - elif time1 < start and time2 <= stop: - start_bound, stop_bound = start, time2 - elif time1 <= start and stop < time2: - start_bound, stop_bound = start, stop - short_range_exposure = interaction.normed_exposure_between_bounds(self.concentration_model, start_bound, stop_bound) + start, stop = interaction.extract_between_bounds(time1, time2) + short_range_exposure = interaction.normed_exposure_between_bounds( + self.concentration_model, start, stop) fdep = interaction.expiration.particle.fraction_deposited(evaporation_factor=1.0) diameter = interaction.expiration.particle.diameter From b11dd90e928f316701a41113b3555c353b25a5c4 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 22:21:15 +0200 Subject: [PATCH 03/10] Adding full algo tests with distributions, in particular on the concentration/dose/probability of infection --- cara/tests/test_full_algorithm.py | 247 ++++++++++++++++++++++++------ 1 file changed, 197 insertions(+), 50 deletions(-) diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index fb6ade81..6726f3a4 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -244,7 +244,7 @@ class SimpleShortRangeModel: return dilution - @method_cache +# @method_cache def jet_concentration(self,conc_model: SimpleConcentrationModel) -> _VectorisedFloat: """ virion concentration at the origin of the jet (close to @@ -267,13 +267,13 @@ class SimpleShortRangeModel: def concentration(self, conc_model: SimpleConcentrationModel, time: float) -> _VectorisedFloat: """ compute the short-range part of the concentration, and add it - to the background concentration + to the long-range concentration """ if self.interaction_interval.triggered(time): - background_concentration = conc_model.concentration(time) + lr_concentration = conc_model.concentration(time) S = self.dilution_factor() return (self.jet_concentration(conc_model) - - background_concentration) / S + - lr_concentration) / S else: return 0. @@ -352,8 +352,17 @@ class SimpleExposureModel(SimpleConcentrationModel): epsabs=0.,limit=500)[0] * self.viral_load * self.breathing_rate) + def total_concentration(self, t: float): + """ + total concentration at time t + """ + res = self.concentration(t) + for sr_mod in self.sr_models: + res += sr_mod.concentration(self,t) + return res + @method_cache - def integrated_background_concentration(self,t1: float,t2: float, + def integrated_longrange_concentration(self,t1: float,t2: float, evaporation: float) -> _VectorisedFloat: """ background (long-range) concentration integrated from t1 to t2 @@ -417,7 +426,7 @@ class SimpleExposureModel(SimpleConcentrationModel): epsabs=0.,limit=500)[0] * self.viral_load * 1e-6 * (t2-t1) ) result += sr_model.breathing_rate * ( - res-self.integrated_background_concentration(t1,t2,evaporation) + res-self.integrated_longrange_concentration(t1,t2,evaporation) )/sr_model.dilution_factor() return result @@ -429,7 +438,7 @@ class SimpleExposureModel(SimpleConcentrationModel): """ result = 0. for t1,t2 in self.infected_presence.boundaries(): - result += (self.integrated_background_concentration(t1,t2,self.evaporation) + result += (self.integrated_longrange_concentration(t1,t2,self.evaporation) * self.breathing_rate) result += self.integrated_shortrange_concentration() @@ -468,6 +477,25 @@ def c_model() -> mc.ConcentrationModel: ).build_model(SAMPLE_SIZE) +@pytest.fixture +def c_model_distr() -> mc.ConcentrationModel: + return mc.ConcentrationModel( + room=models.Room(volume=50, humidity=0.3), + ventilation=models.AirChange(active=models.PeriodicInterval( + period=120, duration=120), air_exch=1.), + infected=mc.InfectedPopulation( + number=1, + presence=presence, + virus=virus_distributions['SARS_CoV_2_DELTA'], + mask=models.Mask.types['No mask'], + activity=activity_distributions['Seated'], + expiration=expiration_distributions['Breathing'], + host_immunity=0., + ), + evaporation_factor=0.3, + ).build_model(SAMPLE_SIZE) + + @pytest.fixture def sr_models() -> typing.Tuple[mc.ShortRangeModel, ...]: return ( @@ -516,10 +544,107 @@ def simple_sr_models() -> typing.Tuple[SimpleShortRangeModel, ...]: ) +@pytest.fixture +def expo_sr_model(c_model,sr_models) -> mc.ExposureModel: + return mc.ExposureModel( + concentration_model=c_model, + short_range=sr_models, + exposed=mc.Population( + number=1, + presence=presence, + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + host_immunity=0., + ), + ).build_model(SAMPLE_SIZE) + + +@pytest.fixture +def simple_expo_sr_model(simple_sr_models) -> SimpleExposureModel: + return SimpleExposureModel( + infected_presence = presence, + viral_load = models.Virus.types['SARS_CoV_2_DELTA'].viral_load_in_sputum, + breathing_rate = models.Activity.types['Seated'].exhalation_rate, + room_volume = 50., + lambda_ventilation= 1., + BLO_factors = expiration_BLO_factors['Breathing'], + finf = models.Virus.types['SARS_CoV_2_DELTA'].viable_to_RNA_ratio, + HI = 0., + ID50 = models.Virus.types['SARS_CoV_2_DELTA'].infectious_dose, + transmissibility = models.Virus.types['SARS_CoV_2_DELTA'].transmissibility_factor, + sr_models = simple_sr_models, + ) + + +@pytest.fixture +def expo_sr_model_distr(c_model_distr) -> mc.ExposureModel: + return mc.ExposureModel( + concentration_model=c_model_distr, + short_range=( + mc.ShortRangeModel( + expiration = short_range_expiration_distributions['Breathing'], + activity = activity_distributions['Seated'], + presence = interaction_intervals[0], + distance = short_range_distances, + ).build_model(SAMPLE_SIZE), + mc.ShortRangeModel( + expiration = short_range_expiration_distributions['Speaking'], + activity = activity_distributions['Seated'], + presence = interaction_intervals[1], + distance = short_range_distances, + ).build_model(SAMPLE_SIZE), + ), + exposed=mc.Population( + number=1, + presence=presence, + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + host_immunity=0., + ), + ).build_model(SAMPLE_SIZE) + + +@pytest.fixture +def simple_expo_sr_model_distr(c_model_distr) -> SimpleExposureModel: + return SimpleExposureModel( + infected_presence = presence, + viral_load = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).viral_load_in_sputum, + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + room_volume = 50., + lambda_ventilation= 1., + BLO_factors = expiration_BLO_factors['Breathing'], + finf = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).viable_to_RNA_ratio, + HI = 0., + ID50 = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).infectious_dose, + transmissibility = virus_distributions['SARS_CoV_2_DELTA' + ].transmissibility_factor, + sr_models = ( + SimpleShortRangeModel( + interaction_interval = interaction_intervals[0], + distance = short_range_distances.generate_samples(SAMPLE_SIZE), + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + BLO_factors = expiration_BLO_factors['Breathing'], + ), + SimpleShortRangeModel( + interaction_interval = interaction_intervals[1], + distance = short_range_distances.generate_samples(SAMPLE_SIZE), + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + BLO_factors = expiration_BLO_factors['Speaking'], + ) + ), + ) + + @pytest.mark.parametrize( "time", np.linspace(8.5,17.5,12), ) -def test_background_concentration(time,c_model,simple_c_model): +def test_longrange_concentration(time,c_model,simple_c_model): npt.assert_allclose( c_model.concentration(time).mean(), simple_c_model.concentration(time), rtol=TOLERANCE @@ -542,7 +667,7 @@ def test_shortrange_concentration(time,c_model,simple_c_model, ) -def test_background_exposure(c_model): +def test_longrange_exposure(c_model): simple_expo_model = SimpleExposureModel( infected_presence = presence, viral_load = models.Virus.types['SARS_CoV_2_DELTA'].viral_load_in_sputum, @@ -577,7 +702,27 @@ def test_background_exposure(c_model): ) -def test_background_exposure_with_distributions(): +@pytest.mark.parametrize( + "time", [11., 12.5, 17.] +) +def test_longrange_concentration_with_distributions(c_model_distr,time): + simple_expo_model = SimpleConcentrationModel( + infected_presence = presence, + viral_load = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).viral_load_in_sputum, + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + room_volume = 50., + lambda_ventilation= 1., + BLO_factors = expiration_BLO_factors['Breathing'], + ) + npt.assert_allclose( + c_model_distr.concentration(time).mean(), + simple_expo_model.concentration(time).mean(), rtol=TOLERANCE + ) + + +def test_longrange_exposure_with_distributions(c_model_distr): simple_expo_model = SimpleExposureModel( infected_presence = presence, viral_load = virus_distributions['SARS_CoV_2_DELTA' @@ -597,21 +742,7 @@ def test_background_exposure_with_distributions(): sr_models = (), ) expo_model = mc.ExposureModel( - concentration_model=mc.ConcentrationModel( - room=models.Room(volume=50, humidity=0.3), - ventilation=models.AirChange(active=models.PeriodicInterval( - period=120, duration=120), air_exch=1.), - infected=mc.InfectedPopulation( - number=1, - presence=presence, - virus=virus_distributions['SARS_CoV_2_DELTA'], - mask=models.Mask.types['No mask'], - activity=activity_distributions['Seated'], - expiration=expiration_distributions['Breathing'], - host_immunity=0., - ), - evaporation_factor=0.3, - ), + concentration_model=c_model_distr, short_range=(), exposed=mc.Population( number=1, @@ -631,31 +762,21 @@ def test_background_exposure_with_distributions(): ) -def test_exposure_with_shortrange(c_model,sr_models,simple_sr_models): - simple_expo_sr_model = SimpleExposureModel( - infected_presence = presence, - viral_load = models.Virus.types['SARS_CoV_2_DELTA'].viral_load_in_sputum, - breathing_rate = models.Activity.types['Seated'].exhalation_rate, - room_volume = 50., - lambda_ventilation= 1., - BLO_factors = expiration_BLO_factors['Breathing'], - finf = models.Virus.types['SARS_CoV_2_DELTA'].viable_to_RNA_ratio, - HI = 0., - ID50 = models.Virus.types['SARS_CoV_2_DELTA'].infectious_dose, - transmissibility = models.Virus.types['SARS_CoV_2_DELTA'].transmissibility_factor, - sr_models = simple_sr_models, - ) - expo_sr_model = mc.ExposureModel( - concentration_model=c_model, - short_range=sr_models, - exposed=mc.Population( - number=1, - presence=presence, - mask=models.Mask.types['No mask'], - activity=models.Activity.types['Seated'], - host_immunity=0., - ), - ).build_model(SAMPLE_SIZE) +# tests on the concentration with short-range should be skipped until +# one finds a way to avoid the large variability of the concentration +# with short-range 'Speaking' or 'Shouting' interactions +@pytest.mark.skip +@pytest.mark.parametrize( + "time", [10.75, 14.75, 16.] +) +def test_concentration_with_shortrange(expo_sr_model,simple_expo_sr_model,time): + npt.assert_allclose( + expo_sr_model.concentration(time).mean(), + simple_expo_sr_model.total_concentration(time).mean(), rtol=TOLERANCE + ) + + +def test_exposure_with_shortrange(expo_sr_model,simple_expo_sr_model): npt.assert_allclose( expo_sr_model.deposited_exposure().mean(), simple_expo_sr_model.dose().mean(), rtol=TOLERANCE @@ -665,3 +786,29 @@ def test_exposure_with_shortrange(c_model,sr_models,simple_sr_models): simple_expo_sr_model.probability_infection().mean(), rtol=TOLERANCE ) + +@pytest.mark.skip +@pytest.mark.parametrize( + "time", [10.75, 14.75, 16.] +) +def test_concentration_with_shortrange_and_distributions( + expo_sr_model_distr,simple_expo_sr_model_distr,time): + npt.assert_allclose( + expo_sr_model_distr.concentration(time).mean(), + simple_expo_sr_model_distr.total_concentration(time).mean(), + rtol=TOLERANCE + ) + + +def test_exposure_with_shortrange_and_distributions(expo_sr_model_distr, + simple_expo_sr_model_distr): + npt.assert_allclose( + expo_sr_model_distr.deposited_exposure().mean(), + simple_expo_sr_model_distr.dose().mean(), rtol=0.05 + ) + npt.assert_allclose( + expo_sr_model_distr.infection_probability().mean(), + simple_expo_sr_model_distr.probability_infection().mean(), + rtol=0.03 + ) + From 9a4790d407ead266ead34f6d61f1cecd6d12d9b1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 22:23:49 +0200 Subject: [PATCH 04/10] Decoupling jet exposure and the part subtracted from the background concentration, in the short-range model, in order to treat correctly the diameter integrals. This correct the probability of infection. --- cara/models.py | 70 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 49 insertions(+), 21 deletions(-) diff --git a/cara/models.py b/cara/models.py index 8097fd71..cbb547c9 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1192,29 +1192,46 @@ class ShortRangeModel: elif time1 <= start and stop < time2: return start, stop - def normed_exposure_between_bounds(self, concentration_model: ConcentrationModel, - time1: float, time2: float): + def _normed_jet_exposure_between_bounds(self, + concentration_model: ConcentrationModel, + time1: float, time2: float): """ - Get the integrated short-range concentration of viruses in the air between the times start and stop, - normalized by the virus viral load. + Get the part of the integrated short-range concentration of + viruses in the air, between the times start and stop, coming + from the jet concentration, normalized by the viral load, and + without dilution. """ start, stop = self.extract_between_bounds(time1, time2) - jet_origin = self.expiration.jet_origin_concentration() - dilution = self.dilution_factor() + return jet_origin * (stop - start) - total_normed_concentration_diluted = ( - concentration_model.integrated_concentration(start, stop - )/dilution/ - concentration_model.virus.viral_load_in_sputum + def _normed_interpolated_longrange_exposure_between_bounds( + self, concentration_model: ConcentrationModel, + time1: float, time2: float): + """ + Get the part of the integrated short-range concentration due + to the background concentration, normalized by the viral load + and the breathing rate, and without dilution. + One needs to interpolate the integrated long-range concentration + for the particle diameters defined here. + TODO: make sure any potential extrapolation has a + negligible effect. + """ + start, stop = self.extract_between_bounds(time1, time2) + if stop<=start: + return 0. + + normed_int_concentration = ( + concentration_model.integrated_concentration(start, stop) + /concentration_model.virus.viral_load_in_sputum + /concentration_model.infected.activity.exhalation_rate ) - total_normed_concentration_interpolated = np.interp( + normed_int_concentration_interpolated = np.interp( self.expiration.particle.diameter, concentration_model.infected.particle.diameter, - total_normed_concentration_diluted + normed_int_concentration ) - return (jet_origin/dilution * (stop - start) - ) - total_normed_concentration_interpolated + return normed_int_concentration_interpolated @dataclass(frozen=True) @@ -1292,7 +1309,7 @@ class ExposureModel: # we compute first the mean of all diameter-dependent quantities # to perform properly the Monte-Carlo integration over # particle diameters (doing things in another order would - # lead to wrong results). + # lead to wrong results for the probability of infection). dep_exposure_integrated = np.array(self._long_range_normed_exposure_between_bounds(time1, time2) * aerosols * fdep).mean() @@ -1323,24 +1340,35 @@ class ExposureModel: deposited_exposure = 0. for interaction in self.short_range: start, stop = interaction.extract_between_bounds(time1, time2) - short_range_exposure = interaction.normed_exposure_between_bounds( + short_range_jet_exposure = interaction._normed_jet_exposure_between_bounds( self.concentration_model, start, stop) + short_range_lr_exposure = interaction._normed_interpolated_longrange_exposure_between_bounds( + self.concentration_model, start, stop) + dilution = interaction.dilution_factor() fdep = interaction.expiration.particle.fraction_deposited(evaporation_factor=1.0) diameter = interaction.expiration.particle.diameter - # Aerosols not considered given the formula for the initial concentration at mouth/nose. + # Aerosols not considered given the formula for the initial + # concentration at mouth/nose. if diameter is not None and not np.isscalar(diameter): # we compute first the mean of all diameter-dependent quantities # to perform properly the Monte-Carlo integration over # particle diameters (doing things in another order would - # lead to wrong results). - deposited_exposure += np.array(short_range_exposure * - fdep).mean() + # lead to wrong results for the probability of infection). + deposited_exposure += (np.array(short_range_jet_exposure + * fdep).mean()/dilution + - np.array(short_range_lr_exposure * fdep).mean() + * self.concentration_model.infected.activity.exhalation_rate + /dilution) else: # in the case of a single diameter or no diameter defined, # one should not take any mean at this stage. - deposited_exposure += short_range_exposure*fdep + deposited_exposure += (short_range_jet_exposure + * fdep/dilution + - short_range_lr_exposure * fdep + * self.concentration_model.infected.activity.exhalation_rate + /dilution) # multiply by the (diameter-independent) inhalation rate deposited_exposure *= interaction.activity.inhalation_rate From a86d9a838a766cb670fdd8e3156b0fe34ddaa445 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 22:25:50 +0200 Subject: [PATCH 05/10] In test_short_range_model, removing test on short-range normed exposure (method does not exist anymore), and adding tests on extract_between_bounds method --- cara/tests/models/test_short_range_model.py | 48 +++++++++++++-------- 1 file changed, 29 insertions(+), 19 deletions(-) diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index a07a656b..b1dee050 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -45,7 +45,8 @@ def test_short_range_model_ndarray(concentration_model, short_range_model): model = short_range_model.build_model(250_000) assert isinstance(model._normed_concentration(concentration_model, 10.75), np.ndarray) assert isinstance(model.short_range_concentration(concentration_model, 10.75), np.ndarray) - assert isinstance(model.normed_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) + assert isinstance(model._normed_jet_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) + assert isinstance(model._normed_interpolated_longrange_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) assert isinstance(model.short_range_concentration(concentration_model, 14.0), float) @@ -69,6 +70,32 @@ def test_dilution_factor(activity, expected_dilution): ) +def test_extract_between_bounds_raise_on_wrong_order(short_range_model): + model = short_range_model.build_model(1) + with pytest.raises(ValueError, match='time1 must be less or equal to time2'): + model.extract_between_bounds(11.,10.) + + +@pytest.mark.parametrize( + "time1, time2, expected_start, expected_stop", [ + [10., 12., 10.5, 11.], + [10., 10.7, 10.5, 10.7], + [10., 10.45, 0., 0.], + [11.01, 11.5, 0., 0.], + [10.8, 10.9, 10.8, 10.9], + [10.8, 11.5, 10.8, 11.], + [10.5, 11., 10.5, 11.], + ] +) +def test_extract_between_bounds(short_range_model, time1, time2, + expected_start, expected_stop): + model = short_range_model.build_model(1) + np.testing.assert_equal( + model.extract_between_bounds(time1, time2), + (expected_start, expected_stop), + ) + + @pytest.mark.parametrize( "time, expected_short_range_concentration", [ [8.5, 0.], @@ -83,23 +110,6 @@ def test_short_range_concentration(time, expected_short_range_concentration, con model = short_range_model.build_model(250_000) np.testing.assert_allclose( np.array(model.short_range_concentration(concentration_model, time)).mean(), - expected_short_range_concentration, rtol=0.01 + expected_short_range_concentration, rtol=0.02 ) -@pytest.mark.parametrize( - "start, stop, expected_exposure", [ - [8.5, 12.5, 7.875963317294013e-09], - [10.5, 11.0, 7.875963317294013e-09], - [10.4, 11.1, 7.875963317294013e-09], - [10.5, 11.1, 7.875963317294013e-09], - [10.6, 11.1, 7.66539809488759e-09], - [10.4, 10.9, 7.66539809488759e-09], - - ] -) -def test_normed_exposure_between_bounds(start, stop, expected_exposure, concentration_model, short_range_model): - concentration_model = concentration_model.build_model(250_000) - model = short_range_model.build_model(250_000) - np.testing.assert_almost_equal( - model.normed_exposure_between_bounds(concentration_model, start, stop).mean(), expected_exposure - ) From 6d9210939d2fee299cc78b3d5dcd77db68f0d0ff Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sat, 23 Apr 2022 09:45:03 +0200 Subject: [PATCH 06/10] Updating version number --- cara/apps/calculator/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index e7c1d5d6..8d4e5e45 100644 --- a/cara/apps/calculator/__init__.py +++ b/cara/apps/calculator/__init__.py @@ -33,7 +33,7 @@ from .user import AuthenticatedUser, AnonymousUser # calculator version. If the calculator needs to make breaking changes (e.g. change # form attributes) then it can also increase its MAJOR version without needing to # increase the overall CARA version (found at ``cara.__version__``). -__version__ = "4.1.1" +__version__ = "4.1.2" class BaseRequestHandler(RequestHandler): From ac49e102d4bb3e1ce02673275b8a7a4b11b29c6d Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 25 Apr 2022 15:44:27 +0200 Subject: [PATCH 07/10] Removing commented decorator in test_full_algorithm.py; inverting order Speaking/Breathing in short-range tests in test_full_algorithm.py --- cara/tests/test_full_algorithm.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index 6726f3a4..372ca8b1 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -19,7 +19,7 @@ from cara.monte_carlo.data import (expiration_distributions, # TODO: seed better the random number generators np.random.seed(2000) SAMPLE_SIZE = 1_000_000 -TOLERANCE = 0.02 +TOLERANCE = 0.04 sqrt2pi = np.sqrt(2.*np.pi) sqrt2 = np.sqrt(2.) @@ -244,7 +244,6 @@ class SimpleShortRangeModel: return dilution -# @method_cache def jet_concentration(self,conc_model: SimpleConcentrationModel) -> _VectorisedFloat: """ virion concentration at the origin of the jet (close to @@ -500,14 +499,14 @@ def c_model_distr() -> mc.ConcentrationModel: def sr_models() -> typing.Tuple[mc.ShortRangeModel, ...]: return ( mc.ShortRangeModel( - expiration = short_range_expiration_distributions['Breathing'], + expiration = short_range_expiration_distributions['Speaking'], activity = models.Activity.types['Seated'], presence = interaction_intervals[0], distance = 0.854, ).build_model(SAMPLE_SIZE), mc.ShortRangeModel( - expiration = short_range_expiration_distributions['Speaking'], - activity = models.Activity.types['Seated'], + expiration = short_range_expiration_distributions['Breathing'], + activity = models.Activity.types['Heavy exercise'], presence = interaction_intervals[1], distance = 0.854, ).build_model(SAMPLE_SIZE), @@ -533,14 +532,14 @@ def simple_sr_models() -> typing.Tuple[SimpleShortRangeModel, ...]: interaction_interval = interaction_intervals[0], distance = 0.854, breathing_rate = models.Activity.types['Seated'].exhalation_rate, - BLO_factors = expiration_BLO_factors['Breathing'], + BLO_factors = expiration_BLO_factors['Speaking'], ), SimpleShortRangeModel( interaction_interval = interaction_intervals[1], distance = 0.854, - breathing_rate = models.Activity.types['Seated'].exhalation_rate, - BLO_factors = expiration_BLO_factors['Speaking'], - ) + breathing_rate = models.Activity.types['Heavy exercise'].exhalation_rate, + BLO_factors = expiration_BLO_factors['Breathing'], + ), ) From 2d1f7fa823f79f6a58f71cedd33b81ab814e9868 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 25 Apr 2022 17:25:59 +0200 Subject: [PATCH 08/10] Fixing bug pointed in issue #262 --- cara/models.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cara/models.py b/cara/models.py index cbb547c9..2c65c770 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1356,7 +1356,7 @@ class ExposureModel: # to perform properly the Monte-Carlo integration over # particle diameters (doing things in another order would # lead to wrong results for the probability of infection). - deposited_exposure += (np.array(short_range_jet_exposure + this_deposited_exposure = (np.array(short_range_jet_exposure * fdep).mean()/dilution - np.array(short_range_lr_exposure * fdep).mean() * self.concentration_model.infected.activity.exhalation_rate @@ -1364,14 +1364,15 @@ class ExposureModel: else: # in the case of a single diameter or no diameter defined, # one should not take any mean at this stage. - deposited_exposure += (short_range_jet_exposure + this_deposited_exposure = (short_range_jet_exposure * fdep/dilution - short_range_lr_exposure * fdep * self.concentration_model.infected.activity.exhalation_rate /dilution) # multiply by the (diameter-independent) inhalation rate - deposited_exposure *= interaction.activity.inhalation_rate + deposited_exposure += (this_deposited_exposure * + interaction.activity.inhalation_rate) # then we multiply by diameter-independent quantities: viral load # and fraction of infected virions From faaf1ed086068ee9f845fa289977bea0f9c08e11 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 25 Apr 2022 17:32:30 +0200 Subject: [PATCH 09/10] In models, ExposureModel: multiplication by dilution factor done in a better place, in deposited_exposure_between_bounds method --- cara/models.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/cara/models.py b/cara/models.py index 2c65c770..9d7898fd 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1357,22 +1357,20 @@ class ExposureModel: # particle diameters (doing things in another order would # lead to wrong results for the probability of infection). this_deposited_exposure = (np.array(short_range_jet_exposure - * fdep).mean()/dilution + * fdep).mean() - np.array(short_range_lr_exposure * fdep).mean() - * self.concentration_model.infected.activity.exhalation_rate - /dilution) + * self.concentration_model.infected.activity.exhalation_rate) else: # in the case of a single diameter or no diameter defined, # one should not take any mean at this stage. - this_deposited_exposure = (short_range_jet_exposure - * fdep/dilution + this_deposited_exposure = (short_range_jet_exposure * fdep - short_range_lr_exposure * fdep - * self.concentration_model.infected.activity.exhalation_rate - /dilution) + * self.concentration_model.infected.activity.exhalation_rate) # multiply by the (diameter-independent) inhalation rate deposited_exposure += (this_deposited_exposure * - interaction.activity.inhalation_rate) + interaction.activity.inhalation_rate + /dilution) # then we multiply by diameter-independent quantities: viral load # and fraction of infected virions From f0bdda7c7d991e35ad705c424b1ae68cca7728c1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 26 Apr 2022 11:19:33 +0200 Subject: [PATCH 10/10] Adding mask inhale efficiency factor to short-range exposure; adding corresponding test; some cosmetics in test_short_range_model.py --- cara/models.py | 2 +- cara/tests/models/test_short_range_model.py | 58 ++++++++++++++++++--- 2 files changed, 51 insertions(+), 9 deletions(-) diff --git a/cara/models.py b/cara/models.py index 9d7898fd..c6df3464 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1377,7 +1377,7 @@ class ExposureModel: f_inf = self.concentration_model.infected.fraction_of_infectious_virus() deposited_exposure *= (f_inf * self.concentration_model.virus.viral_load_in_sputum - ) + * (1 - self.exposed.mask.inhale_efficiency())) # long-range concentration deposited_exposure += self.long_range_deposited_exposure_between_bounds(time1, time2) diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index b1dee050..ceef0bb2 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -6,10 +6,13 @@ import pytest from cara import models import cara.monte_carlo as mc_models from cara.apps.calculator.model_generator import build_expiration -from cara.monte_carlo.data import short_range_expiration_distributions, short_range_distances, activity_distributions +from cara.monte_carlo.data import short_range_expiration_distributions,\ + expiration_distributions, short_range_distances, activity_distributions # TODO: seed better the random number generators np.random.seed(2000) +SAMPLE_SIZE = 250_000 + @pytest.fixture def concentration_model() -> mc_models.ConcentrationModel: @@ -41,8 +44,8 @@ def short_range_model(): def test_short_range_model_ndarray(concentration_model, short_range_model): - concentration_model = concentration_model.build_model(250_000) - model = short_range_model.build_model(250_000) + concentration_model = concentration_model.build_model(SAMPLE_SIZE) + model = short_range_model.build_model(SAMPLE_SIZE) assert isinstance(model._normed_concentration(concentration_model, 10.75), np.ndarray) assert isinstance(model.short_range_concentration(concentration_model, 10.75), np.ndarray) assert isinstance(model._normed_jet_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) @@ -60,10 +63,10 @@ def test_short_range_model_ndarray(concentration_model, short_range_model): ] ) def test_dilution_factor(activity, expected_dilution): - model = models.ShortRangeModel(expiration="Breathing", + model = mc_models.ShortRangeModel(expiration=short_range_expiration_distributions['Breathing'], activity=models.Activity.types[activity], presence=models.SpecificInterval(present_times=((10.5, 11.0),)), - distance=0.854) + distance=0.854).build_model(SAMPLE_SIZE) assert isinstance(model.dilution_factor(), np.ndarray) np.testing.assert_almost_equal( model.dilution_factor(), expected_dilution, decimal=10 @@ -105,11 +108,50 @@ def test_extract_between_bounds(short_range_model, time1, time2, [12.0, 0.], ] ) -def test_short_range_concentration(time, expected_short_range_concentration, concentration_model, short_range_model): - concentration_model = concentration_model.build_model(250_000) - model = short_range_model.build_model(250_000) +def test_short_range_concentration(time, expected_short_range_concentration, + concentration_model, short_range_model): + concentration_model = concentration_model.build_model(SAMPLE_SIZE) + model = short_range_model.build_model(SAMPLE_SIZE) np.testing.assert_allclose( np.array(model.short_range_concentration(concentration_model, time)).mean(), expected_short_range_concentration, rtol=0.02 ) + +def test_short_range_exposure_with_ndarray_mask(): + c_model = mc_models.ConcentrationModel( + room=models.Room(volume=50, humidity=0.3), + ventilation=models.AirChange(active=models.PeriodicInterval(period=120, duration=120), + air_exch=10_000_000,), + infected=mc_models.InfectedPopulation( + number=1, + presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + virus=models.Virus.types['SARS_CoV_2_DELTA'], + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + expiration=expiration_distributions['Breathing'], + host_immunity=0., + ), + evaporation_factor=0.3, + ) + sr_model = mc_models.ShortRangeModel(expiration=short_range_expiration_distributions['Shouting'], + activity=models.Activity.types['Heavy exercise'], + presence=models.SpecificInterval(present_times=((10.5, 11.0),)), + distance=0.854) + e_model = mc_models.ExposureModel( + concentration_model = c_model, + short_range = (sr_model,), + exposed = mc_models.Population( + number=1, + presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + mask=models.Mask(η_inhale=np.array([0., 0.3, 0.5])), + activity=models.Activity.types['Light activity'], + host_immunity=0., + ), + ).build_model(SAMPLE_SIZE) + assert isinstance(e_model.deposited_exposure(), np.ndarray) + assert len(e_model.deposited_exposure()) == 3 + np.testing.assert_allclose(e_model.deposited_exposure(), + e_model.deposited_exposure()[0]*np.array([1., 0.7, 0.5]), + rtol=1e-8) +