diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index d8a7ce9b..d72952e6 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, dilution_factor +from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions, short_range_distances from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions, short_range_expiration_distributions @@ -245,29 +245,31 @@ class FormData: else: humidity = 0.5 room = models.Room(volume=volume, humidity=humidity) + + infected_population = self.infected_population() if self.short_range_option == "short_range_yes": + short_range_expirations = tuple(short_range_expiration_distributions[interaction['expiration']] for interaction in self.short_range_interactions) + short_range_activities = tuple([infected_population.activity for _ in self.short_range_interactions]) sr_presence=self.short_range_intervals() - sr_activities=self.short_range_activities() - short_range_expirations = tuple(short_range_expiration_distributions[activity] for activity in sr_activities) - dilutions=dilution_factor(activities=sr_activities) else: - sr_presence=() short_range_expirations=() - dilutions=() - + short_range_activities=() + sr_presence=() + # 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(), + infected=infected_population, evaporation_factor=0.3, ), short_range = mc.ShortRangeModel( - presence=sr_presence, expirations=short_range_expirations, - dilutions=dilutions, + activities=short_range_activities, + presence=sr_presence, + distances=short_range_distances, ), exposed=self.exposed_population(), ) @@ -655,7 +657,7 @@ class FormData: for interaction in self.short_range_interactions: start_time = time_string_to_minutes(interaction['start_time']) duration = float(interaction['duration']) - short_range_intervals.append((interaction['activity'], models.SpecificInterval((start_time/60, (start_time + duration)/60)))) + short_range_intervals.append((interaction['expiration'], models.SpecificInterval((start_time/60, (start_time + duration)/60)))) return tuple(short_range_intervals) @@ -665,10 +667,6 @@ class FormData: breaks=self.exposed_lunch_break_times() + self.exposed_coffee_break_times(), ) - def short_range_activities(self) -> typing.List[str]: - return ([interaction['activity'] for interaction in self.short_range_interactions] - if self.short_range_interactions else []) - def build_expiration(expiration_definition) -> mc._ExpirationBase: if isinstance(expiration_definition, str): diff --git a/cara/apps/calculator/static/js/form.js b/cara/apps/calculator/static/js/form.js index dcaf0662..545f6f99 100644 --- a/cara/apps/calculator/static/js/form.js +++ b/cara/apps/calculator/static/js/form.js @@ -408,7 +408,7 @@ function validate_form(form) { $(".form_field_outer_row").each(function (index, element){ let obj = {}; const $element = $(element); - obj.activity = $element.find("[name='short_range_activity']").val(); + obj.expiration = $element.find("[name='short_range_expiration']").val(); obj.start_time = $element.find("[name='short_range_start_time']").val(); obj.duration = $element.find("[name='short_range_duration']").val(); short_range_interactions.push(JSON.stringify(obj)); @@ -656,8 +656,8 @@ $(document).ready(function () { let index = 1; for (const interaction of JSON.parse(value)) { $("#dialog_sr").append(inject_sr_interaction(index, value = interaction, is_validated="row_validated")) - $('#sr_activity_no_' + String(index)).val(interaction.activity).change(); - document.getElementById('sr_activity_no_' + String(index)).disabled = true; + $('#sr_expiration_no_' + String(index)).val(interaction.expiration).change(); + document.getElementById('sr_expiration_no_' + String(index)).disabled = true; document.getElementById('sr_start_no_' + String(index)).disabled = true; document.getElementById('sr_duration_no_' + String(index)).disabled = true; document.getElementById('edit_row_no_' + String(index)).style.cssText = 'display:inline !important'; @@ -838,35 +838,35 @@ $(document).ready(function () { function inject_sr_interaction(index, value, is_validated) { return `
-
-
+
+
-

-
-
+
-
-
-

-
-
-
-
-
-

+ +
+
+

-
+ +
+
+

+
+ +
-
` } @@ -889,12 +889,12 @@ $(document).ready(function () { // Validate row button (Save button) $("body").on("click", ".validate_node_btn_frm_field", function() { var index = $(this).attr('id').split('_').slice(-1)[0]; - let activity = validate_sr_parameter('#sr_activity_no_' + String(index)[0], "Required input."); + let activity = validate_sr_parameter('#sr_expiration_no_' + String(index)[0], "Required input."); let start = validate_sr_parameter('#sr_start_no_' + String(index)[0], "Required input."); let duration = validate_sr_parameter('#sr_duration_no_' + String(index)[0], "Required input."); if (activity && start && duration) { if (validate_sr_time('#sr_start_no_' + String(index)) && validate_sr_time('#sr_start_no_' + String(index))) { - document.getElementById('sr_activity_no_' + String(index)).disabled = true; + document.getElementById('sr_expiration_no_' + String(index)).disabled = true; document.getElementById('sr_start_no_' + String(index)).disabled = true; document.getElementById('sr_duration_no_' + String(index)).disabled = true; document.getElementById('edit_row_no_' + String(index)).style.cssText = 'display:inline !important'; @@ -915,7 +915,7 @@ $(document).ready(function () { $("body").on("click", ".edit_node_btn_frm_field", function() { $(this).hide(); let id = $(this).attr('id').split('_').slice(-1)[0]; - document.getElementById('sr_activity_no_' + String(id)).disabled = false; + document.getElementById('sr_expiration_no_' + String(id)).disabled = false; document.getElementById('sr_start_no_' + String(id)).disabled = false; document.getElementById('sr_duration_no_' + String(id)).disabled = false; document.getElementById('validate_row_no_' + String(id)).style.cssText = 'display:inline !important'; diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 208b3568..9fab2cd9 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -504,9 +504,10 @@ baseline_model = models.ExposureModel( evaporation_factor=0.3, ), short_range=models.ShortRangeModel( - presence=(), expirations=(), - dilutions=(), + activities=(), + presence=(), + distances=(), ), exposed=models.Population( number=10, diff --git a/cara/apps/templates/base/calculator.form.html.j2 b/cara/apps/templates/base/calculator.form.html.j2 index 92d0d42f..499fc761 100644 --- a/cara/apps/templates/base/calculator.form.html.j2 +++ b/cara/apps/templates/base/calculator.form.html.j2 @@ -391,11 +391,12 @@
-
-
+ +
+
diff --git a/cara/apps/templates/base/calculator.report.html.j2 b/cara/apps/templates/base/calculator.report.html.j2 index 43126bce..d35bf86a 100644 --- a/cara/apps/templates/base/calculator.report.html.j2 +++ b/cara/apps/templates/base/calculator.report.html.j2 @@ -335,7 +335,7 @@

    {% for interaction in form.short_range_interactions %} -
  • Expiratory activity {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.activity }}
    • +
  • Expiratory activity {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.expiration }}
  • Start time {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.start_time }}
  • Duration {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.duration }} {{ "minutes" if interaction.duration|float > 1 else "minute" }}
    • {% endfor %} diff --git a/cara/models.py b/cara/models.py index 5b7a1e1d..fd658861 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1074,14 +1074,52 @@ class ConcentrationModel: @dataclass(frozen=True) class ShortRangeModel: - #: Short-range activities and respective interactions - presence: typing.Tuple[typing.Tuple[str, SpecificInterval], ...] - #: Expiration types expirations: typing.Tuple[_ExpirationBase, ...] - #: The dilution factors for each of the expiratory activity - dilutions: typing.Tuple[_VectorisedFloat, ...] + #: Activity types + activities: typing.Tuple[Activity, ...] + + #: Short-range expirtory activities and respective interactions + presence: typing.Tuple[typing.Tuple[str, SpecificInterval], ...] + + #: Interpersonal distances + distances: _VectorisedFloat + + def dilution_factor(self, BR: _VectorisedFloat) -> _VectorisedFloat: + ''' + The dilution factor for the respective expiratory activity type. + ''' + # Average mouth diameter + D = 0.02 + # Convert Breathing rate from m3/h to m3/s + BR_s = BR/3600 + # Area of the mouth assuming a perfect circle + Am = np.pi*(D**2)/4 + # Initial velocity from the division of the Breathing rate with the area + u0 = BR_s/Am + + 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 = np.mean(Cx1*(Q0*u0)**0.25*(tstar + t01)**0.5 - x01) + # Dilution factor at the transition point xstar + Sxstar = np.mean(2*Cr1*(xstar+x01)/D) + + return np.piecewise(self.distances, [self.distances < xstar, self.distances >= xstar], + [ + lambda distance: 2*Cr1*(distance + x01)/D, + lambda distance: Sxstar*(1 + Cr2*(distance - xstar)/Cr1/(xstar + x01))**3 + ]) def _normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: """ @@ -1095,16 +1133,18 @@ class ShortRangeModel: start, stop = tuple(period.boundaries()) # Verifies if the given time falls within a short-range interaction if start < time <= stop: - dilution = self.dilutions[index] + dilution = self.dilution_factor(self.activities[index].exhalation_rate) jet_origin_concentration = self.expirations[index].jet_origin_concentration() # Long-range concentration normalized by the virus viral load - long_range_normed_concentration = concentration_model.concentration(time) / concentration_model.virus.viral_load_in_sputum + long_range_normed_concentration = (concentration_model.concentration(time) / + concentration_model.virus.viral_load_in_sputum) # The long-range concentration values are then approximated using interpolation: # The set of points where we want the interpolated values are the short-range particle diameters (given the current expiration); # The set of points with a known value are the long-range particle diameters (given the initial expiration); # The set of known values are the long-range concentration values normalized by the viral load. - long_range_normed_concentration_interpolated=np.interp(self.expirations[index].particle.diameter, concentration_model.infected.particle.diameter, long_range_normed_concentration) + long_range_normed_concentration_interpolated=np.interp(self.expirations[index].particle.diameter, + concentration_model.infected.particle.diameter, long_range_normed_concentration) # Short-range concentration formula. The long-range concentration is added in the concentration method (ExposureModel). return ((1/dilution)*(jet_origin_concentration - long_range_normed_concentration_interpolated)) @@ -1146,7 +1186,7 @@ class ShortRangeModel: start_bound, stop_bound = start, stop jet_origin_integrated = self.expirations[index].jet_origin_concentration() - dilution = self.dilutions[index] + dilution = self.dilution_factor(self.activities[index].exhalation_rate) 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) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 98095389..58700c7b 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -182,39 +182,39 @@ def expiration_distribution( ) -def dilution_factor(activities): - D = 0.02 - # From https://www.mdpi.com/1660-4601/17/4/1445/htm - distance = LogNormal(0.8542127255693238, 0.42755967248106513).generate_samples(250_000) +# def dilution_factor(activities): +# D = 0.02 +# # Fit from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm +# distance = LogNormal(-0.269359136417347, 0.4728300188814934).generate_samples(250_000) - factors = [] - for activity in activities: - u0 = 0.98 if activity == "Breathing" else 3.9 - 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 - 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 +# factors = [] +# for activity in activities: +# u0 = 0.98 if activity == "Breathing" else 3.9 +# 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 +# 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 expiration_BLO_factors = { @@ -235,3 +235,6 @@ short_range_expiration_distributions = { exp_type: expiration_distribution(BLO_factors, d_max=100) for exp_type, BLO_factors in expiration_BLO_factors.items() } + + +short_range_distances = LogNormal(-0.269359136417347, 0.4728300188814934)