From e7a1c831b723209371e937a22eee44c9760c4ec1 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 20 Jul 2021 12:36:21 +0200 Subject: [PATCH 01/30] Concentration profile plot with cumulative dose #fixes164 --- cara/apps/calculator/report_generator.py | 38 ++++++++++++++++++------ 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index a899a18f..4d2b8dc0 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -105,18 +105,20 @@ def img2base64(img_data) -> str: # A src suitable for a tag such as f'. return f'data:image/png;base64,{pic_hash}' - def plot(times, concentrations, model: models.ExposureModel): fig = plt.figure() ax = fig.add_subplot(1, 1, 1) datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times] - ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean concentration') + + #Concentration as mean viral concentration (virion m$^{-3}$) + concentrations = [c * model.concentration_model.virus.quantum_infectious_dose for c in concentrations] + + ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean viral concentration') ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - ax.set_xlabel('Time of day') - ax.set_ylabel('Mean concentration ($q/m^3$)') - ax.set_title('Mean concentration of infectious quanta') + ax.set_xlabel('Time of day', fontsize=14) + ax.set_ylabel('Mean viral concentration\n(virion m$^{-3}$)', fontsize=14) + ax.set_title('Concentration profile') ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) # Plot presence of exposed person @@ -128,13 +130,31 @@ def plot(times, concentrations, model: models.ExposureModel): label="Presence of exposed person(s)" if i == 0 else "" ) + #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose + factor = 0.6 * np.mean(model.exposed.activity.inhalation_rate) * (1 - model.exposed.mask.η_inhale) + 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))] + + ax1 = ax.twinx() + ax1.plot(datetimes, qds, label='qD - Mean viral concentration', 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) + ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1]) + ax1.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) + # Place a legend outside of the axes itself. - ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - ax.set_ylim(0) + fig.legend(bbox_to_anchor=(1.05, 0.9), loc='upper left') + ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1]) + # Remove top spines + ax.spines['top'].set_visible(False) + ax1.spines['top'].set_visible(False) return fig - def minutes_to_time(minutes: int) -> str: minute_string = str(minutes % 60) minute_string = "0" * (2 - len(minute_string)) + minute_string From e07bee91ac7bd81e5e653cc13c6ff7230497601f Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 20 Jul 2021 15:34:06 +0200 Subject: [PATCH 02/30] Legend correction and order --- cara/apps/calculator/report_generator.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 4d2b8dc0..ddcf1cd9 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -123,7 +123,7 @@ def plot(times, concentrations, model: models.ExposureModel): # Plot presence of exposed person for i, (presence_start, presence_finish) in enumerate(model.exposed.presence.boundaries()): - plt.fill_between( + ax.fill_between( datetimes, concentrations, 0, where=(np.array(times) > presence_start) & (np.array(times) < presence_finish), color="#1f77b4", alpha=0.1, @@ -139,7 +139,7 @@ def plot(times, concentrations, model: models.ExposureModel): qds = [np.trapz(modified_concentrations[:i + 1], times[:i + 1]) * factor for i in range(len(times))] ax1 = ax.twinx() - ax1.plot(datetimes, qds, label='qD - Mean viral concentration', 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) @@ -147,7 +147,12 @@ def plot(times, concentrations, model: models.ExposureModel): ax1.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) # Place a legend outside of the axes itself. - fig.legend(bbox_to_anchor=(1.05, 0.9), loc='upper left') + ax_handles, ax_labels = ax.get_legend_handles_labels() + ax1_handles, ax1_labels = ax1.get_legend_handles_labels() + handles = ax_handles + ax1_handles + labels = ax_labels + ax1_labels + order = [0, 2, 1] # 0 - Mean viral concentration 1 - Presence of exposed person(s) 2 - Mean cumulative dose + fig.legend(handles = [handles[idx] for idx in order], labels = [labels[idx] for idx in order], bbox_to_anchor=(1.05, 0.9), loc='upper left') ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1]) # Remove top spines ax.spines['top'].set_visible(False) From f7b2f3c567adedd73de7bfbede628c2d6ea73444 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 21 Jul 2021 09:42:46 +0200 Subject: [PATCH 03/30] Alternative Scenarios with cumulative dose plot --- cara/apps/calculator/report_generator.py | 27 ++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index ddcf1cd9..1d8061db 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -238,6 +238,7 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: np.ndarray): fig = plt.figure() ax = fig.add_subplot(1, 1, 1) + ax1 = ax.twinx() dash_styled_scenarios = [ 'Base scenario with FFP2 masks', @@ -245,15 +246,31 @@ def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: np.ndarray) 'Base scenario with HEPA and FFP2 masks', ] - sample_dts = [datetime(1970, 1, 1) + timedelta(hours=time) for time in sample_times] + datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in sample_times] + for name, statistics in scenarios.items(): concentrations = statistics['concentrations'] - + factor = statistics['factor'] + present_indexes = statistics['present_indexes'] + + modified_concentrations = np.array(concentrations) + modified_concentrations[~present_indexes] = 0 + qds = [np.trapz(modified_concentrations[:i + 1], sample_times[:i + 1]) * factor for i in range(len(sample_times))] + if name in dash_styled_scenarios: - ax.plot(sample_dts, concentrations, label=name, linestyle='--') + ax.plot(datetimes, concentrations, label=name, linestyle='--') + ax1.plot(datetimes, qds, label='Mean cumulative dose', linestyle='dotted') else: - ax.plot(sample_dts, concentrations, label=name, linestyle='-', alpha=0.5) + ax.plot(datetimes, concentrations, label=name, linestyle='-', alpha=0.5) + ax1.plot(datetimes, qds, label='Mean cumulative dose', linestyle='dotted', alpha=0.5) + + + ax1.spines["right"].set_linestyle("--") + ax1.spines["right"].set_linestyle((0,(1,5))) + ax1.set_ylabel('Mean cumulative dose\n(virion)', fontsize=14) + ax1.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) + # Place a legend outside of the axes itself. ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') ax.spines['right'].set_visible(False) @@ -276,6 +293,8 @@ def scenario_statistics(mc_model: mc.ExposureModel, sample_times: np.ndarray): np.mean(model.concentration_model.concentration(time)) for time in sample_times ], + 'factor': 0.6 * np.mean(model.exposed.activity.inhalation_rate) * (1 - model.exposed.mask.η_inhale), + 'present_indexes': np.array([model.exposed.person_present(t) for t in sample_times]), } From 85045f02dbaf198b3908618d7e72e7f7b49d2690 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 21 Jul 2021 11:01:41 +0200 Subject: [PATCH 04/30] Alternative Scenarios with cumulative dose fix #227 --- cara/apps/calculator/report_generator.py | 41 ++++++++++--------- .../templates/base/calculator.report.html.j2 | 4 +- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 1d8061db..22aaf932 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -249,37 +249,41 @@ def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: np.ndarray) datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in sample_times] for name, statistics in scenarios.items(): + model = statistics['model'] concentrations = statistics['concentrations'] - factor = statistics['factor'] - present_indexes = statistics['present_indexes'] - + + #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose + factor = 0.6 * np.mean(model.exposed.activity.inhalation_rate) * (1 - model.exposed.mask.η_inhale) + present_indexes = np.array([model.exposed.person_present(t) for t in sample_times]) + modified_concentrations = np.array(concentrations) modified_concentrations[~present_indexes] = 0 qds = [np.trapz(modified_concentrations[:i + 1], sample_times[:i + 1]) * factor for i in range(len(sample_times))] + # Plot concentrations and cumulative dose if name in dash_styled_scenarios: ax.plot(datetimes, concentrations, label=name, linestyle='--') - ax1.plot(datetimes, qds, label='Mean cumulative dose', linestyle='dotted') + ax1.plot(datetimes, qds, label=f'Mean cumulative dose:\n{name}', linestyle='dotted') else: ax.plot(datetimes, concentrations, label=name, linestyle='-', alpha=0.5) - ax1.plot(datetimes, qds, label='Mean cumulative dose', linestyle='dotted', alpha=0.5) - + ax1.plot(datetimes, qds, label=f'Mean cumulative dose:\n{name}', linestyle='dotted', alpha=0.5) + + # Place a legend outside of the axes itself. + fig.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left') + + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + ax.set_xlabel('Time of day', fontsize=14) + ax.set_ylabel('Mean viral concentration\n(virion m$^{-3}$)', fontsize=14) + ax.set_title('Concentration profile') + ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) + ax1.spines['top'].set_visible(False) ax1.spines["right"].set_linestyle("--") ax1.spines["right"].set_linestyle((0,(1,5))) ax1.set_ylabel('Mean cumulative dose\n(virion)', fontsize=14) ax1.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) - - # Place a legend outside of the axes itself. - ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) - - ax.set_xlabel('Time of day') - ax.set_ylabel('Mean concentration ($q/m^3$)') - ax.set_title('Mean concentration of infectious quanta') return fig @@ -287,14 +291,13 @@ def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: np.ndarray) def scenario_statistics(mc_model: mc.ExposureModel, sample_times: np.ndarray): model = mc_model.build_model(size=_DEFAULT_MC_SAMPLE_SIZE) return { + 'model': model, 'probability_of_infection': np.mean(model.infection_probability()), 'expected_new_cases': np.mean(model.expected_new_cases()), 'concentrations': [ - np.mean(model.concentration_model.concentration(time)) + np.mean(model.concentration_model.concentration(time)) * model.concentration_model.virus.quantum_infectious_dose for time in sample_times ], - 'factor': 0.6 * np.mean(model.exposed.activity.inhalation_rate) * (1 - model.exposed.mask.η_inhale), - 'present_indexes': np.array([model.exposed.person_present(t) for t in sample_times]), } diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 64a44cfa..98655d8c 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -219,7 +219,7 @@

Alternative scenarios:

-

+ {% block report_scenarios_summary_table %} @@ -242,7 +242,7 @@ {% endblock report_scenarios_summary_table %} -

+

Notes for alternative scenarios:
From bc38f5ec7a4950ccba4b16c175014e1fef4b9457 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 23 Jul 2021 09:15:25 +0200 Subject: [PATCH 05/30] cumulated_exposure method in models.py --- cara/apps/calculator/report_generator.py | 4 ++-- cara/models.py | 9 ++++++--- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 22aaf932..217e0902 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -131,7 +131,7 @@ def plot(times, concentrations, model: models.ExposureModel): ) #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose - factor = 0.6 * np.mean(model.exposed.activity.inhalation_rate) * (1 - model.exposed.mask.η_inhale) + 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 @@ -139,7 +139,7 @@ def plot(times, concentrations, model: models.ExposureModel): qds = [np.trapz(modified_concentrations[:i + 1], times[:i + 1]) * factor for i in range(len(times))] ax1 = ax.twinx() - ax1.plot(datetimes, qds, label='Mean cumulative dose', color='#1f77b4', linestyle='dotted') + ax1.plot(datetimes, cumulated_exposure, 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 585fb824..b762b451 100644 --- a/cara/models.py +++ b/cara/models.py @@ -876,15 +876,18 @@ class ExposureModel: return exposure * self.repeats - def infection_probability(self) -> _VectorisedFloat: + def cumulated_exposure(self) -> _VectorisedFloat: exposure = self.quanta_exposure() - - inf_aero = ( + + return ( self.exposed.activity.inhalation_rate * (1 - self.exposed.mask.inhale_efficiency()) * exposure * self.fraction_deposited ) + def infection_probability(self) -> _VectorisedFloat: + inf_aero = self.cumulated_exposure() + # Probability of infection. return (1 - np.exp(-inf_aero)) * 100 From 460c5cf520af94ff2143e83f3d14d2f3c53a4998 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 23 Jul 2021 10:57:52 +0200 Subject: [PATCH 06/30] 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 From 680db38cce62da76ed5f4513676f04c8fab09f75 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 23 Jul 2021 12:15:50 +0200 Subject: [PATCH 07/30] pipelined fixed for returning value models.py --- cara/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index b115d2f4..93785d1d 100644 --- a/cara/models.py +++ b/cara/models.py @@ -100,7 +100,7 @@ class Interval: return True return False - def generate_truncated_interval(self, time_stop: float) -> tuple: + def generate_truncated_interval(self, time_stop: float) -> "Interval": truncated_boundaries = [] for start, end in self.boundaries(): if start < time_stop <= end: @@ -112,7 +112,7 @@ class Interval: else: truncated_boundaries.append((start, end)) - return SpecificInterval(tuple(truncated_boundaries)) + return SpecificInterval(present_times = tuple(truncated_boundaries)) @dataclass(frozen=True) class SpecificInterval(Interval): From 56e8714633c25c5bd25d1ca964beee7ba2ab34ba Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 23 Jul 2021 14:31:46 +0200 Subject: [PATCH 08/30] Concentration plot for alternative scenarios fix --- cara/apps/calculator/report_generator.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index f2933735..3a8d3197 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -248,12 +248,7 @@ def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: np.ndarray) concentrations = statistics['concentrations'] #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose - factor = 0.6 * np.mean(model.exposed.activity.inhalation_rate) * (1 - model.exposed.mask.η_inhale) - present_indexes = np.array([model.exposed.person_present(t) for t in sample_times]) - - modified_concentrations = np.array(concentrations) - modified_concentrations[~present_indexes] = 0 - qds = [np.trapz(modified_concentrations[:i + 1], sample_times[:i + 1]) * factor for i in range(len(sample_times))] + qds = [np.mean(dataclass_utils.nested_replace(model, {'exposed.presence': model.exposed.presence.generate_truncated_interval(t)}).cumulated_exposure()) for t in sample_times] # Plot concentrations and cumulative dose if name in dash_styled_scenarios: From 7bab71653262e394e98307b8586e28a236f03359 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 23 Jul 2021 15:11:19 +0200 Subject: [PATCH 09/30] added @cached to integrated_concentrations method --- cara/models.py | 1 + 1 file changed, 1 insertion(+) diff --git a/cara/models.py b/cara/models.py index 93785d1d..27bd965a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -839,6 +839,7 @@ class ConcentrationModel: fac = np.exp(-IVRR * delta_time) return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac + @cached() def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: """ Get the integrated concentration dose between the times start and stop. From 5983fca92ebcba13745db4f8f004d62d7fcf79ea Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 26 Jul 2021 09:55:56 +0200 Subject: [PATCH 10/30] test exposure model fix --- cara/apps/calculator/report_generator.py | 7 ++-- cara/models.py | 47 ++++++++++++++---------- cara/tests/models/test_exposure_model.py | 29 +++++++++------ cara/tests/models/test_interval.py | 22 ----------- 4 files changed, 48 insertions(+), 57 deletions(-) delete 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 3a8d3197..b4db8d28 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -10,6 +10,7 @@ import zlib import loky import jinja2 import matplotlib +from numpy.lib.function_base import quantile matplotlib.use('agg') import matplotlib.pyplot as plt import numpy as np @@ -131,8 +132,8 @@ def plot(times, concentrations, model: models.ExposureModel): ) #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose - qds = [np.mean(dataclass_utils.nested_replace(model, {'exposed.presence': model.exposed.presence.generate_truncated_interval(t)}).cumulated_exposure()) for t in times] - + qds = [np.mean(model.cumulated_exposure_vs_time(t)) for t in times] + ax1 = ax.twinx() ax1.plot(datetimes, qds, label='Mean cumulative dose', color='#1f77b4', linestyle='dotted') ax1.spines["right"].set_linestyle("--") @@ -248,7 +249,7 @@ def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: np.ndarray) concentrations = statistics['concentrations'] #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose - qds = [np.mean(dataclass_utils.nested_replace(model, {'exposed.presence': model.exposed.presence.generate_truncated_interval(t)}).cumulated_exposure()) for t in sample_times] + qds = [np.mean(model.cumulated_exposure_vs_time(t)) for t in sample_times] # Plot concentrations and cumulative dose if name in dash_styled_scenarios: diff --git a/cara/models.py b/cara/models.py index 27bd965a..924d91cc 100644 --- a/cara/models.py +++ b/cara/models.py @@ -100,20 +100,6 @@ class Interval: return True return False - def generate_truncated_interval(self, time_stop: float) -> "Interval": - 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(present_times = tuple(truncated_boundaries)) - @dataclass(frozen=True) class SpecificInterval(Interval): #: A sequence of times (start, stop), in hours, that the infected person @@ -881,17 +867,32 @@ class ExposureModel: #: The fraction of viruses actually deposited in the respiratory tract fraction_deposited: _VectorisedFloat = 0.6 - def quanta_exposure(self) -> _VectorisedFloat: - """The number of virus quanta per meter^3.""" + def quanta_exposure_vs_time(self, time: float) -> _VectorisedFloat: + """The number of virus quanta per meter^3 integrated until time.""" exposure = 0.0 for start, stop in self.exposed.presence.boundaries(): - exposure += self.concentration_model.integrated_concentration(start, stop) - + if start > time: + break + elif time <= stop: + stop = time + exposure += self.concentration_model.integrated_concentration(start, stop) + break + else: + exposure += self.concentration_model.integrated_concentration(start, stop) + return exposure * self.repeats - def cumulated_exposure(self) -> _VectorisedFloat: - exposure = self.quanta_exposure() + def quanta_exposure(self) -> _VectorisedFloat: + """The number of virus quanta per meter^3 for the full simulation time.""" + if self.exposed.presence.transition_times(): + return self.quanta_exposure_vs_time(max(self.exposed.presence.transition_times())) + else: + return 0 + + def cumulated_exposure_vs_time(self, time: float) -> _VectorisedFloat: + + exposure = self.quanta_exposure_vs_time(time) return ( self.exposed.activity.inhalation_rate * @@ -899,6 +900,12 @@ class ExposureModel: exposure * self.fraction_deposited ) + def cumulated_exposure(self) -> _VectorisedFloat: + if self.exposed.presence.transition_times(): + return self.cumulated_exposure_vs_time(max(self.exposed.presence.transition_times())) + else: + return 0 + def infection_probability(self) -> _VectorisedFloat: inf_aero = self.cumulated_exposure() diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 298a9c93..e3e35e87 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -4,26 +4,31 @@ import typing import numpy as np import numpy.testing import pytest +from dataclasses import dataclass from cara import models from cara.models import ExposureModel +@dataclass(frozen=True) class KnownConcentrations(models.ConcentrationModel): """ A ConcentrationModel which is based on pre-known quanta concentrations and which therefore doesn't need other components. Useful for testing. """ - def __init__(self, concentration_function: typing.Callable) -> None: - self._func = concentration_function + #def __init__(self, concentration_function: typing.Callable) -> None: + # self._func = concentration_function + + + concentration_function: typing.Callable def infectious_virus_removal_rate(self, time: float) -> models._VectorisedFloat: # very large decay constant -> same as constant concentration return 1.e50 def _concentration_limit(self, time: float) -> models._VectorisedFloat: - return self._func(time) + return self.concentration_function(time) def state_change_times(self): return [0, 24] @@ -32,7 +37,7 @@ class KnownConcentrations(models.ConcentrationModel): return 24 def concentration(self, time: float) -> models._VectorisedFloat: # noqa - return self._func(time) + return self.concentration_function(time) halftime = models.PeriodicInterval(120, 60) @@ -57,19 +62,19 @@ populations = [ @pytest.mark.parametrize( "population, cm, f_dep, expected_exposure, expected_cumulated_exposure, expected_probability",[ - [populations[1], KnownConcentrations(lambda t: 1.2), 1., + [populations[1], KnownConcentrations(None, None, None, lambda t: 1.2), 1., 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., + [populations[2], KnownConcentrations(None, None, None, lambda t: 1.2), 1., 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., + [populations[0], KnownConcentrations(None, None, None,lambda t: np.array([1.2, 2.4])), 1., 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., + [populations[1], KnownConcentrations(None, None, None,lambda t: np.array([1.2, 2.4])), 1., 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.]), + [populations[0], KnownConcentrations(None, None, None,lambda t: 2.4), np.array([0.5, 1.]), 28.8, np.array([4.104, 8.208]), np.array([98.3493482895, 99.9727534893])], ]) def test_exposure_model_ndarray(population, cm, f_dep, @@ -95,7 +100,7 @@ def test_exposure_model_ndarray(population, cm, f_dep, @pytest.mark.parametrize("population", populations) def test_exposure_model_ndarray_and_float_mix(population): - cm = KnownConcentrations(lambda t: 0 if np.floor(t) % 2 else np.array([1.2, 1.2])) + cm = KnownConcentrations(None, None, None, lambda t: 0 if np.floor(t) % 2 else np.array([1.2, 1.2])) model = ExposureModel(cm, population) expected_exposure = np.array([14.4, 14.4]) @@ -109,8 +114,8 @@ def test_exposure_model_ndarray_and_float_mix(population): @pytest.mark.parametrize("population", populations) def test_exposure_model_compare_scalar_vector(population): - cm_scalar = KnownConcentrations(lambda t: 1.2) - cm_array = KnownConcentrations(lambda t: np.array([1.2, 1.2])) + cm_scalar = KnownConcentrations(None, None, None,lambda t: 1.2) + cm_array = KnownConcentrations(None, None, None, lambda t: np.array([1.2, 1.2])) model_scalar = ExposureModel(cm_scalar, population) model_array = ExposureModel(cm_array, population) expected_exposure = 14.4 diff --git a/cara/tests/models/test_interval.py b/cara/tests/models/test_interval.py deleted file mode 100644 index 5aa73603..00000000 --- a/cara/tests/models/test_interval.py +++ /dev/null @@ -1,22 +0,0 @@ -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 From 1803d4e919ee826304deecae9cd670282f5ca549 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 26 Jul 2021 14:16:51 +0200 Subject: [PATCH 11/30] mypy errors - added some dummy variables to test_exposure_model.py --- cara/tests/models/test_exposure_model.py | 32 +++++++++++++++--------- 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index e3e35e87..1b06c23f 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -17,10 +17,6 @@ class KnownConcentrations(models.ConcentrationModel): which therefore doesn't need other components. Useful for testing. """ - #def __init__(self, concentration_function: typing.Callable) -> None: - # self._func = concentration_function - - concentration_function: typing.Callable def infectious_virus_removal_rate(self, time: float) -> models._VectorisedFloat: @@ -58,23 +54,35 @@ populations = [ models.Activity(np.array([0.51,0.57]), 0.57), ), ] +dummyRoom = models.Room(50, 0.5) +dummyVentilation = models._VentilationBase() +dummyInfPopulation = models.InfectedPopulation( + number=1, + presence=halftime, + mask=models.Mask.types['Type I'], + activity=models.Activity.types['Standing'], + virus=models.Virus.types['SARS_CoV_2_B117'], + expiration=models.Expiration.types['Talking'] +) +def known_concentrations(func): + return KnownConcentrations(dummyRoom, dummyVentilation, dummyInfPopulation, func) @pytest.mark.parametrize( "population, cm, f_dep, expected_exposure, expected_cumulated_exposure, expected_probability",[ - [populations[1], KnownConcentrations(None, None, None, lambda t: 1.2), 1., + [populations[1], known_concentrations(lambda t: 1.2), 1., 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(None, None, None, lambda t: 1.2), 1., + [populations[2], known_concentrations(lambda t: 1.2), 1., 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(None, None, None,lambda t: np.array([1.2, 2.4])), 1., + [populations[0], known_concentrations(lambda t: np.array([1.2, 2.4])), 1., 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(None, None, None,lambda t: np.array([1.2, 2.4])), 1., + [populations[1], known_concentrations(lambda t: np.array([1.2, 2.4])), 1., 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(None, None, None,lambda t: 2.4), np.array([0.5, 1.]), + [populations[0], known_concentrations(lambda t: 2.4), np.array([0.5, 1.]), 28.8, np.array([4.104, 8.208]), np.array([98.3493482895, 99.9727534893])], ]) def test_exposure_model_ndarray(population, cm, f_dep, @@ -100,7 +108,7 @@ def test_exposure_model_ndarray(population, cm, f_dep, @pytest.mark.parametrize("population", populations) def test_exposure_model_ndarray_and_float_mix(population): - cm = KnownConcentrations(None, None, None, lambda t: 0 if np.floor(t) % 2 else np.array([1.2, 1.2])) + cm = known_concentrations(lambda t: 0 if np.floor(t) % 2 else np.array([1.2, 1.2])) model = ExposureModel(cm, population) expected_exposure = np.array([14.4, 14.4]) @@ -114,8 +122,8 @@ def test_exposure_model_ndarray_and_float_mix(population): @pytest.mark.parametrize("population", populations) def test_exposure_model_compare_scalar_vector(population): - cm_scalar = KnownConcentrations(None, None, None,lambda t: 1.2) - cm_array = KnownConcentrations(None, None, None, lambda t: np.array([1.2, 1.2])) + cm_scalar = known_concentrations(lambda t: 1.2) + cm_array = known_concentrations(lambda t: np.array([1.2, 1.2])) model_scalar = ExposureModel(cm_scalar, population) model_array = ExposureModel(cm_array, population) expected_exposure = 14.4 From f329d76e2224c721ae46016bf232458b5e3f24d7 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 26 Jul 2021 14:20:16 +0200 Subject: [PATCH 12/30] Performed some calculations missed on the code --- cara/tests/models/test_exposure_model.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 1b06c23f..acf85d3f 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -71,16 +71,16 @@ def known_concentrations(func): @pytest.mark.parametrize( "population, cm, f_dep, expected_exposure, expected_cumulated_exposure, expected_probability",[ [populations[1], known_concentrations(lambda t: 1.2), 1., - np.array([14.4, 14.4]), np.array([3.44736/0.6, 3.20112/0.6]), np.array([99.6803184113, 99.5181053773])], + np.array([14.4, 14.4]), np.array([5.7456, 5.3352]), np.array([99.6803184113, 99.5181053773])], [populations[2], known_concentrations(lambda t: 1.2), 1., - np.array([14.4, 14.4]), np.array([2.2032/0.6, 2.4624/0.6]), np.array([97.4574432074, 98.3493482895])], + np.array([14.4, 14.4]), np.array([3.672, 4.104]), np.array([97.4574432074, 98.3493482895])], [populations[0], known_concentrations(lambda t: np.array([1.2, 2.4])), 1., - np.array([14.4, 28.8]), np.array([2.4624/0.6, 4.9248/0.6]), np.array([98.3493482895, 99.9727534893])], + np.array([14.4, 28.8]), np.array([4.104, 8.208]), np.array([98.3493482895, 99.9727534893])], [populations[1], known_concentrations(lambda t: np.array([1.2, 2.4])), 1., - np.array([14.4, 28.8]), np.array([3.44736/0.6, 6.40224/0.6]), np.array([99.6803184113, 99.9976777757])], + np.array([14.4, 28.8]), np.array([5.7456, 10.6704]), np.array([99.6803184113, 99.9976777757])], [populations[0], known_concentrations(lambda t: 2.4), np.array([0.5, 1.]), 28.8, np.array([4.104, 8.208]), np.array([98.3493482895, 99.9727534893])], From 072c47d70a01cd2fd6b55ce02acdbcf2c6ab13fb Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 17 Aug 2021 15:31:32 +0200 Subject: [PATCH 13/30] updated variable names and type --- cara/apps/calculator/report_generator.py | 7 +++++-- cara/apps/calculator/static/js/report.js | 10 +++++----- .../templates/base/calculator.report.html.j2 | 4 ++-- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 66d66910..5124f94a 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -44,12 +44,15 @@ def calculate_report_data(model: models.ExposureModel): er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean() exposed_occupants = model.exposed.number expected_new_cases = np.array(model.expected_new_cases()).mean() - cumulative_dose = np.array([model.cumulated_exposure_vs_time(t) for t in times]).mean() + cumulative_doses = [ + np.array(model.cumulated_exposure_vs_time(float(time))).mean() + for time in times + ] return { "times": list(times), "exposed_presence_intervals": [list(interval) for interval in model.exposed.presence.boundaries()], - "cumulative_dose": cumulative_dose, + "cumulative_doses": cumulative_doses, "concentrations": concentrations, "highest_const": highest_const, "prob_inf": prob, diff --git a/cara/apps/calculator/static/js/report.js b/cara/apps/calculator/static/js/report.js index 58ab40ae..69ef2ee2 100644 --- a/cara/apps/calculator/static/js/report.js +++ b/cara/apps/calculator/static/js/report.js @@ -1,5 +1,5 @@ /* Generate the concentration plot using d3 library. */ -function draw_concentration_plot(svg_id, times, concentrations, cumulative_dose, exposed_presence_intervals) { +function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses, exposed_presence_intervals) { var visBoundingBox = d3.select(svg_id) .node() .getBoundingClientRect(); @@ -7,7 +7,7 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_dose, var time_format = d3.timeFormat('%H:%M'); var data = [] - times.map((time, index) => data.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': concentrations[index], 'cumulative_dose': cumulative_dose[index] })) + times.map((time, index) => data.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': concentrations[index], 'cumulative_doses': cumulative_doses[index] })) var vis = d3.select(svg_id), width = visBoundingBox.width - 400, @@ -20,7 +20,7 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_dose, bisecHour = d3.bisector((d) => { return d.hour; }).left, yRange = d3.scaleLinear().range([height - margins.bottom, margins.top]).domain([0., Math.max(...concentrations)]), - yCumulatedRange = d3.scaleLinear().range([height - margins.bottom, margins.top]).domain([0., Math.max(...cumulative_dose)]), + yCumulatedRange = d3.scaleLinear().range([height - margins.bottom, margins.top]).domain([0., Math.max(...cumulative_doses)]), xAxis = d3.axisBottom(xRange).tickFormat(d => time_format(d)), yAxis = d3.axisLeft(yRange), @@ -49,9 +49,9 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_dose, // Line representing the cumulative concentration. var lineCumulativeFunc = d3.line() - .defined(d => !isNaN(d.cumulative_dose)) + .defined(d => !isNaN(d.cumulative_doses)) .x(d => xTimeRange(d.time)) - .y(d => yCumulatedRange(d.cumulative_dose)) + .y(d => yCumulatedRange(d.cumulative_doses)) .curve(d3.curveBasis); vis.append('svg:path') diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index be127ece..b4cf7fb7 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -224,9 +224,9 @@

Alternative scenarios:

From ab181b534543b2ca293daf2ec45a2f1c58335769 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 20 Aug 2021 16:44:27 +0200 Subject: [PATCH 14/30] removed unused import --- cara/tests/models/test_exposure_model.py | 1 - 1 file changed, 1 deletion(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 3c32ccd3..67238bc2 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -1,4 +1,3 @@ -import decimal import typing import numpy as np From 67b60442199f52ba811297f2dcbcfa78b92fc217 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 20 Aug 2021 16:45:36 +0200 Subject: [PATCH 15/30] changed virus to virions --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 319c9782..d48263d0 100644 --- a/cara/models.py +++ b/cara/models.py @@ -894,7 +894,7 @@ class ExposureModel: return exposure * self.repeats def exposure(self) -> _VectorisedFloat: - """The number of virus per meter^3 for the full simulation time.""" + """The number of virions per meter^3 for the full simulation time.""" if self.exposed.presence.transition_times(): return self.exposure_vs_time(max(self.exposed.presence.transition_times())) else: From cb334cbafa16b0832a047edd2493c9bdf712f5c8 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 2 Sep 2021 16:19:32 +0200 Subject: [PATCH 16/30] variable renaming for the exposure and inhaled_exposure --- cara/apps/calculator/report_generator.py | 4 ++-- cara/models.py | 16 ++++++++-------- cara/tests/models/test_exposure_model.py | 10 +++++----- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 93e3f695..8d40cd41 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -114,7 +114,7 @@ def calculate_report_data(model: models.ExposureModel): exposed_occupants = model.exposed.number expected_new_cases = np.array(model.expected_new_cases()).mean() cumulative_doses = [ - np.array(model.cumulated_exposure_vs_time(float(time))).mean() + np.array(model.inhaled_exposure_between_bounds(float(time))).mean() for time in times ] @@ -276,7 +276,7 @@ def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: typing.List concentrations = statistics['concentrations'] #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose - qds = [np.mean(model.cumulated_exposure_vs_time(t)) for t in sample_times] + qds = [np.mean(model.inhaled_exposure_between_bounds(t)) for t in sample_times] # Plot concentrations and cumulative dose if name in dash_styled_scenarios: diff --git a/cara/models.py b/cara/models.py index d48263d0..fc165ccb 100644 --- a/cara/models.py +++ b/cara/models.py @@ -877,8 +877,8 @@ class ExposureModel: #: The fraction of viruses actually deposited in the respiratory tract fraction_deposited: _VectorisedFloat = 0.6 - def exposure_vs_time(self, time: float) -> _VectorisedFloat: - """The number of virus per meter^3 integrated until time.""" + def exposure_between_bounds(self, time: float) -> _VectorisedFloat: + """The cumulative number of virions per meter^3 from model start to the given time.""" exposure = 0.0 for start, stop in self.exposed.presence.boundaries(): @@ -896,13 +896,13 @@ class ExposureModel: def exposure(self) -> _VectorisedFloat: """The number of virions per meter^3 for the full simulation time.""" if self.exposed.presence.transition_times(): - return self.exposure_vs_time(max(self.exposed.presence.transition_times())) + return self.exposure_between_bounds(max(self.exposed.presence.transition_times())) else: return 0 - def cumulated_exposure_vs_time(self, time: float) -> _VectorisedFloat: + def inhaled_exposure_between_bounds(self, time: float) -> _VectorisedFloat: - exposure = self.exposure_vs_time(time) + exposure = self.exposure_between_bounds(time) return ( self.exposed.activity.inhalation_rate * @@ -910,14 +910,14 @@ class ExposureModel: exposure * self.fraction_deposited ) - def cumulated_exposure(self) -> _VectorisedFloat: + def inhaled_exposure(self) -> _VectorisedFloat: if self.exposed.presence.transition_times(): - return self.cumulated_exposure_vs_time(max(self.exposed.presence.transition_times())) + return self.inhaled_exposure_between_bounds(max(self.exposed.presence.transition_times())) else: return 0 def infection_probability(self) -> _VectorisedFloat: - inf_aero = self.cumulated_exposure() + inf_aero = self.inhaled_exposure() # Probability of infection. return (1 - np.exp(-(inf_aero/self.concentration_model.virus.infectious_dose))) * 100 diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 67238bc2..df530f7e 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -70,7 +70,7 @@ def known_concentrations(func): @pytest.mark.parametrize( - "population, cm, f_dep, expected_exposure, expected_cumulated_exposure, expected_probability",[ + "population, cm, f_dep, expected_exposure, expected_inhaled_exposure, expected_probability",[ [populations[1], known_concentrations(lambda t: 36.), 1., np.array([432, 432]), np.array([172.368, 160.056]), np.array([99.6803184113, 99.5181053773])], @@ -87,7 +87,7 @@ def known_concentrations(func): 864, np.array([123.12, 246.24]), np.array([98.3493482895, 99.9727534893])], ]) def test_exposure_model_ndarray(population, cm, f_dep, - expected_exposure, expected_cumulated_exposure, expected_probability): + expected_exposure, expected_inhaled_exposure, expected_probability): model = ExposureModel(cm, population, fraction_deposited=f_dep) np.testing.assert_almost_equal( model.exposure(), expected_exposure @@ -96,15 +96,15 @@ def test_exposure_model_ndarray(population, cm, f_dep, model.infection_probability(), expected_probability, decimal=10 ) np.testing.assert_almost_equal( - model.cumulated_exposure(), expected_cumulated_exposure, decimal=10 + model.inhaled_exposure(), expected_inhaled_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 isinstance(model.inhaled_exposure(), np.ndarray) assert model.infection_probability().shape == (2,) assert model.expected_new_cases().shape == (2,) - assert model.cumulated_exposure().shape == (2,) + assert model.inhaled_exposure().shape == (2,) @pytest.mark.parametrize("population", populations) From 485518f3f409c99fa23d30bb80f4011875c5972f Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 3 Sep 2021 09:35:57 +0200 Subject: [PATCH 17/30] new time with breaks included --- cara/apps/calculator/report_generator.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 8d40cd41..5966d6f7 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -10,7 +10,7 @@ import zlib import loky import jinja2 import matplotlib -from numpy.lib.function_base import quantile +from numpy.lib.function_base import append, quantile matplotlib.use('agg') import matplotlib.pyplot as plt import numpy as np @@ -113,6 +113,24 @@ def calculate_report_data(model: models.ExposureModel): er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean() exposed_occupants = model.exposed.number expected_new_cases = np.array(model.expected_new_cases()).mean() + print(len(times)) + times2 = np.array([]) + for start, stop in model.exposed.presence.boundaries(): + times2 = np.concatenate((times2, np.linspace(start, stop, int(len(times)/len(model.exposed.presence.boundaries()))))) + + + # for start, stop in model.exposed.presence.boundaries(): + # if start > time: + # break + # elif time <= stop: + # stop = time + # times2.append(time) + # break + # else: + # times2.append(time) + + + print(len(times2)) cumulative_doses = [ np.array(model.inhaled_exposure_between_bounds(float(time))).mean() for time in times From d3e56f69b2ab314c29dfc0598eb3d72bb932595f Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 6 Sep 2021 10:00:53 +0200 Subject: [PATCH 18/30] get boundaries from time function --- cara/apps/calculator/report_generator.py | 42 +++++++++++++++--------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 5966d6f7..1789d865 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -101,6 +101,20 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L return nice_times +def get_boundaries_from_time(model: models.ExposureModel, time: float): + x = [] + for start, stop in model.exposed.presence.boundaries(): + if start > time: + break + elif time <= stop: + stop = time + x.append([start, stop]) + break + else: + x.append([start, stop]) + return x + + def calculate_report_data(model: models.ExposureModel): times = interesting_times(model) @@ -113,24 +127,20 @@ def calculate_report_data(model: models.ExposureModel): er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean() exposed_occupants = model.exposed.number expected_new_cases = np.array(model.expected_new_cases()).mean() - print(len(times)) - times2 = np.array([]) - for start, stop in model.exposed.presence.boundaries(): - times2 = np.concatenate((times2, np.linspace(start, stop, int(len(times)/len(model.exposed.presence.boundaries()))))) + # present_indexes = [model.exposed.person_present(t) for t in times] + # filtered_times = np.array(times)[present_indexes] - # for start, stop in model.exposed.presence.boundaries(): - # if start > time: - # break - # elif time <= stop: - # stop = time - # times2.append(time) - # break - # else: - # times2.append(time) - - - print(len(times2)) + b = [ + np.array(model.cenas(list(get_boundaries_from_time(model, time)))).mean() + for time in times + ] + print(b) + + # cumulative_doses = [ + # np.array(model.inhaled_exposure_between_bounds(list(get_boundaries_from_time(model, time)))).mean() + # for time in times + # ] cumulative_doses = [ np.array(model.inhaled_exposure_between_bounds(float(time))).mean() for time in times From 0b83c6491e2391ead72ded8bff4aea4346e118e9 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 7 Sep 2021 15:30:54 +0200 Subject: [PATCH 19/30] Changed methods names and variables. New logic to calculate exposure between bounds --- cara/apps/calculator/report_generator.py | 105 +++++++++-------------- cara/models.py | 15 ++-- 2 files changed, 51 insertions(+), 69 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 1789d865..10bd67ee 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -1,3 +1,10 @@ +from ... import dataclass_utils +from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE +from ... import monte_carlo as mc +from cara import models +import qrcode +import numpy as np +import matplotlib.pyplot as plt import concurrent.futures import base64 import dataclasses @@ -12,14 +19,6 @@ import jinja2 import matplotlib from numpy.lib.function_base import append, quantile matplotlib.use('agg') -import matplotlib.pyplot as plt -import numpy as np -import qrcode - -from cara import models -from ... import monte_carlo as mc -from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE -from ... import dataclass_utils def model_start_end(model: models.ExposureModel): @@ -97,24 +96,11 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L # Expand the times list to ensure that we have a maximum gap size between # the key times. - nice_times = fill_big_gaps(times, gap_size=(max(times) - min(times)) / approx_n_pts) + nice_times = fill_big_gaps(times, gap_size=( + max(times) - min(times)) / approx_n_pts) return nice_times -def get_boundaries_from_time(model: models.ExposureModel, time: float): - x = [] - for start, stop in model.exposed.presence.boundaries(): - if start > time: - break - elif time <= stop: - stop = time - x.append([start, stop]) - break - else: - x.append([start, stop]) - return x - - def calculate_report_data(model: models.ExposureModel): times = interesting_times(model) @@ -124,23 +110,10 @@ def calculate_report_data(model: models.ExposureModel): ] highest_const = max(concentrations) prob = np.array(model.infection_probability()).mean() - er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean() + er = np.array( + model.concentration_model.infected.emission_rate_when_present()).mean() exposed_occupants = model.exposed.number expected_new_cases = np.array(model.expected_new_cases()).mean() - - # present_indexes = [model.exposed.person_present(t) for t in times] - # filtered_times = np.array(times)[present_indexes] - - b = [ - np.array(model.cenas(list(get_boundaries_from_time(model, time)))).mean() - for time in times - ] - print(b) - - # cumulative_doses = [ - # np.array(model.inhaled_exposure_between_bounds(list(get_boundaries_from_time(model, time)))).mean() - # for time in times - # ] cumulative_doses = [ np.array(model.inhaled_exposure_between_bounds(float(time))).mean() for time in times @@ -199,7 +172,8 @@ def _img2bytes(figure): def _figure2bytes(figure): # Draw the image img_data = io.BytesIO() - figure.savefig(img_data, format='png', bbox_inches="tight", transparent=True) + figure.savefig(img_data, format='png', + bbox_inches="tight", transparent=True) return img_data @@ -250,14 +224,18 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp scenarios = {} # Two special option cases - HEPA and/or FFP2 masks. - FFP2_being_worn = bool(form.mask_wearing_option == 'mask_on' and form.mask_type == 'FFP2') + FFP2_being_worn = bool(form.mask_wearing_option == + 'mask_on' and form.mask_type == 'FFP2') if FFP2_being_worn and form.hepa_option: - FFP2andHEPAalternative = dataclass_utils.replace(form, mask_type='Type I') + FFP2andHEPAalternative = dataclass_utils.replace( + form, mask_type='Type I') scenarios['Base scenario with HEPA filter and Type I masks'] = FFP2andHEPAalternative.build_mc_model() if not FFP2_being_worn and form.hepa_option: - noHEPAalternative = dataclass_utils.replace(form, mask_type = 'FFP2') - noHEPAalternative = dataclass_utils.replace(noHEPAalternative, mask_wearing_option = 'mask_on') - noHEPAalternative = dataclass_utils.replace(noHEPAalternative, hepa_option=False) + noHEPAalternative = dataclass_utils.replace(form, mask_type='FFP2') + noHEPAalternative = dataclass_utils.replace( + noHEPAalternative, mask_wearing_option='mask_on') + noHEPAalternative = dataclass_utils.replace( + noHEPAalternative, hepa_option=False) scenarios['Base scenario without HEPA filter, with FFP2 masks'] = noHEPAalternative.build_mc_model() # The remaining scenarios are based on Type I masks (possibly not worn) @@ -267,7 +245,8 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp form = dataclass_utils.replace(form, hepa_option=False) with_mask = dataclass_utils.replace(form, mask_wearing_option='mask_on') - without_mask = dataclass_utils.replace(form, mask_wearing_option='mask_off') + without_mask = dataclass_utils.replace( + form, mask_wearing_option='mask_off') if form.ventilation_type == 'mechanical_ventilation': #scenarios['Mechanical ventilation with Type I masks'] = with_mask.build_mc_model() @@ -278,8 +257,10 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp scenarios['Windows open without masks'] = without_mask.build_mc_model() # No matter the ventilation scheme, we include scenarios which don't have any ventilation. - with_mask_no_vent = dataclass_utils.replace(with_mask, ventilation_type='no_ventilation') - without_mask_or_vent = dataclass_utils.replace(without_mask, ventilation_type='no_ventilation') + with_mask_no_vent = dataclass_utils.replace( + with_mask, ventilation_type='no_ventilation') + without_mask_or_vent = dataclass_utils.replace( + without_mask, ventilation_type='no_ventilation') scenarios['No ventilation with Type I masks'] = with_mask_no_vent.build_mc_model() scenarios['Neither ventilation nor masks'] = without_mask_or_vent.build_mc_model() @@ -297,27 +278,23 @@ def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: typing.List 'Base scenario with HEPA and FFP2 masks', ] - datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in sample_times] + datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) + for time in sample_times] for name, statistics in scenarios.items(): model = statistics['model'] concentrations = statistics['concentrations'] - #See CERN-OPEN-2021-004, p. 15, eq. 16. - Cumulative Dose - qds = [np.mean(model.inhaled_exposure_between_bounds(t)) for t in sample_times] - - # Plot concentrations and cumulative dose + # Plot concentrations if name in dash_styled_scenarios: ax.plot(datetimes, concentrations, label=name, linestyle='--') - ax1.plot(datetimes, qds, label=f'Mean cumulative dose:\n{name}', linestyle='dotted') else: - ax.plot(datetimes, concentrations, label=name, linestyle='-', alpha=0.5) - ax1.plot(datetimes, qds, label=f'Mean cumulative dose:\n{name}', linestyle='dotted', alpha=0.5) - - + ax.plot(datetimes, concentrations, + label=name, linestyle='-', alpha=0.5) + # Place a legend outside of the axes itself. fig.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left') - + ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.set_xlabel('Time of day', fontsize=14) @@ -379,7 +356,8 @@ class ReportGenerator: executor_factory: typing.Callable[[], concurrent.futures.Executor], ) -> str: model = form.build_model() - context = self.prepare_context(base_url, model, form, executor_factory=executor_factory) + context = self.prepare_context( + base_url, model, form, executor_factory=executor_factory) return self.render(context) def prepare_context( @@ -405,14 +383,15 @@ class ReportGenerator: context['alternative_scenarios'] = comparison_report( alternative_scenarios, scenario_sample_times, executor_factory=executor_factory, ) - context['qr_code'] = generate_qr_code(base_url, self.calculator_prefix, form) + context['qr_code'] = generate_qr_code( + base_url, self.calculator_prefix, form) context['calculator_prefix'] = self.calculator_prefix context['scale_warning'] = { - 'level': 'yellow-2', + 'level': 'yellow-2', 'incidence_rate': 'lower than 25 new cases per 100 000 inhabitants', - 'onsite_access': 'of about 8000', + 'onsite_access': 'of about 8000', 'threshold': '' - } + } return context def _template_environment(self) -> jinja2.Environment: diff --git a/cara/models.py b/cara/models.py index fc165ccb..772c0324 100644 --- a/cara/models.py +++ b/cara/models.py @@ -878,20 +878,23 @@ class ExposureModel: fraction_deposited: _VectorisedFloat = 0.6 def exposure_between_bounds(self, time: float) -> _VectorisedFloat: - """The cumulative number of virions per meter^3 from model start to the given time.""" - exposure = 0.0 - + """The number of virions per meter^3 from model start to the given time.""" + boundaries = [] for start, stop in self.exposed.presence.boundaries(): if start > time: break elif time <= stop: stop = time - exposure += self.concentration_model.integrated_concentration(start, stop) + boundaries.append([start, stop]) break else: - exposure += self.concentration_model.integrated_concentration(start, stop) + boundaries.append([start, stop]) + + exposure = 0.0 + for start, stop in boundaries: + exposure += self.concentration_model.integrated_concentration(start, stop) - return exposure * self.repeats + return exposure def exposure(self) -> _VectorisedFloat: """The number of virions per meter^3 for the full simulation time.""" From c27f567f3ab5aa936b528eb3460cfb0de5d48ae7 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 7 Sep 2021 17:00:51 +0200 Subject: [PATCH 20/30] merge with master --- cara/apps/calculator/model_generator.py | 7 +- cara/apps/calculator/report_generator.py | 111 ++------ cara/apps/calculator/static/js/report.js | 257 ++++++++++++++---- .../templates/base/calculator.report.html.j2 | 20 +- cara/apps/templates/common_text.md.j2 | 3 + cara/models.py | 1 - cara/monte_carlo/data.py | 12 +- cara/monte_carlo/sampleable.py | 12 + 8 files changed, 264 insertions(+), 159 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index a6a41fab..06c29229 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -11,7 +11,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 +from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions LOG = logging.getLogger(__name__) @@ -359,7 +359,10 @@ class FormData: def mask(self) -> models.Mask: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as # the "No mask"-mask - mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] + if self.mask_wearing_option == 'mask_on': + mask = mask_distributions[self.mask_type] + else: + mask = models.Mask.types['No mask'] return mask def infected_population(self) -> mc.InfectedPopulation: diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 10bd67ee..23ad9fb1 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -1,10 +1,3 @@ -from ... import dataclass_utils -from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE -from ... import monte_carlo as mc -from cara import models -import qrcode -import numpy as np -import matplotlib.pyplot as plt import concurrent.futures import base64 import dataclasses @@ -16,9 +9,14 @@ import zlib import loky import jinja2 -import matplotlib -from numpy.lib.function_base import append, quantile -matplotlib.use('agg') +import numpy as np +import qrcode +import json + +from cara import models +from ... import monte_carlo as mc +from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE +from ... import dataclass_utils def model_start_end(model: models.ExposureModel): @@ -96,8 +94,7 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L # Expand the times list to ensure that we have a maximum gap size between # the key times. - nice_times = fill_big_gaps(times, gap_size=( - max(times) - min(times)) / approx_n_pts) + nice_times = fill_big_gaps(times, gap_size=(max(times) - min(times)) / approx_n_pts) return nice_times @@ -110,8 +107,7 @@ def calculate_report_data(model: models.ExposureModel): ] highest_const = max(concentrations) prob = np.array(model.infection_probability()).mean() - er = np.array( - model.concentration_model.infected.emission_rate_when_present()).mean() + er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean() exposed_occupants = model.exposed.number expected_new_cases = np.array(model.expected_new_cases()).mean() cumulative_doses = [ @@ -169,21 +165,12 @@ def _img2bytes(figure): return img_data -def _figure2bytes(figure): - # Draw the image - img_data = io.BytesIO() - figure.savefig(img_data, format='png', - bbox_inches="tight", transparent=True) - return img_data - - def img2base64(img_data) -> str: - plt.close() img_data.seek(0) pic_hash = base64.b64encode(img_data.read()).decode('ascii') # A src suitable for a tag such as f'. return f'data:image/png;base64,{pic_hash}' - + def minutes_to_time(minutes: int) -> str: minute_string = str(minutes % 60) @@ -224,18 +211,14 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp scenarios = {} # Two special option cases - HEPA and/or FFP2 masks. - FFP2_being_worn = bool(form.mask_wearing_option == - 'mask_on' and form.mask_type == 'FFP2') + FFP2_being_worn = bool(form.mask_wearing_option == 'mask_on' and form.mask_type == 'FFP2') if FFP2_being_worn and form.hepa_option: - FFP2andHEPAalternative = dataclass_utils.replace( - form, mask_type='Type I') + FFP2andHEPAalternative = dataclass_utils.replace(form, mask_type='Type I') scenarios['Base scenario with HEPA filter and Type I masks'] = FFP2andHEPAalternative.build_mc_model() if not FFP2_being_worn and form.hepa_option: - noHEPAalternative = dataclass_utils.replace(form, mask_type='FFP2') - noHEPAalternative = dataclass_utils.replace( - noHEPAalternative, mask_wearing_option='mask_on') - noHEPAalternative = dataclass_utils.replace( - noHEPAalternative, hepa_option=False) + noHEPAalternative = dataclass_utils.replace(form, mask_type = 'FFP2') + noHEPAalternative = dataclass_utils.replace(noHEPAalternative, mask_wearing_option = 'mask_on') + noHEPAalternative = dataclass_utils.replace(noHEPAalternative, hepa_option=False) scenarios['Base scenario without HEPA filter, with FFP2 masks'] = noHEPAalternative.build_mc_model() # The remaining scenarios are based on Type I masks (possibly not worn) @@ -245,8 +228,7 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp form = dataclass_utils.replace(form, hepa_option=False) with_mask = dataclass_utils.replace(form, mask_wearing_option='mask_on') - without_mask = dataclass_utils.replace( - form, mask_wearing_option='mask_off') + without_mask = dataclass_utils.replace(form, mask_wearing_option='mask_off') if form.ventilation_type == 'mechanical_ventilation': #scenarios['Mechanical ventilation with Type I masks'] = with_mask.build_mc_model() @@ -257,62 +239,17 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp scenarios['Windows open without masks'] = without_mask.build_mc_model() # No matter the ventilation scheme, we include scenarios which don't have any ventilation. - with_mask_no_vent = dataclass_utils.replace( - with_mask, ventilation_type='no_ventilation') - without_mask_or_vent = dataclass_utils.replace( - without_mask, ventilation_type='no_ventilation') + with_mask_no_vent = dataclass_utils.replace(with_mask, ventilation_type='no_ventilation') + without_mask_or_vent = dataclass_utils.replace(without_mask, ventilation_type='no_ventilation') scenarios['No ventilation with Type I masks'] = with_mask_no_vent.build_mc_model() scenarios['Neither ventilation nor masks'] = without_mask_or_vent.build_mc_model() return scenarios -def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: typing.List[float]): - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1) - ax1 = ax.twinx() - - dash_styled_scenarios = [ - 'Base scenario with FFP2 masks', - 'Base scenario with HEPA filter', - 'Base scenario with HEPA and FFP2 masks', - ] - - datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) - for time in sample_times] - - for name, statistics in scenarios.items(): - model = statistics['model'] - concentrations = statistics['concentrations'] - - # Plot concentrations - if name in dash_styled_scenarios: - ax.plot(datetimes, concentrations, label=name, linestyle='--') - else: - ax.plot(datetimes, concentrations, - label=name, linestyle='-', alpha=0.5) - - # Place a legend outside of the axes itself. - fig.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left') - - ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - ax.set_xlabel('Time of day', fontsize=14) - ax.set_ylabel('Mean viral concentration\n(virion m$^{-3}$)', fontsize=14) - ax.set_title('Concentration profile') - ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) - - ax.set_xlabel('Time of day') - ax.set_ylabel('Mean concentration ($virions/m^{3}$)') - ax.set_title('Mean concentration of virions') - - return fig - - -def scenario_statistics(mc_model: mc.ExposureModel, sample_times: typing.List[float]): +def scenario_statistics(mc_model: mc.ExposureModel, sample_times: np.ndarray): model = mc_model.build_model(size=_DEFAULT_MC_SAMPLE_SIZE) return { - 'model': model, 'probability_of_infection': np.mean(model.infection_probability()), 'expected_new_cases': np.mean(model.expected_new_cases()), 'concentrations': [ @@ -339,7 +276,6 @@ def comparison_report( for (name, model), model_stats in zip(scenarios.items(), results): statistics[name] = model_stats return { - 'plot': img2base64(_figure2bytes(comparison_plot(statistics, sample_times))), 'stats': statistics, } @@ -356,8 +292,7 @@ class ReportGenerator: executor_factory: typing.Callable[[], concurrent.futures.Executor], ) -> str: model = form.build_model() - context = self.prepare_context( - base_url, model, form, executor_factory=executor_factory) + context = self.prepare_context(base_url, model, form, executor_factory=executor_factory) return self.render(context) def prepare_context( @@ -383,8 +318,7 @@ class ReportGenerator: context['alternative_scenarios'] = comparison_report( alternative_scenarios, scenario_sample_times, executor_factory=executor_factory, ) - context['qr_code'] = generate_qr_code( - base_url, self.calculator_prefix, form) + context['qr_code'] = generate_qr_code(base_url, self.calculator_prefix, form) context['calculator_prefix'] = self.calculator_prefix context['scale_warning'] = { 'level': 'yellow-2', @@ -404,6 +338,7 @@ class ReportGenerator: env.filters['minutes_to_time'] = minutes_to_time env.filters['float_format'] = "{0:.2f}".format env.filters['int_format'] = "{:0.0f}".format + env.filters['JSONify'] = json.dumps return env def render(self, context: dict) -> str: diff --git a/cara/apps/calculator/static/js/report.js b/cara/apps/calculator/static/js/report.js index 00e85892..6266a09b 100644 --- a/cara/apps/calculator/static/js/report.js +++ b/cara/apps/calculator/static/js/report.js @@ -27,70 +27,18 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses yCumulatedAxis = d3.axisRight(yCumulatedRange); // Plot tittle. - vis.append('svg:foreignObject') - .attr("background-color", "transparent") - .attr('width', width) - .attr('height', margins.top) - .style('text-align', 'center') - .html('Mean concentration of virions'); + plot_title(vis, width, margins.top, 'Mean concentration of virions'); // Line representing the mean concentration. - var lineFunc = d3.line() - .defined(d => !isNaN(d.concentration)) - .x(d => xTimeRange(d.time)) - .y(d => yRange(d.concentration)) - .curve(d3.curveBasis); - - vis.append('svg:path') - .attr('d', lineFunc(data)) - .attr('stroke', '#1f77b4') - .attr('stroke-width', 2) - .attr('fill', 'none'); - + plot_scenario_data(vis, data, xTimeRange, yRange, '#1f77b4'); // Line representing the cumulative concentration. - var lineCumulativeFunc = d3.line() - .defined(d => !isNaN(d.cumulative_doses)) - .x(d => xTimeRange(d.time)) - .y(d => yCumulatedRange(d.cumulative_doses)) - .curve(d3.curveBasis); + plot_cumulative_data(vis, data, xTimeRange, yCumulatedRange, '#1f77b4'); - vis.append('svg:path') - .attr('d', lineCumulativeFunc(data)) - .attr('stroke', '#1f77b4') - .attr('stroke-width', 2) - .style("stroke-dasharray", "5 5") - .attr('fill', 'none'); + // X axis. + plot_x_axis(vis, height, width, margins, xAxis, 'Time of day'); - // X axis declaration. - vis.append('svg:g') - .attr('class', 'x axis') - .attr('transform', 'translate(0,' + (height - margins.bottom) + ')') - .call(xAxis); - - // X axis label. - vis.append('text') - .attr('class', 'x label') - .attr('fill', 'black') - .attr('text-anchor', 'middle') - .attr('x', (width + margins.right) / 2) - .attr('y', height * 0.97) - .text('Time of day') - - // Y concentration axis declaration. - vis.append('svg:g') - .attr('class', 'y axis') - .attr('transform', 'translate(' + margins.left + ',0)') - .call(yAxis); - - // Y concentration axis label. - vis.append('svg:text') - .attr('class', 'y label') - .attr('fill', 'black') - .attr('transform', 'rotate(-90, 0,' + height + ')') - .attr('text-anchor', 'middle') - .attr('x', (height + margins.bottom) / 2) - .attr('y', (height + margins.left) * 0.92) - .text('Mean concentration (virions/m³)'); + // Y axis + plot_y_axis(vis, height, width, margins, yAxis, 'Mean concentration (virions/m³)') // Y cumulative concentration axis declaration. vis.append('svg:g') @@ -230,4 +178,195 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses focus.select('#tooltip-time').text('x = ' + time_format(d.hour)); focus.select('#tooltip-concentration').text('y = ' + d.concentration.toFixed(2)); } +} + +// 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 +function draw_alternative_scenarios_plot(svg_id, width, height, alternative_scenarios, times) { + // H:M format + var time_format = d3.timeFormat('%H:%M'); + // D3 array of ten categorical colors represented as RGB hexadecimal strings. + var colors = d3.schemeAccent; + + // Variable for the highest concentration for all the scenarios + var highest_concentration = 0. + + var data_for_scenarios = {} + for (scenario in alternative_scenarios) { + scenario_concentrations = alternative_scenarios[scenario].concentrations + + highest_concentration = Math.max(highest_concentration, Math.max(...scenario_concentrations)) + + var data = [] + times.map((time, index) => data.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': scenario_concentrations[index] })) + + // Add data into lines dictionary + data_for_scenarios[scenario] = data + } + + // We need one scenario to get the time range + var first_scenario = Object.values(data_for_scenarios)[0] + + var vis = d3.select(svg_id), + width = width, + height = height, + margins = { top: 30, right: 20, bottom: 50, left: 50 }, + + // H:M time format for x axis. + xRange = d3.scaleTime().range([margins.left, width - margins.right]).domain([first_scenario[0].hour, first_scenario[first_scenario.length - 1].hour]), + xTimeRange = d3.scaleLinear().range([margins.left, width - margins.right]).domain([times[0], times[times.length - 1]]), + + yRange = d3.scaleLinear().range([height - margins.bottom, margins.top]).domain([0., highest_concentration]), + + xAxis = d3.axisBottom(xRange).tickFormat(d => time_format(d)), + yAxis = d3.axisLeft(yRange); + + // Plot title. + plot_title(vis, width, margins.top, 'Mean concentration of virions'); + + // Line representing the mean concentration for each scenario. + for (const [scenario_name, data] of Object.entries(data_for_scenarios)) { + var scenario_index = Object.keys(data_for_scenarios).indexOf(scenario_name) + + // Line representing the mean concentration. + plot_scenario_data(vis, data, xTimeRange, yRange, colors[scenario_index]) + + // Legend for the plot elements - lines. + var size = 20 * (scenario_index + 1) + vis.append('rect') + .attr('x', width + 20) + .attr('y', margins.top + size) + .attr('width', 20) + .attr('height', 3) + .style('fill', colors[scenario_index]); + + vis.append('text') + .attr('x', width + 3 * 20) + .attr('y', margins.top + size) + .text(scenario_name) + .style('font-size', '15px') + .attr('alignment-baseline', 'central'); + + } + + // X axis. + plot_x_axis(vis, height, width, margins, xAxis, "Time of day"); + + // Y axis declaration. + vis.append('svg:g') + .attr('class', 'y axis') + .attr('transform', 'translate(' + margins.left + ',0)') + .call(yAxis); + + // Y axis label. + vis.append('svg:text') + .attr('class', 'y label') + .attr('fill', 'black') + .attr('transform', 'rotate(-90, 0,' + height + ')') + .attr('text-anchor', 'middle') + .attr('x', (height + margins.bottom) / 2) + .attr('y', (height + margins.left) * 0.92) + .text('Mean concentration (virions/m³)'); + + // Legend bounding box. + vis.append('rect') + .attr('width', 275) + .attr('height', 25 * (Object.keys(data_for_scenarios).length)) + .attr('x', width * 1.005) + .attr('y', margins.top + 5) + .attr('stroke', 'lightgrey') + .attr('stroke-width', '2') + .attr('rx', '5px') + .attr('ry', '5px') + .attr('stroke-linejoin', 'round') + .attr('fill', 'none'); +} + + +// Functions used to build the plots' components + +function plot_title(vis, width, margin_top, title) { + vis.append('svg:foreignObject') + .attr('width', width) + .attr('height', margin_top) + .attr('fill', 'none') + .append('xhtml:div') + .style('text-align', 'center') + .html(title); + + return vis; +} + +function plot_x_axis(vis, height, width, margins, xAxis, label) { + // X axis declaration + vis.append('svg:g') + .attr('class', 'x axis') + .attr('transform', 'translate(0,' + (height - margins.bottom) + ')') + .call(xAxis); + + // X axis label. + vis.append('text') + .attr('class', 'x label') + .attr('fill', 'black') + .attr('text-anchor', 'middle') + .attr('x', (width + margins.right) / 2) + .attr('y', height * 0.97) + .text(label); + + return vis; +} + +function plot_y_axis(vis, height, width, margins, yAxis, label) { + // Y axis declaration. + vis.append('svg:g') + .attr('class', 'y axis') + .attr('transform', 'translate(' + margins.left + ',0)') + .call(yAxis); + + // Y axis label. + vis.append('svg:text') + .attr('class', 'y label') + .attr('fill', 'black') + .attr('transform', 'rotate(-90, 0,' + height + ')') + .attr('text-anchor', 'middle') + .attr('x', (height + margins.bottom) / 2) + .attr('y', (height + margins.left) * 0.92) + .text(label); + + return vis; + +} + +function plot_scenario_data(vis, data, xTimeRange, yRange, line_color) { + var lineFunc = d3.line() + .defined(d => !isNaN(d.concentration)) + .x(d => xTimeRange(d.time)) + .y(d => yRange(d.concentration)) + .curve(d3.curveBasis); + + vis.append('svg:path') + .attr('d', lineFunc(data)) + .attr("stroke", line_color) + .attr('stroke-width', 2) + .attr('fill', 'none'); + + return vis; +} + +function plot_cumulative_data(vis, data, xTimeRange, yCumulativeRange, line_color) { + var lineCumulativeFunc = d3.line() + .defined(d => !isNaN(d.cumulative_doses)) + .x(d => xTimeRange(d.time)) + .y(d => yCumulativeRange(d.cumulative_doses)) + .curve(d3.curveBasis); + + vis.append('svg:path') + .attr('d', lineCumulativeFunc(data)) + .attr('stroke', line_color) + .attr('stroke-width', 2) + .style("stroke-dasharray", "5 5") + .attr('fill', 'none'); + + return vis; } \ No newline at end of file diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 371788c0..108a9a9c 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -86,12 +86,12 @@ {% endblock report_summary_footnote %}

* The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004.

- +

@@ -109,8 +109,12 @@
- - + + {% block report_scenarios_summary_table %} @@ -437,4 +441,4 @@ - \ No newline at end of file + diff --git a/cara/apps/templates/common_text.md.j2 b/cara/apps/templates/common_text.md.j2 index 63512dee..ac8829d7 100644 --- a/cara/apps/templates/common_text.md.j2 +++ b/cara/apps/templates/common_text.md.j2 @@ -71,3 +71,6 @@ We wish to thank CERN’s HSE Unit, Beams Department, Experimental Physics Depar [54] Leung, N.H.L et al. Respiratory virus shedding in exhaled breath and efficacy of face masks. Nat Med (2020). 10.1038/s41591-020-0843-2.
[55] Asadi, S., Cappa, C.D., Barreda, S. et al. Efficacy of masks and face coverings in controlling outward aerosol particle emission from expiratory activities. Sci Rep 10, 15665 (2020). https://doi.org/10.1038/s41598-020-72798-7.
[56] Endo A, Abbott S et al. Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China [version 3; peer review: 2 approved]. Wellcome Open Res 2020, 5:67. doi:10.12688/wellcomeopenres.15842.3.
+[57] Jin Pan, Charbel Harb, Weinan Leng & Linsey C. Marr (2021) Inward and outward effectiveness of cloth masks, a surgical mask, and a face shield, Aerosol Science and Technology, 55:6, 718-733, doi: 10.1080/02786826.2021.1890687.
+[58] C. Makison Booth, M. Clayton, B. Crook, J.M. Gawn, Effectiveness of surgical masks against influenza bioaerosols, Journal of Hospital Infection, Volume 84, Issue 1, 2013, Pages 22-26, https://doi.org/10.1016/j.jhin.2013.02.007.
+[59] Riediker, M., Monn, C. (2021). Simulation of SARS-CoV-2 Aerosol Emissions in the Infected Population and Resulting Airborne Exposures in Different Indoor Scenarios. Aerosol Air Qual. Res. 21, 200531. https://doi.org/10.4209/aaqr.2020.08.0531.
\ No newline at end of file diff --git a/cara/models.py b/cara/models.py index 772c0324..a86457eb 100644 --- a/cara/models.py +++ b/cara/models.py @@ -904,7 +904,6 @@ class ExposureModel: return 0 def inhaled_exposure_between_bounds(self, time: float) -> _VectorisedFloat: - exposure = self.exposure_between_bounds(time) return ( diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 613b830c..4415efbb 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -1,7 +1,7 @@ import numpy as np import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel +from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel, Uniform # From CERN-OPEN-2021-04 and refererences therein @@ -58,3 +58,13 @@ virus_distributions = { infectious_dose=60/1.6, ), } + + +# From: +# https://doi.org/10.1080/02786826.2021.1890687 +# https://doi.org/10.1016/j.jhin.2013.02.007 +# https://doi.org/10.4209/aaqr.2020.08.0531 +mask_distributions = { + 'Type I': mc.Mask(Uniform(0.25, 0.80)), + 'FFP2': mc.Mask(Uniform(0.83, 0.91)), +} \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 4333c93f..27907e49 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -28,6 +28,18 @@ class Normal(SampleableDistribution): return np.random.normal(self.mean, self.standard_deviation, size=size) +class Uniform(SampleableDistribution): + """ + Defines a continuous uniform distribution + """ + def __init__(self, low: float, high: float): + self.low = low + self.high = high + + def generate_samples(self, size: int) -> float_array_size_n: + return np.random.uniform(self.low, self.high, size=size) + + class LogNormal(SampleableDistribution): """ Defines a lognormal distribution (i.e. Gaussian distribution vs. the From a689fe8e743096d05bb1f10f061a3d01ad16d884 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 7 Sep 2021 17:06:41 +0200 Subject: [PATCH 21/30] bug with width fixed --- cara/apps/calculator/templates/base/calculator.report.html.j2 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 108a9a9c..6e7056a6 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -86,7 +86,7 @@ {% endblock report_summary_footnote %}

* The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004.

- +