From 1760a518de9bfb4f79630e30e236b81ed2215372 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 8 Feb 2022 09:57:16 +0000 Subject: [PATCH] Implemented ShortRangeModel dataclass --- cara/apps/calculator/model_generator.py | 25 ++++++++++++---- cara/apps/expert.py | 3 +- cara/models.py | 39 +++++++++++++++++++++++-- cara/monte_carlo/data.py | 11 ++++--- 4 files changed, 65 insertions(+), 13 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index d0b1aaaa..da43e39c 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -13,7 +13,7 @@ from cara import data import cara.data.weather import cara.monte_carlo as mc from .. import calculator -from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions, initial_concentrations_mouth +from cara.monte_carlo.data import activity_distributions, dilution_factor, virus_distributions, mask_distributions, initial_concentrations_mouth from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions @@ -246,12 +246,29 @@ class FormData: humidity = 0.5 room = models.Room(volume=volume, humidity=humidity) + if self.short_range_option == "short_range_no": + sr_presence=[] + sr_activities=[] + else: + sr_presence=self.short_range_intervals() + sr_activities=self.short_range_activities() + + jet_origin_concentration = [initial_concentrations_mouth[activity] for activity in sr_activities] + + short_range = models.ShortRangeModel( + presence=sr_presence, + activities=sr_activities, + dilutions=dilution_factor(activities=sr_activities, distance=np.random.uniform(0.5, 1.5, 250000)), + jet_origin_concentrations=jet_origin_concentration, + ) + # Initializes and returns a model with the attributes defined above return mc.ExposureModel( concentration_model=mc.ConcentrationModel( room=room, ventilation=self.ventilation(), infected=self.infected_population(), + short_range=short_range, evaporation_factor=0.3, ), exposed=self.exposed_population(), @@ -647,10 +664,8 @@ class FormData: ) def short_range_activities(self) -> typing.List[str]: - if (self.short_range_interactions): - return [interaction['activity'] for interaction in self.short_range_interactions] - else: - return [] + return ([interaction['activity'] for interaction in self.short_range_interactions] + if self.short_range_interactions else []) def build_expiration(expiration_definition) -> models._ExpirationBase: diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 29df1d1c..a7f82243 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -500,9 +500,8 @@ baseline_model = models.ExposureModel( activity=models.Activity.types['Seated'], expiration=models.Expiration.types['Speaking'], host_immunity=0., - short_range_presence=[], - short_range_activities=[], ), + short_range=[], evaporation_factor=0.3, ), exposed=models.Population( diff --git a/cara/models.py b/cara/models.py index f667ba1f..89f3e3d3 100644 --- a/cara/models.py +++ b/cara/models.py @@ -882,11 +882,35 @@ class InfectedPopulation(_PopulationWithVirus): return self.expiration.particle +@dataclass(frozen=True) +class ShortRangeModel: + #: Short range interactions + presence: typing.List[SpecificInterval] + + #: The type of expiractory activities in the short range interactions + activities: typing.List[str] + + #: The dilution factors for each of the expiratory activity in the short range interactions + dilutions: _VectorisedFloat + + #: The concentration on the jet origin for each of the expiratory activity in the short range interactions + jet_origin_concentrations: _VectorisedFloat + + def short_range_concentration(self, time: float, background_concentration: _VectorisedFloat, viral_load: _VectorisedFloat) -> _VectorisedFloat: + for index, period in enumerate(self.presence): + start, finish = tuple(period.boundaries()) + if time >= start and time <= finish: + dilution = self.dilutions[index] + jet_origin_concentration = self.jet_origin_concentrations[index] * 1e-6 * viral_load + return background_concentration + ((1/dilution)*(jet_origin_concentration - background_concentration)) + + @dataclass(frozen=True) class ConcentrationModel: room: Room ventilation: _VentilationBase infected: _PopulationWithVirus + short_range: ShortRangeModel #: evaporation factor: the particles' diameter is multiplied by this # factor as soon as they are in the air (but AFTER going out of the, @@ -1018,8 +1042,13 @@ class ConcentrationModel: Note that time is not vectorised. You can only pass a single float to this method. """ - return (self._normed_concentration(time) * - self.infected.emission_rate_when_present()) + background_concentration = self._normed_concentration(time) * self.infected.emission_rate_when_present() + for period in self.short_range.presence: + start, finish = tuple(period.boundaries()) + if time >= start and time <= finish: + return self.short_range.short_range_concentration(time, background_concentration, self.virus.viral_load_in_sputum) + + return background_concentration @method_cache def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: @@ -1027,6 +1056,12 @@ class ConcentrationModel: Get the integrated concentration of viruses in the air between the times start and stop, normalized by the emission rate. """ + for period in self.short_range.presence: + time1, time2 = tuple(period.boundaries()) + if (time1 <= start <= time2 and time1 <= stop <= time2): + # Check if the given times are within the short range interactions + return scipy.integrate.quad_vec(lambda t: self.concentration(t), start, stop)[0] + if stop <= self._first_presence_time(): return 0.0 state_change_times = self.state_change_times() diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 9dee5d49..ec58d6ee 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -176,8 +176,10 @@ def expiration_distribution(BLO_factors): cn=BLOmodel(BLO_factors).integrate(0.1, 30.)) -def dilution_factor(distance, D=0.02): - u0 = 0.6 +def dilution_factor(activities, distance, D=0.02): + factors = [] + for activity in activities: + u0 = 0.98 if activity == "Breathing" else 3.9 tstar = 2.0 Cr1 = 0.18 Cr2 = 0.2 @@ -192,13 +194,14 @@ def dilution_factor(distance, D=0.02): xstar = Cx1*(Q0*u0)**0.25*(tstar + t01)**0.5 - x01 # Dilution factor at the transition point xstar Sxstar = 2*Cr1*(xstar+x01)/D - return np.mean(np.piecewise(distance, [distance < xstar, distance >= xstar], + factors.append(np.piecewise(distance, [distance < xstar, distance >= xstar], [lambda distance : 2*Cr1*(distance + x01)/D, lambda distance : Sxstar*(1 + Cr2*(distance - xstar)/Cr1/(xstar + x01))**3])) + return factors def initial_concentration_mouth(BLO_factors): value, error = scipy.integrate.quad(lambda d: BLOmodel(BLO_factors).distribution(d) * BLOmodel(BLO_factors).volume(d), 0.1, 1000) - return value * 1e-6 #result in mL/m^3 + return value expiration_BLO_factors = {