diff --git a/caimira/apps/calculator/__init__.py b/caimira/apps/calculator/__init__.py index 0b014d59..e0a2fd78 100644 --- a/caimira/apps/calculator/__init__.py +++ b/caimira/apps/calculator/__init__.py @@ -38,7 +38,7 @@ from .user import AuthenticatedUser, AnonymousUser # calculator version. If the calculator needs to make breaking changes (e.g. change # form attributes) then it can also increase its MAJOR version without needing to # increase the overall CAiMIRA version (found at ``caimira.__version__``). -__version__ = "4.12" +__version__ = "4.12.1" LOG = logging.getLogger(__name__) diff --git a/caimira/apps/calculator/report_generator.py b/caimira/apps/calculator/report_generator.py index b64b64b4..b4510b3f 100644 --- a/caimira/apps/calculator/report_generator.py +++ b/caimira/apps/calculator/report_generator.py @@ -143,8 +143,11 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing prob_dist_count, prob_dist_bins = np.histogram(prob/100, bins=100, density=True) prob_probabilistic_exposure = np.array(model.total_probability_rule()).mean() expected_new_cases = np.array(model.expected_new_cases()).mean() - uncertainties_plot_src = img2base64(_figure2bytes(uncertainties_plot(model))) if form.conditional_probability_plot else None + uncertainties_plot_src = img2base64(_figure2bytes(uncertainties_plot(model, prob))) if form.conditional_probability_plot else None exposed_presence_intervals = [list(interval) for interval in model.exposed.presence_interval().boundaries()] + conditional_probability_data = {key: value for key, value in + zip(('viral_loads', 'pi_means', 'lower_percentiles', 'upper_percentiles'), + manufacture_conditional_probability_data(model, prob))} return { "model_repr": repr(model), @@ -166,6 +169,7 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing "uncertainties_plot_src": uncertainties_plot_src, "CO2_concentrations": CO2_concentrations, "vl_dist": list(np.log10(model.concentration_model.virus.viral_load_in_sputum)), + "conditional_probability_data": conditional_probability_data, } @@ -188,24 +192,38 @@ def generate_permalink(base_url, get_root_url, get_root_calculator_url, form: F } -def uncertainties_plot(exposure_model: models.ExposureModel): +def conditional_prob_inf_given_vl_dist(infection_probability: models._VectorisedFloat, + viral_loads: np.ndarray, specific_vl: float, step: models._VectorisedFloat): + pi_means = [] + lower_percentiles = [] + upper_percentiles = [] + + for vl_log in viral_loads: + specific_prob = infection_probability[np.where((vl_log-step/2-specific_vl)*(vl_log+step/2-specific_vl)<0)[0]] #type: ignore + pi_means.append(specific_prob.mean()) + lower_percentiles.append(np.quantile(specific_prob, 0.05)) + upper_percentiles.append(np.quantile(specific_prob, 0.95)) + + return pi_means, lower_percentiles, upper_percentiles + + +def manufacture_conditional_probability_data(exposure_model: models.ExposureModel, + infection_probability: models._VectorisedFloat): + + min_vl, max_vl, step = 2, 10, 8/100 + viral_loads = np.arange(min_vl, max_vl, step) + specific_vl = np.log10(exposure_model.concentration_model.virus.viral_load_in_sputum) + pi_means, lower_percentiles, upper_percentiles = conditional_prob_inf_given_vl_dist(infection_probability, viral_loads, + specific_vl, step) + + return list(viral_loads), list(pi_means), list(lower_percentiles), list(upper_percentiles) + + +def uncertainties_plot(exposure_model: models.ExposureModel, prob: models._VectorisedFloat): fig = plt.figure(figsize=(4, 7), dpi=110) - - viral_loads = np.linspace(2, 10, 600) - pi_means, lower_percentiles, upper_percentiles = [], [], [] - for vl in viral_loads: - model_vl = dataclass_utils.nested_replace( - exposure_model, { - 'concentration_model.infected.virus.viral_load_in_sputum' : 10**vl, - } - ) - pi = model_vl.infection_probability()/100 - - pi_means.append(np.mean(pi)) - lower_percentiles.append(np.quantile(pi, 0.05)) - upper_percentiles.append(np.quantile(pi, 0.95)) - - histogram_data = exposure_model.infection_probability() / 100 + + infection_probability = prob / 100 + viral_loads, pi_means, lower_percentiles, upper_percentiles = manufacture_conditional_probability_data(exposure_model, infection_probability) fig, axs = plt.subplots(2, 3, gridspec_kw={'width_ratios': [5, 0.5] + [1], @@ -221,7 +239,7 @@ def uncertainties_plot(exposure_model: models.ExposureModel): axs[0, 0].plot(viral_loads, pi_means, label='Predictive total probability') axs[0, 0].fill_between(viral_loads, lower_percentiles, upper_percentiles, alpha=0.1, label='5ᵗʰ and 95ᵗʰ percentile') - axs[0, 2].hist(histogram_data, bins=30, orientation='horizontal') + axs[0, 2].hist(infection_probability, bins=30, orientation='horizontal') axs[0, 2].set_xticks([]) axs[0, 2].set_xticklabels([]) axs[0, 2].set_facecolor("lightgrey") @@ -230,7 +248,7 @@ def uncertainties_plot(exposure_model: models.ExposureModel): axs[0, 2].set_xlim(0, highest_bar) axs[0, 2].text(highest_bar * 0.5, 0.5, - rf"$\bf{np.round(np.mean(histogram_data) * 100, 1)}$%", ha='center', va='center') + rf"$\bf{np.round(np.mean(infection_probability) * 100, 1)}$%", ha='center', va='center') axs[1, 0].hist(np.log10(exposure_model.concentration_model.infected.virus.viral_load_in_sputum), bins=150, range=(2, 10), color='grey') axs[1, 0].set_facecolor("lightgrey") @@ -310,9 +328,9 @@ def non_zero_percentage(percentage: int) -> str: return "99.9%" else: return "{:0.1f}%".format(percentage) - -def manufacture_viral_load_scenarios(model: mc.ExposureModel) -> typing.Dict[str, mc.ExposureModel]: + +def manufacture_viral_load_scenarios_percentiles(model: mc.ExposureModel) -> typing.Dict[str, mc.ExposureModel]: viral_load = model.concentration_model.infected.virus.viral_load_in_sputum scenarios = {} for percentil in (0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99): @@ -320,7 +338,7 @@ def manufacture_viral_load_scenarios(model: mc.ExposureModel) -> typing.Dict[str specific_vl_scenario = dataclass_utils.nested_replace(model, {'concentration_model.infected.virus.viral_load_in_sputum': vl} ) - scenarios[round(np.log10(vl))] = np.mean(specific_vl_scenario.infection_probability()) + scenarios[vl] = np.mean(specific_vl_scenario.infection_probability()) return scenarios @@ -471,7 +489,7 @@ class ReportGenerator: context.update(report_data) alternative_scenarios = manufacture_alternative_scenarios(form) - context['alternative_viral_load'] = manufacture_viral_load_scenarios(model) if form.conditional_probability_viral_loads else None + context['alternative_viral_load'] = manufacture_viral_load_scenarios_percentiles(model) if form.conditional_probability_viral_loads else None context['alternative_scenarios'] = comparison_report( form, report_data, alternative_scenarios, scenario_sample_times, executor_factory=executor_factory, ) diff --git a/caimira/apps/templates/base/calculator.form.html.j2 b/caimira/apps/templates/base/calculator.form.html.j2 index 2ab873b3..8277f829 100644 --- a/caimira/apps/templates/base/calculator.form.html.j2 +++ b/caimira/apps/templates/base/calculator.form.html.j2 @@ -408,7 +408,6 @@
-
diff --git a/caimira/tests/test_conditional_probability.py b/caimira/tests/test_conditional_probability.py new file mode 100644 index 00000000..85f1c666 --- /dev/null +++ b/caimira/tests/test_conditional_probability.py @@ -0,0 +1,77 @@ +import numpy as np +import pytest +from retry import retry + +import caimira.monte_carlo as mc +from caimira import models +from caimira.dataclass_utils import nested_replace +from caimira.apps.calculator import report_generator +from caimira.monte_carlo.data import activity_distributions, virus_distributions, expiration_distributions + + +@pytest.fixture +def baseline_exposure_model(): + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=50, inside_temp=models.PiecewiseConstant((0., 24.), (298,)), humidity=0.5), + ventilation=models.MultipleVentilation( + ventilations=( + models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25), + ) + ), + infected=mc.InfectedPopulation( + number=1, + presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))), + virus=virus_distributions['SARS_CoV_2_DELTA'], + mask=models.Mask.types['No mask'], + activity=activity_distributions['Seated'], + expiration=expiration_distributions['Breathing'], + host_immunity=0., + ), + evaporation_factor=0.3, + ) + return mc.ExposureModel( + concentration_model=concentration_mc, + short_range=(), + exposed=mc.Population( + number=3, + presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))), + activity=activity_distributions['Seated'], + mask=models.Mask.types['No mask'], + host_immunity=0., + ), + geographical_data=models.Cases(), + ) + + +@retry(tries=3) +def test_conditional_prob_inf_given_vl_dist(baseline_exposure_model): + + viral_loads = np.array([3., 5., 7., 9.,]) + mc_model: models.ExposureModel = baseline_exposure_model.build_model(2_000_000) + + expected_pi_means = [] + expected_lower_percentiles = [] + expected_upper_percentiles = [] + + for vl in viral_loads: + model_vl: models.ExposureModel = nested_replace( + mc_model, { + 'concentration_model.infected.virus.viral_load_in_sputum' : 10**vl, + } + ) + pi = model_vl.infection_probability()/100 + + expected_pi_means.append(np.mean(pi)) + expected_lower_percentiles.append(np.quantile(pi, 0.05)) + expected_upper_percentiles.append(np.quantile(pi, 0.95)) + + infection_probability = mc_model.infection_probability() / 100 + specific_vl = np.log10(mc_model.concentration_model.infected.virus.viral_load_in_sputum) + step = 8/100 + actual_pi_means, actual_lower_percentiles, actual_upper_percentiles = ( + report_generator.conditional_prob_inf_given_vl_dist(infection_probability, viral_loads, specific_vl, step) + ) + + assert np.allclose(actual_pi_means, expected_pi_means, atol=0.002) + assert np.allclose(actual_lower_percentiles, expected_lower_percentiles, atol=0.002) + assert np.allclose(actual_upper_percentiles, expected_upper_percentiles, atol=0.002)