From 7e1dfa2684478f3293a9e2ff555dc8f5bcf175d9 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 16 Mar 2022 16:58:30 +0000 Subject: [PATCH] Replaced scipy method with analytical method and updated short range tests --- cara/apps/calculator/model_generator.py | 5 +-- cara/apps/calculator/report_generator.py | 6 ++- cara/apps/expert.py | 1 - cara/models.py | 48 +++++++++++++++------ cara/tests/conftest.py | 1 - cara/tests/models/test_exposure_model.py | 1 - cara/tests/models/test_short_range_model.py | 43 ++++++++++++------ cara/tests/test_monte_carlo.py | 1 - cara/tests/test_monte_carlo_full_models.py | 1 - 9 files changed, 70 insertions(+), 37 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index ea69b39c..d4147be2 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -259,7 +259,6 @@ class FormData: evaporation_factor=0.3, ), short_range = mc.ShortRangeModel( - activities=sr_activities, presence=sr_presence, expirations=short_range_expirations, dilutions=dilution_factor(activities=sr_activities, distance=np.random.uniform(0.5, 1.5, 250000)), @@ -639,13 +638,13 @@ class FormData: breaks=self.infected_lunch_break_times() + self.infected_coffee_break_times(), ) - def short_range_intervals(self) -> typing.Tuple[models.SpecificInterval, ...]: + def short_range_intervals(self) -> typing.Tuple[typing.Tuple[str, models.SpecificInterval], ...]: if (self.short_range_interactions): short_range_intervals = [] for interaction in self.short_range_interactions: start_time = time_string_to_minutes(interaction['start_time']) duration = float(interaction['duration']) - short_range_intervals.append(models.SpecificInterval((start_time/60, (start_time + duration)/60))) + short_range_intervals.append((interaction['activity'], models.SpecificInterval((start_time/60, (start_time + duration)/60)))) return tuple(short_range_intervals) else: return () diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 6fecaf98..6f35c88a 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -98,8 +98,10 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L def calculate_report_data(model: models.ExposureModel): times = interesting_times(model) + short_range_activities = [] short_range_intervals = [] - for interval in model.short_range.presence: + for (activity, interval) in model.short_range.presence: + short_range_activities.append(activity) short_range_intervals.append(list(interval.boundaries())) short_range_concentrations = [ @@ -111,7 +113,7 @@ def calculate_report_data(model: models.ExposureModel): if len(short_range_intervals) != 0: for index, (start, stop) in enumerate(short_range_intervals): # For visualization issues, add short range breathing activity to the initial long range concentrations - if model.short_range.activities[index] == 'Breathing': + if short_range_activities[index] == 'Breathing': sr_breathing_concentrations.append(np.array(model.concentration(float(stop))).mean()) break diff --git a/cara/apps/expert.py b/cara/apps/expert.py index b0a7d0f5..208b3568 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -504,7 +504,6 @@ baseline_model = models.ExposureModel( evaporation_factor=0.3, ), short_range=models.ShortRangeModel( - activities=(), presence=(), expirations=(), dilutions=(), diff --git a/cara/models.py b/cara/models.py index a56322c2..9f1ac43d 100644 --- a/cara/models.py +++ b/cara/models.py @@ -37,7 +37,6 @@ import typing import numpy as np from scipy.interpolate import interp1d -import scipy.integrate if not typing.TYPE_CHECKING: from memoization import cached @@ -1075,11 +1074,8 @@ class ConcentrationModel: @dataclass(frozen=True) class ShortRangeModel: - #: Short range activities - activities: typing.Tuple[str, ...] - - #: Short range interactions - presence: typing.Tuple[SpecificInterval, ...] + #: Short range activities and respective interactions + presence: typing.Tuple[typing.Tuple[str, SpecificInterval], ...] #: Expiration types expirations: typing.Tuple[_ExpirationBase, ...] @@ -1095,7 +1091,7 @@ class ShortRangeModel: short range concentration normalized by the virus viral load. Otherwise it returns 0. """ - for index, period in enumerate(self.presence): + for index, (activity, period) in enumerate(self.presence): start, stop = tuple(period.boundaries()) # Verifies if the given time falls within a short range interaction if start < time <= stop: @@ -1122,15 +1118,41 @@ class ShortRangeModel: return (self._normed_concentration(concentration_model, time) * concentration_model.virus.viral_load_in_sputum) - def normed_exposure_between_bounds(self, concentration_model, time1, time2): + @method_cache + def _normed_short_range_concentration_cached(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: + # A cached version of the _normed_concentration method. Use this + # method if you expect that there may be multiple short range concentration + # 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): """ Get the integrated short range concentration of viruses in the air between the times start and stop, normalized by the virus viral load. """ - # TODO Implement the integration without using scipy. - return scipy.integrate.quad_vec(lambda t: self._normed_concentration(concentration_model, t), time1, time2)[0] + for index, (activity, period) in enumerate(self.presence): + start, stop = tuple(period.boundaries()) + 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 + + jet_origin_integrated = self.expirations[index].jet_origin_concentration() + dilution = self.dilutions[index] + + total_normed_concentration = -(concentration_model.integrated_concentration(start_bound, stop_bound)/concentration_model.virus.viral_load_in_sputum/dilution) + total_normed_concentration_interpolated = np.interp(self.expirations[index].particle.diameter, concentration_model.infected.particle.diameter, total_normed_concentration) + return (jet_origin_integrated/dilution * (stop_bound - start_bound)) + total_normed_concentration_interpolated + - @dataclass(frozen=True) class ExposureModel: """ @@ -1201,7 +1223,7 @@ class ExposureModel: initial deposited exposure. """ deposited_exposure = 0. - for index, period in enumerate(self.short_range.presence): + for index, (activity, period) in enumerate(self.short_range.presence): start, stop = tuple(period.boundaries()) if stop < time1: continue @@ -1217,7 +1239,7 @@ class ExposureModel: start_bound, stop_bound = start, stop short_range_exposure = self.short_range.normed_exposure_between_bounds(self.concentration_model, start_bound, stop_bound) - fdep = self.short_range.expirations[index].particle.fraction_deposited(evaporation_factor=1.0) #ASK for evaporation factor + fdep = self.short_range.expirations[index].particle.fraction_deposited(evaporation_factor=1.0) diameter = self.short_range.expirations[index].particle.diameter # Aerosols not considered given the formula for the initial concentration at mouth/nose. diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py index 1784c383..3b7b4b88 100644 --- a/cara/tests/conftest.py +++ b/cara/tests/conftest.py @@ -32,7 +32,6 @@ def baseline_concentration_model(): @pytest.fixture def baseline_sr_model(): return models.ShortRangeModel( - activities=(), presence=(), expirations=(), dilutions=(), diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 5b3eeb55..f5a98fbd 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -172,7 +172,6 @@ def conc_model(): @pytest.fixture def sr_model(): return models.ShortRangeModel( - activities=(), presence=(), expirations=(), dilutions=(), diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index 2cee767b..3e40cc56 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -30,39 +30,54 @@ def concentration_model() -> mc_models.ConcentrationModel: ) -activities = ['Breathing', 'Speaking', 'Shouting'] - - @pytest.fixture -def presences() -> typing.Tuple[models.SpecificInterval, ...]: +def presences() -> typing.Tuple[typing.Tuple[str, models.SpecificInterval], ...]: return ( - models.SpecificInterval((10.5, 11.0)), - models.SpecificInterval((14.5, 15.0)), - models.SpecificInterval((16.5, 17.5)), + ('Breathing', models.SpecificInterval((10.5, 11.0))), + ('Speaking', models.SpecificInterval((14.5, 15.0))), + ('Shouting', models.SpecificInterval((16.5, 17.5))), ) @pytest.fixture -def expirations() -> typing.Tuple[mc_models.Expiration, ...]: +def expirations(presences) -> typing.Tuple[mc_models.Expiration, ...]: # Monte carlo expirations! So build model. - return tuple(short_range_expiration_distributions[activity] for activity in activities) + return tuple(short_range_expiration_distributions[activity] for (activity, presence) in presences) @pytest.fixture -def dilutions(): +def dilutions(presences): return dilution_factor( - activities=activities, + activities=[activity for (activity, presence) in presences], distance=np.random.uniform(0.5, 1.5, 250_000), ) def test_short_range_model_ndarray(concentration_model, presences, expirations, dilutions): concentration_model = concentration_model.build_model(250_000) - model = mc_models.ShortRangeModel(activities, presences, expirations, dilutions) + model = mc_models.ShortRangeModel(presences, expirations, dilutions) model = model.build_model(250_000) assert isinstance(model._normed_concentration(concentration_model, 10.75), np.ndarray) assert isinstance(model.short_range_concentration(concentration_model, 14.75), np.ndarray) - assert isinstance(model.normed_exposure_between_bounds(concentration_model, 16.6, 17.7), np.ndarray) + assert isinstance(model.normed_exposure_between_bounds(concentration_model, 16.6, 17.5), np.ndarray) + + +@pytest.mark.parametrize( + "start, stop, expected_exposure", [ + [8.5, 12.5, 5.844666077067048e-09], + [10.5, 11.0, 5.830120846251791e-09], + [10.6, 11.9, 4.6397748633454945e-09], + ] +) +def test_normed_exposure_between_bounds( + start, stop, expected_exposure, concentration_model, + presences, expirations, dilutions): + concentration_model = concentration_model.build_model(250_000) + model = mc_models.ShortRangeModel(presences, expirations, dilutions) + model = model.build_model(250_000) + np.testing.assert_almost_equal( + model.normed_exposure_between_bounds(concentration_model, start, stop).mean(), expected_exposure + ) @pytest.mark.parametrize( @@ -77,7 +92,7 @@ def test_short_range_model( concentration_model, presences, expirations, dilutions, ): concentration_model = concentration_model.build_model(250_000) - model = mc_models.ShortRangeModel(activities, presences, expirations, dilutions) + model = mc_models.ShortRangeModel(presences, expirations, dilutions) model = model.build_model(250_000) np.testing.assert_almost_equal( model._normed_concentration(concentration_model, time).mean(), expected_sr_normed_concentration, decimal=0 diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index e7a70120..0e4d2bb0 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -64,7 +64,6 @@ def baseline_mc_concentration_model() -> cara.monte_carlo.ConcentrationModel: @pytest.fixture def baseline_mc_sr_model() -> cara.monte_carlo.ShortRangeModel: return cara.monte_carlo.ShortRangeModel( - activities=(), presence=(), expirations=(), dilutions=(), diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index a7338ca5..26911a5f 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -43,7 +43,6 @@ TorontoTemperatures = { @pytest.fixture def sr_model_mc() -> mc.ShortRangeModel: return mc.ShortRangeModel( - activities=(), presence=(), expirations=(), dilutions=(),