diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 859f8e49..28f99b85 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, monte_carlo -import qrcode -import numpy as np -import matplotlib.pyplot as plt import concurrent.futures import base64 import dataclasses @@ -18,8 +11,14 @@ import loky import jinja2 import matplotlib matplotlib.use('agg') +import matplotlib.pyplot as plt +import numpy as np +import qrcode -from cara.monte_carlo.data import activity_distributions, virus_distributions +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): @@ -39,8 +38,7 @@ def calculate_report_data(model: models.ExposureModel): for time in times] 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() @@ -111,64 +109,28 @@ def img2base64(img_data) -> str: 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') + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) - points = 600 - viral_loads = np.linspace(2, 12, points) + ax.set_xlabel('Time of day') + ax.set_ylabel('Mean concentration ($q/m^3$)') + ax.set_title('Mean concentration of virions') + ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) - er_means = [] - er_medians = [] - lower_percentiles = [] - upper_percentiles = [] - - for vl in viral_loads: - exposure_mc = mc.ExposureModel( - concentration_model=mc.ConcentrationModel( - room=models.Room(volume=100, humidity=0.5), - ventilation=models.AirChange( - active=models.SpecificInterval(((0, 24),)), - air_exch=0.25, - ), - infected=mc.InfectedPopulation( - number=1, - virus=models.Virus( - viral_load_in_sputum = 10**vl, - infectious_dose = 50., - ), - presence=mc.SpecificInterval(((0, 2),)), - mask=models.Mask.types["No mask"], - activity=activity_distributions['Seated'], - expiration=models.MultipleExpiration( - expirations=(models.Expiration.types['Talking'], - models.Expiration.types['Breathing']), - weights=(0.3, 0.7)), - ), - ), - exposed=mc.Population( - number=14, - presence=mc.SpecificInterval(((0, 2),)), - activity=models.Activity.types['Seated'], - mask=models.Mask.types["No mask"], - ), + # Plot presence of exposed person + for i, (presence_start, presence_finish) in enumerate(model.exposed.presence.boundaries()): + plt.fill_between( + datetimes, concentrations, 0, + where=(np.array(times) > presence_start) & (np.array(times) < presence_finish), + color="#1f77b4", alpha=0.1, + label="Presence of exposed person(s)" if i == 0 else "" ) - exposure_model = exposure_mc.build_model(size=_DEFAULT_MC_SAMPLE_SIZE) - print(exposure_model) - - emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present() - er_means.append(np.mean(emission_rate)) - er_medians.append(np.median(emission_rate)) - lower_percentiles.append(np.quantile(emission_rate, 0.01)) - upper_percentiles.append(np.quantile(emission_rate, 0.99)) - - ax.plot(viral_loads, er_means) - ax.fill_between(viral_loads, lower_percentiles, - upper_percentiles, alpha=0.2) - ax.set_yscale('log') - plt.title('Emission rate vs Viral Load') - plt.ylabel('ER (Virions/h)', fontsize=14) - plt.xticks(ticks=[i for i in range(2, 13)], labels=[ - '$10^{' + str(i) + '}$' for i in range(2, 13)]) - plt.xlabel('Viral load (RNA copies mL$^{-1}$)', fontsize=14) + # Place a legend outside of the axes itself. + ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + ax.set_ylim(0) return fig @@ -212,18 +174,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) @@ -233,8 +191,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() @@ -245,10 +202,8 @@ 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() @@ -265,16 +220,14 @@ 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] + sample_dts = [datetime(1970, 1, 1) + timedelta(hours=time) for time in sample_times] for name, statistics in scenarios.items(): concentrations = statistics['concentrations'] if name in dash_styled_scenarios: ax.plot(sample_dts, concentrations, label=name, linestyle='--') else: - ax.plot(sample_dts, concentrations, - label=name, linestyle='-', alpha=0.5) + ax.plot(sample_dts, concentrations, label=name, linestyle='-', alpha=0.5) # Place a legend outside of the axes itself. ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') @@ -335,8 +288,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( @@ -363,15 +315,14 @@ 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', - 'threshold': '' - } + 'onsite_access': 'of about 8000', + 'threshold' : '' + } return context def _template_environment(self) -> jinja2.Environment: @@ -389,3 +340,4 @@ class ReportGenerator: def render(self, context: dict) -> str: template = self._template_environment().get_template("calculator.report.html.j2") return template.render(**context) +