From 460c5cf520af94ff2143e83f3d14d2f3c53a4998 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 23 Jul 2021 10:57:52 +0200 Subject: [PATCH] Cumulative dose calculation code in models fixes #164 --- cara/apps/calculator/report_generator.py | 9 ++------- cara/models.py | 13 +++++++++++++ cara/tests/models/test_exposure_model.py | 22 ++++++++++++++-------- cara/tests/models/test_interval.py | 22 ++++++++++++++++++++++ 4 files changed, 51 insertions(+), 15 deletions(-) create mode 100644 cara/tests/models/test_interval.py diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 217e0902..f2933735 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -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) diff --git a/cara/models.py b/cara/models.py index b762b451..b115d2f4 100644 --- a/cara/models.py +++ b/cara/models.py @@ -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): diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 1c6c6d3d..298a9c93 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -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) diff --git a/cara/tests/models/test_interval.py b/cara/tests/models/test_interval.py new file mode 100644 index 00000000..5aa73603 --- /dev/null +++ b/cara/tests/models/test_interval.py @@ -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