diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index d72952e6..5131c56d 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -248,14 +248,15 @@ class FormData: infected_population = self.infected_population() + short_range = [] 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() - else: - short_range_expirations=() - short_range_activities=() - sr_presence=() + for interaction in self.short_range_interactions: + short_range.append(mc.ShortRangeModel( + expiration=short_range_expiration_distributions[interaction['expiration']], + activity=infected_population.activity, + presence=self.short_range_interval(interaction), + distances=short_range_distances, + )) # Initializes and returns a model with the attributes defined above return mc.ExposureModel( @@ -265,12 +266,7 @@ class FormData: infected=infected_population, evaporation_factor=0.3, ), - short_range = mc.ShortRangeModel( - expirations=short_range_expirations, - activities=short_range_activities, - presence=sr_presence, - distances=short_range_distances, - ), + short_range = tuple(short_range), exposed=self.exposed_population(), ) @@ -652,14 +648,10 @@ class FormData: breaks=self.infected_lunch_break_times() + self.infected_coffee_break_times(), ) - def short_range_intervals(self) -> typing.Tuple[typing.Tuple[str, models.SpecificInterval], ...]: - 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((interaction['expiration'], models.SpecificInterval((start_time/60, (start_time + duration)/60)))) - return tuple(short_range_intervals) - + def short_range_interval(self, interaction) -> models.SpecificInterval: + start_time = time_string_to_minutes(interaction['start_time']) + duration = float(interaction['duration']) + return models.SpecificInterval(present_times=((start_time/60, (start_time + duration)/60),),) def exposed_present_interval(self) -> models.Interval: return self.present_interval( diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index f2926b0b..b526d5bf 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -97,37 +97,28 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L return nice_times -def interesting_sr_interactions(model: models.ExposureModel) -> typing.Tuple[typing.List, typing.List]: - short_range_activities = [] - short_range_intervals = [] - for (activity, interval) in model.short_range.presence: - short_range_activities.append(activity) - short_range_intervals.append(list(interval.boundaries())) - return short_range_activities, short_range_intervals - - -def concentrations_with_sr_breathing(model: models.ExposureModel, times: typing.List[float], short_range_intervals: typing.List, - short_range_activities:typing.List) -> typing.List[float]: +def concentrations_with_sr_breathing(form: FormData, model: models.ExposureModel, times: typing.List[float], short_range_intervals: typing.List) -> typing.List[float]: lower_concentrations = [] for time in times: for index, (start, stop) in enumerate(short_range_intervals): # For visualization issues, add short-range breathing activity to the initial long-range concentrations - if start <= time <= stop and short_range_activities[index] == 'Breathing': + if start <= time <= stop and form.short_range_interactions[index]['expiration'] == 'Breathing': lower_concentrations.append(np.array(model.concentration(float(time))).mean()) break lower_concentrations.append(np.array(model.concentration_model.concentration(float(time))).mean()) return lower_concentrations -def calculate_report_data(model: models.ExposureModel): +def calculate_report_data(form: FormData, model: models.ExposureModel): times = interesting_times(model) - short_range_activities, short_range_intervals = interesting_sr_interactions(model) + short_range_intervals = [interaction.presence.boundaries()[0] for interaction in model.short_range] + short_range_expirations = [interaction['expiration'] for interaction in form.short_range_interactions] if form.short_range_option == "short_range_yes" else [] concentrations = [ np.array(model.concentration(float(time))).mean() for time in times ] - lower_concentrations = concentrations_with_sr_breathing(model, times, short_range_intervals, short_range_activities) + lower_concentrations = concentrations_with_sr_breathing(form, model, times, short_range_intervals) highest_const = max(concentrations) cumulative_doses = np.cumsum([ @@ -148,7 +139,7 @@ def calculate_report_data(model: models.ExposureModel): "times": list(times), "exposed_presence_intervals": [list(interval) for interval in model.exposed.presence.boundaries()], "short_range_intervals": short_range_intervals, - "short_range_activities": short_range_activities, + "short_range_expirations": short_range_expirations, "concentrations": concentrations, "concentrations_zoomed": lower_concentrations, "highest_const": highest_const, @@ -272,9 +263,6 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp def scenario_statistics(mc_model: mc.ExposureModel, sample_times: typing.List[float]): model = mc_model.build_model(size=_DEFAULT_MC_SAMPLE_SIZE) - short_range_activities, short_range_intervals = interesting_sr_interactions(model) - lower_concentrations = concentrations_with_sr_breathing(model, sample_times, short_range_intervals, short_range_activities) - return { 'probability_of_infection': np.mean(model.infection_probability()), 'expected_new_cases': np.mean(model.expected_new_cases()), @@ -282,7 +270,6 @@ def scenario_statistics(mc_model: mc.ExposureModel, sample_times: typing.List[fl np.mean(model.concentration(time)) for time in sample_times ], - 'concentrations_zoomed': lower_concentrations, } @@ -340,7 +327,7 @@ class ReportGenerator: scenario_sample_times = interesting_times(model) - context.update(calculate_report_data(model)) + context.update(calculate_report_data(form, model)) alternative_scenarios = manufacture_alternative_scenarios(form) context['alternative_scenarios'] = comparison_report( alternative_scenarios, scenario_sample_times, executor_factory=executor_factory, diff --git a/cara/apps/calculator/static/js/report.js b/cara/apps/calculator/static/js/report.js index a7c59bcc..56ba86bc 100644 --- a/cara/apps/calculator/static/js/report.js +++ b/cara/apps/calculator/static/js/report.js @@ -5,7 +5,7 @@ function draw_plot(svg_id) { let button_full_exposure = document.getElementById("button_full_exposure"); let button_hide_high_concentration = document.getElementById("button_hide_high_concentration"); let long_range_checkbox = document.getElementById('long_range_cumulative_checkbox') - let show_sr_legend = short_range_activities.length > 0; + let show_sr_legend = short_range_expirations.length > 0; var data_for_graphs = { 'concentrations': [], @@ -101,7 +101,7 @@ function draw_plot(svg_id) { .text('Presence of exposed person(s)') .style('font-size', '15px'); - sr_unique_activities = [...new Set(short_range_activities)] + sr_unique_activities = [...new Set(short_range_expirations)] if (show_sr_legend) { // Long range cumulative dose line legend - line and area var legendLongCumulativeIcon = vis.append('line') @@ -196,8 +196,8 @@ function draw_plot(svg_id) { shortRangeArea[index] = d3.area(); drawShortRangeArea[index] = draw_area.append('svg:path'); - if (short_range_activities[index] == 'Breathing') drawShortRangeArea[index].attr('fill', 'red').attr('fill-opacity', '0.2'); - else if (short_range_activities[index] == 'Speaking') drawShortRangeArea[index].attr('fill', 'green').attr('fill-opacity', '0.1'); + if (short_range_expirations[index] == 'Breathing') drawShortRangeArea[index].attr('fill', 'red').attr('fill-opacity', '0.2'); + else if (short_range_expirations[index] == 'Speaking') drawShortRangeArea[index].attr('fill', 'green').attr('fill-opacity', '0.1'); else drawShortRangeArea[index].attr('fill', 'blue').attr('fill-opacity', '0.1'); }); @@ -541,6 +541,7 @@ function draw_plot(svg_id) { // Generate the alternative scenarios plot using d3 library. // 'alternative_scenarios' is a dictionary with all the alternative scenarios // 'times' is a list of times for all the scenarios +// The method is prepared to consider short range interactions if needed. function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_plot_svg_id) { // H:M format var time_format = d3.timeFormat('%H:%M'); @@ -847,7 +848,7 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ window.addEventListener("resize", e => { redraw(); if (button_full_exposure && button_full_exposure.disabled) update_alternative_concentration_plot('concentrations'); - else update_alternative_concentration_plot('concentrations_zoomed') + else update_alternative_concentration_plot('concentrations') }); } diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 9fab2cd9..84beed27 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -503,12 +503,7 @@ baseline_model = models.ExposureModel( ), evaporation_factor=0.3, ), - short_range=models.ShortRangeModel( - expirations=(), - activities=(), - presence=(), - distances=(), - ), + short_range=(), exposed=models.Population( number=10, presence=models.SpecificInterval(((8., 12.), (13., 17.))), diff --git a/cara/apps/templates/base/calculator.report.html.j2 b/cara/apps/templates/base/calculator.report.html.j2 index d35bf86a..45a8ff0d 100644 --- a/cara/apps/templates/base/calculator.report.html.j2 +++ b/cara/apps/templates/base/calculator.report.html.j2 @@ -105,7 +105,7 @@ let long_range_cumulative_doses = {{ long_range_cumulative_doses | JSONify }} let exposed_presence_intervals = {{ exposed_presence_intervals | JSONify }} let short_range_intervals = {{ short_range_intervals | JSONify }} - let short_range_activities = {{ short_range_activities | JSONify }} + let short_range_expirations = {{ short_range_expirations | JSONify }} draw_plot("concentration_plot")

diff --git a/cara/models.py b/cara/models.py index fd658861..68617e7a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1074,30 +1074,30 @@ class ConcentrationModel: @dataclass(frozen=True) class ShortRangeModel: - #: Expiration types - expirations: typing.Tuple[_ExpirationBase, ...] + #: Expiration type + expiration: _ExpirationBase - #: Activity types - activities: typing.Tuple[Activity, ...] + #: Activity type + activity: Activity - #: Short-range expirtory activities and respective interactions - presence: typing.Tuple[typing.Tuple[str, SpecificInterval], ...] + #: Short-range expiration and respective presence + presence: SpecificInterval #: Interpersonal distances distances: _VectorisedFloat - def dilution_factor(self, BR: _VectorisedFloat) -> _VectorisedFloat: + def dilution_factor(self) -> _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 + BR = self.activity.exhalation_rate/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 + u0 = BR/Am tstar = 2.0 Cr1 = 0.18 @@ -1111,15 +1111,18 @@ class ShortRangeModel: # 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) + xstar = 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) + Sxstar = 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 - ]) + distances = np.array(self.distances) + intermediate_range1 = distances < xstar + intermediate_range2 = distances >= xstar + + factors = np.empty(distances.shape, dtype=np.float64) + factors[intermediate_range1] = 2*Cr1*(distances[intermediate_range1] + x01)/D + factors[intermediate_range2] = Sxstar[intermediate_range2]*(1 + Cr2*(distances[intermediate_range2] - xstar[intermediate_range2])/Cr1/(xstar[intermediate_range2] + x01))**3 + return factors def _normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: """ @@ -1129,29 +1132,27 @@ class ShortRangeModel: short-range concentration normalized by the virus viral load. Otherwise it returns 0. """ - 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: - 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) - - # 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) - - # 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)) - + start, stop = self.presence.boundaries()[0] + # Verifies if the given time falls within a short-range interaction + if start <= time <= stop: + dilution = self.dilution_factor() + jet_origin_concentration = self.expiration.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) + + # 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.expiration.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)) return 0. - def short_range_concentration(self, concentration_model: ConcentrationModel, time: float): + def short_range_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: """ Virus short-range exposure concentration, as a function of time. """ @@ -1170,27 +1171,14 @@ class ShortRangeModel: Get the integrated short-range concentration of viruses in the air between the times start and stop, normalized by the virus viral load. """ - 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.dilution_factor(self.activities[index].exhalation_rate) + start_bound, stop_bound = self.presence.boundaries()[0] + + jet_origin_integrated = self.expiration.jet_origin_concentration() + dilution = self.dilution_factor() - 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 + 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.expiration.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) @@ -1205,8 +1193,8 @@ class ExposureModel: #: The virus concentration model which this exposure model should consider. concentration_model: ConcentrationModel - #: The short-range model which this exposure model should consider. - short_range: ShortRangeModel + #: The list of short-range models which this exposure model should consider. + short_range: typing.Tuple[ShortRangeModel, ...] #: The population of non-infected people to be used in the model. exposed: Population @@ -1248,9 +1236,11 @@ class ExposureModel: It considers the long-range concentration with the contribution of the short-range concentration. - """ - return (self.concentration_model.concentration(time) + - self.short_range.short_range_concentration(self.concentration_model, time)) + """ + concentration = self.concentration_model.concentration(time) + for interaction in self.short_range: + concentration += interaction.short_range_concentration(self.concentration_model, time) + return concentration def long_range_deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: deposited_exposure = 0. @@ -1295,8 +1285,8 @@ class ExposureModel: initial deposited exposure. """ deposited_exposure = 0. - for index, (activity, period) in enumerate(self.short_range.presence): - start, stop = tuple(period.boundaries()) + for interaction in self.short_range: + start, stop = interaction.presence.boundaries()[0] if stop < time1: continue elif start > time2: @@ -1309,10 +1299,10 @@ class ExposureModel: start_bound, stop_bound = start, time2 elif time1 <= start and stop < time2: start_bound, stop_bound = start, stop - short_range_exposure = self.short_range.normed_exposure_between_bounds(self.concentration_model, start_bound, stop_bound) + short_range_exposure = interaction.normed_exposure_between_bounds(self.concentration_model, start_bound, stop_bound) - fdep = self.short_range.expirations[index].particle.fraction_deposited(evaporation_factor=1.0) - diameter = self.short_range.expirations[index].particle.diameter + 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. if diameter is not None and not np.isscalar(diameter): diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 58700c7b..2c6268e6 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -182,41 +182,6 @@ def expiration_distribution( ) -# 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 - - expiration_BLO_factors = { 'Breathing': (1., 0., 0.), 'Speaking': (1., 1., 1.), @@ -237,4 +202,5 @@ short_range_expiration_distributions = { } +# Fit from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm short_range_distances = LogNormal(-0.269359136417347, 0.4728300188814934) diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py index 3b7b4b88..3cce3eb3 100644 --- a/cara/tests/conftest.py +++ b/cara/tests/conftest.py @@ -31,11 +31,7 @@ def baseline_concentration_model(): @pytest.fixture def baseline_sr_model(): - return models.ShortRangeModel( - presence=(), - expirations=(), - dilutions=(), - ) + return () @pytest.fixture diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index f5a98fbd..44979024 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -171,11 +171,8 @@ def conc_model(): @pytest.fixture def sr_model(): - return models.ShortRangeModel( - presence=(), - expirations=(), - dilutions=(), - ) + return () + # Expected deposited exposure were computed with a trapezoidal integration, using # a mesh of 10'000 pts per exposed presence interval. diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index bc1546da..94aaba0b 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -6,7 +6,7 @@ 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 dilution_factor, short_range_expiration_distributions +from cara.monte_carlo.data import short_range_expiration_distributions, short_range_distances, activity_distributions @pytest.fixture @@ -15,7 +15,7 @@ def concentration_model() -> mc_models.ConcentrationModel: room=models.Room(volume=75), ventilation=models.AirChange( active=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), - air_exch=30., + air_exch=10_000_000., ), infected=mc_models.InfectedPopulation( number=1, @@ -31,69 +31,72 @@ def concentration_model() -> mc_models.ConcentrationModel: @pytest.fixture -def presences() -> typing.Tuple[typing.Tuple[str, models.SpecificInterval], ...]: - return ( - ('Breathing', models.SpecificInterval((10.5, 11.0))), - ('Speaking', models.SpecificInterval((14.5, 15.0))), - ('Shouting', models.SpecificInterval((16.5, 17.5))), +def short_range_model(): + return mc_models.ShortRangeModel(expiration=short_range_expiration_distributions['Breathing'], + activity=activity_distributions['Seated'], + presence=models.SpecificInterval(present_times=((10.5, 11.0),),), + distances=short_range_distances) + + +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) + 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.short_range_concentration(concentration_model, 14.0), float) + + +@pytest.mark.parametrize( + "activity, expected_dilution", [ + ["Seated", 1016.558562890949], + ["Standing", 258.848724384952], + ["Light activity", 101.71457048643958], + ["Moderate activity", 77.93123076916304], + ["Heavy exercise", 145.70265022304739], + ] +) +def test_dilution_factor(activity, expected_dilution): + model = mc_models.ShortRangeModel(expiration="Breathing", + activity=activity_distributions[activity], + presence=models.SpecificInterval(present_times=((10.5, 11.0),),), + distances=short_range_distances).build_model(10) + assert isinstance(model.dilution_factor(), np.ndarray) + np.testing.assert_almost_equal( + model.dilution_factor().mean(), expected_dilution, decimal=10 ) -@pytest.fixture -def expirations(presences) -> typing.Tuple[mc_models.Expiration, ...]: - # Monte carlo expirations! So build model. - return tuple(short_range_expiration_distributions[activity] for (activity, presence) in presences) - - -@pytest.fixture -def dilutions(presences): - return dilution_factor(activities=[activity for (activity, presence) in presences]) - - -def test_short_range_model_ndarray(concentration_model, presences, expirations, dilutions): +@pytest.mark.parametrize( + "time, expected_short_range_concentration", [ + [8.5, 0.], + [10.5, 15.124713216706837], + [10.6, 15.193703513876928], + [11.0, 15.265132134078277], + [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 = 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.5), np.ndarray) - + model = short_range_model.build_model(250_000) + np.testing.assert_almost_equal( + np.array(model.short_range_concentration(concentration_model, time)).mean(), expected_short_range_concentration + ) @pytest.mark.parametrize( "start, stop, expected_exposure", [ - [8.5, 12.5, 5.963061627547172e-10], - [10.5, 11.0, 5.934552264225482e-10], - [10.6, 11.9, 4.709684623109963e-10], + [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, - presences, expirations, dilutions): +def test_normed_exposure_between_bounds(start, stop, expected_exposure, concentration_model, short_range_model): concentration_model = concentration_model.build_model(250_000) - model = mc_models.ShortRangeModel(presences, expirations, dilutions) - model = 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, decimal=10 - ) - - -@pytest.mark.parametrize( - "time, expected_sr_normed_concentration, expected_concentration", [ - [10.75, 1.1670056689678455e-08, 1.1731466540816526], - # [14.75, 3.6414877020308386e-06, 474.36376505412574], - # [16.75, 1.973757599365769e-05, 2850.6219623225024], - ] -) -def test_short_range_model( - time, expected_sr_normed_concentration, expected_concentration, - 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_concentration(concentration_model, time).mean(), expected_sr_normed_concentration, decimal=0 - ) - np.testing.assert_almost_equal( - model.short_range_concentration(concentration_model, time).mean(), expected_concentration, decimal=0 + model.normed_exposure_between_bounds(concentration_model, start, stop).mean(), expected_exposure ) diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index 0e4d2bb0..a398e716 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -63,11 +63,7 @@ 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( - presence=(), - expirations=(), - dilutions=(), - ) + return () @pytest.fixture diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 26911a5f..e9d8ac7c 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -39,17 +39,8 @@ TorontoTemperatures = { # references values for infection_probability and expected new cases # in the following tests, were obtained from the feature/mc branch - @pytest.fixture -def sr_model_mc() -> mc.ShortRangeModel: - return mc.ShortRangeModel( - presence=(), - expirations=(), - dilutions=(), - ) - -@pytest.fixture -def shared_office_mc(sr_model_mc): +def shared_office_mc(): """ Corresponds to the 1st line of Table 4 in https://doi.org/10.1101/2021.10.14.21264988 """ @@ -80,7 +71,7 @@ def shared_office_mc(sr_model_mc): ) return mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=3, presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))), @@ -92,7 +83,7 @@ def shared_office_mc(sr_model_mc): @pytest.fixture -def classroom_mc(sr_model_mc): +def classroom_mc(): """ Corresponds to the 2nd line of Table 4 in https://doi.org/10.1101/2021.10.14.21264988 """ @@ -123,7 +114,7 @@ def classroom_mc(sr_model_mc): ) return mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=19, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), @@ -135,7 +126,7 @@ def classroom_mc(sr_model_mc): @pytest.fixture -def ski_cabin_mc(sr_model_mc): +def ski_cabin_mc(): """ Corresponds to the 3rd line of Table 4 in https://doi.org/10.1101/2021.10.14.21264988 """ @@ -157,7 +148,7 @@ def ski_cabin_mc(sr_model_mc): ) return mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=3, presence=models.SpecificInterval(((0, 20/60),)), @@ -169,7 +160,7 @@ def ski_cabin_mc(sr_model_mc): @pytest.fixture -def skagit_chorale_mc(sr_model_mc): +def skagit_chorale_mc(): """ Corresponds to the 4th line of Table 4 in https://doi.org/10.1101/2021.10.14.21264988, assuming viral is 10**9 instead of a LogCustomKernel distribution. @@ -197,7 +188,7 @@ def skagit_chorale_mc(sr_model_mc): ) return mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=60, presence=models.SpecificInterval(((0, 2.5), )), @@ -209,7 +200,7 @@ def skagit_chorale_mc(sr_model_mc): @pytest.fixture -def bus_ride_mc(sr_model_mc): +def bus_ride_mc(): """ Corresponds to the 5th line of Table 4 in https://doi.org/10.1101/2021.10.14.21264988, assuming viral is 5*10**8 instead of a LogCustomKernel distribution. @@ -237,7 +228,7 @@ def bus_ride_mc(sr_model_mc): ) return mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=67, presence=models.SpecificInterval(((0, 1.67), )), @@ -249,7 +240,7 @@ def bus_ride_mc(sr_model_mc): @pytest.fixture -def gym_mc(sr_model_mc): +def gym_mc(): """ Gym model for testing """ @@ -272,7 +263,7 @@ def gym_mc(sr_model_mc): ) return mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=28, presence=concentration_mc.infected.presence, @@ -284,7 +275,7 @@ def gym_mc(sr_model_mc): @pytest.fixture -def waiting_room_mc(sr_model_mc): +def waiting_room_mc(): """ Waiting room model for testing """ @@ -307,7 +298,7 @@ def waiting_room_mc(sr_model_mc): ) return mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=14, presence=concentration_mc.infected.presence, @@ -355,7 +346,7 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, - expected_dose, expected_ER, sr_model_mc): + expected_dose, expected_ER): concentration_mc = mc.ConcentrationModel( room=models.Room(volume=33, humidity=0.5), ventilation=models.MultipleVentilation( @@ -385,7 +376,7 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi, ) exposure_mc = mc.ExposureModel( concentration_model=concentration_mc, - short_range=sr_model_mc, + short_range=(), exposed=mc.Population( number=1, presence=concentration_mc.infected.presence,