From 1ea01c12745c6190c62e8fab3d07213281fc6797 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 21 Jan 2022 17:11:05 +0100 Subject: [PATCH] Short range presence and short range activities as part of the backend model --- cara/apps/calculator/model_generator.py | 50 +++++++++++++++++++++++- cara/apps/calculator/report_generator.py | 36 ++++++++++++++++- cara/apps/expert.py | 3 ++ cara/models.py | 6 +++ cara/monte_carlo/data.py | 37 +++++++++++++++++- 5 files changed, 128 insertions(+), 4 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 595094bc..4ded107d 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -3,6 +3,8 @@ import datetime import html import logging import typing +import ast +import json import numpy as np @@ -11,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 +from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions, initial_concentrations_mouth from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions @@ -76,6 +78,8 @@ class FormData: window_width: float windows_number: int window_opening_regime: str + short_range_option: str + short_range_interactions: list #: The default values for undefined fields. Note that the defaults here #: and the defaults in the html form must not be contradictory. @@ -127,6 +131,8 @@ class FormData: 'windows_frequency': 0., 'windows_number': 0, 'window_opening_regime': 'windows_open_permanently', + 'short_range_option': False, + 'short_range_interactions': [], } @classmethod @@ -416,6 +422,8 @@ class FormData: number=infected_occupants, virus=virus, presence=self.infected_present_interval(), + short_range_presence=self.short_range_intervals(), + short_range_activities=self.short_range_activities(), mask=self.mask(), activity=activity, expiration=expiration, @@ -448,6 +456,7 @@ class FormData: exposed = mc.Population( number=exposed_occupants, presence=self.exposed_present_interval(), + short_range_presence=self.short_range_interactions, activity=activity, mask=self.mask(), host_immunity=0., @@ -623,12 +632,30 @@ class FormData: breaks=self.infected_lunch_break_times() + self.infected_coffee_break_times(), ) + def short_range_intervals(self) -> models.Interval: + 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))) + + return short_range_intervals + else: + return [] + def exposed_present_interval(self) -> models.Interval: return self.present_interval( self.exposed_start, self.exposed_finish, breaks=self.exposed_lunch_break_times() + self.exposed_coffee_break_times(), ) + def short_range_activities(self): + if (self.short_range_interactions): + return [interaction['activity'] for interaction in self.short_range_interactions] + else: + return [] + def build_expiration(expiration_definition) -> models._ExpirationBase: if isinstance(expiration_definition, str): @@ -642,6 +669,12 @@ def build_expiration(expiration_definition) -> models._ExpirationBase: return expiration_distribution(tuple(BLO_factors)) +def build_concentration_mouth(short_range_def) -> models._ExpirationBase: + expiration_definition = short_range_def[1]['activity'] + if isinstance(expiration_definition, str): + return initial_concentrations_mouth[expiration_definition] + + def baseline_raw_form_data(): # Note: This isn't a special "baseline". It can be updated as required. return { @@ -691,7 +724,8 @@ def baseline_raw_form_data(): 'window_type': 'window_sliding', 'window_width': '2', 'windows_number': '1', - 'window_opening_regime': 'windows_open_permanently' + 'window_opening_regime': 'windows_open_permanently', + 'short_range_option': 'short_range_no', } @@ -736,6 +770,15 @@ def time_minutes_to_string(time: int) -> str: return "{0:0=2d}".format(int(time/60)) + ":" + "{0:0=2d}".format(time%60) +def string_to_list(l: str) -> list: + if (l != []): + return list(ast.literal_eval(l.replace(""", "\""))) + + +def list_to_string(s: list) -> str: + return json.dumps(s) + + def _safe_int_cast(value) -> int: if isinstance(value, int): return value @@ -767,3 +810,6 @@ for _field in dataclasses.fields(FormData): elif _field.type is bool: _CAST_RULES_FORM_ARG_TO_NATIVE[_field.name] = lambda v: v == '1' _CAST_RULES_NATIVE_TO_FORM_ARG[_field.name] = int + elif _field.type is list: + _CAST_RULES_FORM_ARG_TO_NATIVE[_field.name] = string_to_list + _CAST_RULES_NATIVE_TO_FORM_ARG[_field.name] = list_to_string diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index e703e225..ec8230c2 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -12,8 +12,9 @@ import jinja2 import numpy as np from cara import models +from cara.monte_carlo.data import initial_concentrations_mouth, dilution_factor from ... import monte_carlo as mc -from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE +from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE, build_expiration from ... import dataclass_utils @@ -96,13 +97,45 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L return nice_times +def short_range_interesting_times(model: models.ExposureModel, times: typing.List[float]) -> typing.List[float]: + short_range_times = [] + for period in model.concentration_model.infected.short_range_presence: + start, finish = period.boundaries() + short_range_times = short_range_times + [time for time in times if time > start and time < finish] + return short_range_times + + +def jet_origin_concentrations(model: models.ExposureModel) -> models._ExpirationBase: + return [initial_concentrations_mouth[activity] for activity in model.concentration_model.infected.short_range_activities] + + +def short_range_initial_concentrations(model: models.ExposureModel, time: float): + dilution = dilution_factor(np.linspace(0.1, 2., 1000)) + jet_origin_initial_concentrations = jet_origin_concentrations(model) + for index, interaction in enumerate(model.concentration_model.infected.short_range_presence): + start, finish = interaction.boundaries() + if start <= time <= finish: + expiration = build_expiration(model.concentration_model.infected.short_range_activities[index]) + single_exposure_model = dataclass_utils.nested_replace( + model, {'concentration_model.infected.expiration': expiration} + ) + concentration = single_exposure_model.concentration_model._normed_concentration(float(time)) + jet_origin_concentration = jet_origin_initial_concentrations[index] * model.concentration_model.infected.virus.viral_load_in_sputum + return concentration + ((1/dilution)*(jet_origin_concentration - concentration)) + + def calculate_report_data(model: models.ExposureModel): times = interesting_times(model) + times_short_range = short_range_interesting_times(model, times) concentrations = [ np.array(model.concentration_model.concentration(float(time))).mean() for time in times ] + short_range_concentrations = [ + np.array(short_range_initial_concentrations(model, float(time))).mean() + for time in times_short_range + ] highest_const = max(concentrations) prob = np.array(model.infection_probability()).mean() er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean() @@ -118,6 +151,7 @@ def calculate_report_data(model: models.ExposureModel): "exposed_presence_intervals": [list(interval) for interval in model.exposed.presence.boundaries()], "cumulative_doses": list(cumulative_doses), "concentrations": concentrations, + "short_range_concentrations": short_range_concentrations, "highest_const": highest_const, "prob_inf": prob, "emission_rate": er, diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 261afe75..e779558a 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -496,6 +496,8 @@ baseline_model = models.ExposureModel( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((8., 12.), (13., 17.))), + short_range_presence=[], + short_range_activities=[], mask=models.Mask.types['No mask'], activity=models.Activity.types['Seated'], expiration=models.Expiration.types['Speaking'], @@ -506,6 +508,7 @@ baseline_model = models.ExposureModel( exposed=models.Population( number=10, presence=models.SpecificInterval(((8., 12.), (13., 17.))), + short_range_presence=[], activity=models.Activity.types['Seated'], mask=models.Mask.types['No mask'], host_immunity=0., diff --git a/cara/models.py b/cara/models.py index 3cedbf22..3092a6ea 100644 --- a/cara/models.py +++ b/cara/models.py @@ -678,6 +678,9 @@ class Population: #: The times in which the people are in the room. presence: Interval + #: Short range interactions + short_range_presence: list + #: The kind of mask being worn by the people. mask: Mask @@ -749,6 +752,9 @@ class EmittingPopulation(_PopulationWithVirus): class InfectedPopulation(_PopulationWithVirus): #: The type of expiration that is being emitted whilst doing the activity. expiration: _ExpirationBase + + #: The type of expiractory activities in the short range interactions + short_range_activities: list @method_cache def fraction_of_infectious_virus(self) -> _VectorisedFloat: diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 3b9172d2..c4e1dfb2 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -3,6 +3,7 @@ import typing import numpy as np from scipy import special as sp +import scipy.integrate import cara.monte_carlo as mc from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel,CustomKernel,Uniform @@ -51,6 +52,9 @@ class BLOmodel: np.exp(-(np.log(d) - mu) ** 2 / (2 * sigma ** 2)) for A,cn,mu,sigma in zip(self.BLO_factors, self.cn, self.mu, self.sigma) ) + + def volume(self, d): + return(np.pi * d**3) / 6 def integrate(self, dmin, dmax): """ @@ -162,7 +166,7 @@ def expiration_distribution(BLO_factors): Returns an Expiration with an aerosol diameter distribution, defined by the BLO factors (a length-3 tuple). The total concentration of aerosols is computed by integrating - the distribution between 0.1 and 30 microns - these boundaries are + the distribution between 0.1 and D microns - these boundaries are an historical choice based on previous implementations of the model (it limits the influence of the O-mode). """ @@ -172,6 +176,31 @@ def expiration_distribution(BLO_factors): BLOmodel(BLO_factors).integrate(0.1, 30.)) +def dilution_factor(distance, D=0.02): + u0 = 0.6 + tstar = 2.0 + Cr1 = 0.18 + Cr2 = 0.2 + Cx1 = 2.4 + # The expired flow rate during the expiration period, m^3/s + Q0 = u0 * np.pi/4*D**2 + # Parameters in the jet-like stage + x01 = D/2/Cr1 + # Time of virtual origin + t01 = (x01/Cx1)**2 * (Q0*u0)**(-0.5) + # The transition point, m + 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], + [lambda distance : 2*Cr1*(distance + x01)/D, lambda distance : Sxstar*(1 + Cr2*(distance - xstar)/Cr1/(xstar + x01))**3])) + + +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 + + expiration_BLO_factors = { 'Breathing': (1., 0., 0.), 'Speaking': (1., 1., 1.), @@ -184,3 +213,9 @@ expiration_distributions = { exp_type: expiration_distribution(BLO_factors) for exp_type,BLO_factors in expiration_BLO_factors.items() } + + +initial_concentrations_mouth = { + exp_type: initial_concentration_mouth(BLO_factors) + for exp_type,BLO_factors in expiration_BLO_factors.items() +}