From 646b005baa0e218f9f249b4663dd7783d5f40900 Mon Sep 17 00:00:00 2001 From: markus Date: Wed, 24 Feb 2021 13:18:20 +0100 Subject: [PATCH 1/4] fix compare_concentration_curves --- cara/montecarlo.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 7f72b870..2cf996f6 100644 --- a/cara/montecarlo.py +++ b/cara/montecarlo.py @@ -1036,16 +1036,19 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label times = np.arange(start, stop, TIME_STEP) - concentrations = [[np.mean(model.concentration_model.concentration(t)) for t in times] for model in exp_models] + concentrations = [[np.median(model.concentration_model.concentration(t)) for t in times] for model in exp_models] for c, label, color in zip(concentrations, labels, colors): plt.plot(times, c, label=label, color=color) plt.legend() + factors = [0.6 * model.exposed.activity.inhalation_rate * (1 - model.exposed.mask.η_inhale) for model in exp_models] + qds = [[np.trapz(c[:i + 1], times[:i + 1]) * factor for i in range(len(times))] + for c, factor in zip(concentrations, factors)] if show_qd: plt.twinx() - for c, label, color in zip(concentrations, labels, colors): - plt.plot(times, np.cumsum(c), label='qD - ' + label, color=color, linestyle='dotted') + for qd, label, color in zip(qds, labels, colors): + plt.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted') plt.legend() plt.show() From fdf4f38bbfb9e15913e9213fff634b38182e69db Mon Sep 17 00:00:00 2001 From: markus Date: Wed, 24 Feb 2021 13:19:56 +0100 Subject: [PATCH 2/4] update usage example --- cara/mc-output.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/cara/mc-output.py b/cara/mc-output.py index 9651e468..1a60e9bb 100644 --- a/cara/mc-output.py +++ b/cara/mc-output.py @@ -1,7 +1,8 @@ from cara.montecarlo import * from cara.model_scenarios import * -# compare_concentration_curves([classroom_model, classroom_model_with_hepa], ['Just window', 'Window and HEPA']) +classroom_model.infection_probability() +compare_concentration_curves([classroom_model, classroom_model_with_hepa], ['Just window', 'Window and HEPA']) #print(np.mean(chorale_model.infection_probability())) #print(np.mean(chorale_model.infection_probability())+np.std(chorale_model.infection_probability())) @@ -11,12 +12,12 @@ from cara.model_scenarios import * -#print(np.mean(classroom_model_with_hepa.infection_probability())) -composite_plot_pi_vs_viral_load([shared_office_worst_model[1], shared_office_model[1],shared_office_better_model[1]], - labels=['No mask &\nwindows closed', 'Baseline', 'Baseline +\nHEPA filter'], - colors=['tomato', '#1f77b4', 'limegreen'], - title='$P(I|qID)$ vs $vl$ - Shared office scenario', - vl_points=200) +# #print(np.mean(classroom_model_with_hepa.infection_probability())) +# composite_plot_pi_vs_viral_load([shared_office_worst_model[1], shared_office_model[1],shared_office_better_model[1]], +# labels=['No mask &\nwindows closed', 'Baseline', 'Baseline +\nHEPA filter'], +# colors=['tomato', '#1f77b4', 'limegreen'], +# title='$P(I|qID)$ vs $vl$ - Shared office scenario', +# vl_points=200) #plot_pi_vs_viral_load([shared_office_model[1]], labels=['Baseline, qID=60', 'HEPA, qID=60', 'No mask + windows closed, qID=60'],title='$P(I|qID)$ - Shared office scenario') From fac33a5cf12eea8b0e96030179483a71b5fb5868 Mon Sep 17 00:00:00 2001 From: markus Date: Wed, 24 Feb 2021 14:32:37 +0100 Subject: [PATCH 3/4] add tertiary axis and labels --- cara/montecarlo.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 2cf996f6..02e04f8c 100644 --- a/cara/montecarlo.py +++ b/cara/montecarlo.py @@ -1037,18 +1037,32 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label times = np.arange(start, stop, TIME_STEP) concentrations = [[np.median(model.concentration_model.concentration(t)) for t in times] for model in exp_models] + fig, ax = plt.subplots() for c, label, color in zip(concentrations, labels, colors): - plt.plot(times, c, label=label, color=color) + ax.plot(times, c, label=label, color=color) + + ax.legend(loc='upper left') - plt.legend() factors = [0.6 * model.exposed.activity.inhalation_rate * (1 - model.exposed.mask.η_inhale) for model in exp_models] qds = [[np.trapz(c[:i + 1], times[:i + 1]) * factor for i in range(len(times))] for c, factor in zip(concentrations, factors)] + plt.suptitle("SUPTITLE HERE") + plt.xlabel("XLABEL HERE") + plt.ylabel("YLABEL HERE") if show_qd: - plt.twinx() - for qd, label, color in zip(qds, labels, colors): - plt.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted') - plt.legend() + ax1 = ax.twinx() + for qd, label, color in zip(qds, labels, colors): + ax1.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted') + ax1.set_ylabel('qD LABEL HERE') + + ax2 = ax.twinx() + ax2.spines["right"].set_position(("axes", 1.2)) + ax2.set_ylabel('RNA LABEL HERE') + ax2.set_ylim(tuple(y * exp_models[0].concentration_model.virus.qID for y in ax1.get_ylim())) + ax1.legend(loc='upper right') + + plt.tight_layout() + plt.show() From f19d85ae8b452b40ecfc4186efba2025f5414fa0 Mon Sep 17 00:00:00 2001 From: markus Date: Wed, 24 Feb 2021 14:53:24 +0100 Subject: [PATCH 4/4] scale y-axes --- cara/montecarlo.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 02e04f8c..fe0f4c27 100644 --- a/cara/montecarlo.py +++ b/cara/montecarlo.py @@ -1042,6 +1042,7 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label ax.plot(times, c, label=label, color=color) ax.legend(loc='upper left') + ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2) factors = [0.6 * model.exposed.activity.inhalation_rate * (1 - model.exposed.mask.η_inhale) for model in exp_models] qds = [[np.trapz(c[:i + 1], times[:i + 1]) * factor for i in range(len(times))] @@ -1056,6 +1057,7 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label for qd, label, color in zip(qds, labels, colors): ax1.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted') ax1.set_ylabel('qD LABEL HERE') + ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.2) ax2 = ax.twinx() ax2.spines["right"].set_position(("axes", 1.2))