pingu_98/cara
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Replaced scipy method with analytical method and updated short range tests

Browse source
This commit is contained in:
Luis Aleixo 2022-03-16 16:58:30 +00:00
parent e40120b1bd
commit 7e1dfa2684
9 changed files with 70 additions and 37 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
cara/apps/calculator/model_generator.py
View file

@ -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 ()

6
cara/apps/calculator/report_generator.py
View file

@ -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

1
cara/apps/expert.py
View file

@ -504,7 +504,6 @@ baseline_model = models.ExposureModel(
evaporation_factor=0.3,
),
short_range=models.ShortRangeModel(
activities=(),
presence=(),
expirations=(),
dilutions=(),

48
cara/models.py
View file

@ -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.

1
cara/tests/conftest.py
View file

@ -32,7 +32,6 @@ def baseline_concentration_model():
@pytest.fixture
def baseline_sr_model():
return models.ShortRangeModel(
activities=(),
presence=(),
expirations=(),
dilutions=(),

1
cara/tests/models/test_exposure_model.py
View file

@ -172,7 +172,6 @@ def conc_model():
@pytest.fixture
def sr_model():
return models.ShortRangeModel(
activities=(),
presence=(),
expirations=(),
dilutions=(),

43
cara/tests/models/test_short_range_model.py
View file

@ -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

1
cara/tests/test_monte_carlo.py
View file

@ -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=(),

1
cara/tests/test_monte_carlo_full_models.py
View file

@ -43,7 +43,6 @@ TorontoTemperatures = {
@pytest.fixture
def sr_model_mc() -> mc.ShortRangeModel:
return mc.ShortRangeModel(
activities=(),
presence=(),
expirations=(),
dilutions=(),