Cumulative dose calculation code in models fixes #164

This commit is contained in:
Luis Aleixo 2021-07-23 10:57:52 +02:00
parent bc38f5ec7a
commit 460c5cf520
4 changed files with 51 additions and 15 deletions

View file

@ -131,15 +131,10 @@ def plot(times, concentrations, model: models.ExposureModel):
)
#See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose
cumulated_exposure = model.cumulated_exposure()
present_indexes = np.array([model.exposed.person_present(t) for t in times])
modified_concentrations = np.array(concentrations)
modified_concentrations[~present_indexes] = 0
qds = [np.trapz(modified_concentrations[:i + 1], times[:i + 1]) * factor for i in range(len(times))]
qds = [np.mean(dataclass_utils.nested_replace(model, {'exposed.presence': model.exposed.presence.generate_truncated_interval(t)}).cumulated_exposure()) for t in times]
ax1 = ax.twinx()
ax1.plot(datetimes, cumulated_exposure, label='Mean cumulative dose', color='#1f77b4', linestyle='dotted')
ax1.plot(datetimes, qds, label='Mean cumulative dose', color='#1f77b4', linestyle='dotted')
ax1.spines["right"].set_linestyle("--")
ax1.spines["right"].set_linestyle((0,(1,5)))
ax1.set_ylabel('Mean cumulative dose\n(virion)', fontsize=14)

View file

@ -100,6 +100,19 @@ class Interval:
return True
return False
def generate_truncated_interval(self, time_stop: float) -> tuple:
truncated_boundaries = []
for start, end in self.boundaries():
if start < time_stop <= end:
end = time_stop
truncated_boundaries.append((start, end))
break
elif time_stop <= start:
break
else:
truncated_boundaries.append((start, end))
return SpecificInterval(tuple(truncated_boundaries))
@dataclass(frozen=True)
class SpecificInterval(Interval):

View file

@ -1,3 +1,4 @@
import decimal
import typing
import numpy as np
@ -55,24 +56,24 @@ populations = [
@pytest.mark.parametrize(
"population, cm, f_dep, expected_exposure, expected_probability",[
"population, cm, f_dep, expected_exposure, expected_cumulated_exposure, expected_probability",[
[populations[1], KnownConcentrations(lambda t: 1.2), 1.,
np.array([14.4, 14.4]), np.array([99.6803184113, 99.5181053773])],
np.array([14.4, 14.4]), np.array([3.44736/0.6, 3.20112/0.6]), np.array([99.6803184113, 99.5181053773])],
[populations[2], KnownConcentrations(lambda t: 1.2), 1.,
np.array([14.4, 14.4]), np.array([97.4574432074, 98.3493482895])],
np.array([14.4, 14.4]), np.array([2.2032/0.6, 2.4624/0.6]), np.array([97.4574432074, 98.3493482895])],
[populations[0], KnownConcentrations(lambda t: np.array([1.2, 2.4])), 1.,
np.array([14.4, 28.8]), np.array([98.3493482895, 99.9727534893])],
np.array([14.4, 28.8]), np.array([2.4624/0.6, 4.9248/0.6]), np.array([98.3493482895, 99.9727534893])],
[populations[1], KnownConcentrations(lambda t: np.array([1.2, 2.4])), 1.,
np.array([14.4, 28.8]), np.array([99.6803184113, 99.9976777757])],
np.array([14.4, 28.8]), np.array([3.44736/0.6, 6.40224/0.6]), np.array([99.6803184113, 99.9976777757])],
[populations[0], KnownConcentrations(lambda t: 2.4), np.array([0.5, 1.]),
28.8, np.array([98.3493482895, 99.9727534893])],
28.8, np.array([4.104, 8.208]), np.array([98.3493482895, 99.9727534893])],
])
def test_exposure_model_ndarray(population, cm, f_dep,
expected_exposure, expected_probability):
def test_exposure_model_ndarray(population, cm, f_dep,
expected_exposure, expected_cumulated_exposure, expected_probability):
model = ExposureModel(cm, population, fraction_deposited = f_dep)
np.testing.assert_almost_equal(
model.quanta_exposure(), expected_exposure
@ -80,11 +81,16 @@ def test_exposure_model_ndarray(population, cm, f_dep,
np.testing.assert_almost_equal(
model.infection_probability(), expected_probability, decimal=10
)
np.testing.assert_almost_equal(
model.cumulated_exposure(), expected_cumulated_exposure, decimal=10
)
assert isinstance(model.infection_probability(), np.ndarray)
assert isinstance(model.expected_new_cases(), np.ndarray)
assert isinstance(model.cumulated_exposure(), np.ndarray)
assert model.infection_probability().shape == (2,)
assert model.expected_new_cases().shape == (2,)
assert model.cumulated_exposure().shape == (2,)
@pytest.mark.parametrize("population", populations)

View file

@ -0,0 +1,22 @@
import numpy as np
import pytest
from cara import models
@pytest.mark.parametrize(
"stop_time, expected_boundaries",
[
[1.05, ((0, 1),)],
[1.8, ((0, 1), (1.1, 1.8))],
[2., ((0, 1), (1.1, 1.999))],
[3., ((0, 1), (1.1, 1.999), (2, 3))],
[-1, ()],
[4, ((0, 1), (1.1, 1.999), (2, 3))],
],
)
def test_interval_truncation(stop_time, expected_boundaries):
interesting_times = models.SpecificInterval(
([0, 1], [1.1, 1.999], [2, 3]), )
assert interesting_times.generate_truncated_interval(
stop_time).boundaries() == expected_boundaries