From 473f50d1e47982b3950d94fd48d71f9fdc0c01fd Mon Sep 17 00:00:00 2001 From: markus Date: Wed, 17 Feb 2021 16:49:33 +0100 Subject: [PATCH] add generate_cdf_curves_vs_qr --- cara/montecarlo.py | 107 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 82 insertions(+), 25 deletions(-) diff --git a/cara/montecarlo.py b/cara/montecarlo.py index d1555a58..53b88c13 100644 --- a/cara/montecarlo.py +++ b/cara/montecarlo.py @@ -9,6 +9,7 @@ from scipy.integrate import quad import typing import matplotlib.pyplot as plt import matplotlib.patches as patches +import matplotlib.lines as mlines from sklearn.neighbors import KernelDensity # This is a test comment @@ -635,6 +636,61 @@ def generate_boxplot(masked: bool = False, samples: int = 200000, qid: int = 100 plt.show() +def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid: int = 100) -> None: + """ + Generates and displays a boxplot for comparing the qR-values in scenarios with various combinations of breathing + rates and expiratory activities + :param masked: A bool indicating whether or not the infected subject is wearing a mask + :param samples: The number of samples to use for the monte carlo simulation + :param qid: The number of "viral copies per quantum" + :return: Nothing, a graph is displayed + """ + fig, axs = plt.subplots(3, 1, sharex='all') + + # TODO: Insert title and y-label + plt.suptitle("INSERT TITLE HERE") + fig.text(0.02, 0.5, 'INSERT Y-LABEL HERE', va='center', rotation='vertical') + + scenarios = [MCInfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 4), (5, 9))), + masked=masked, + virus=MCVirus(halflife=1.1, qID=qid), + expiratory_activity=e, + samples=samples, + breathing_category=b, + ) for e, b in product((1, 2, 3), (1, 3, 5))] + qr_values = [np.log10(scenario.emission_rate_when_present()) for scenario in scenarios] + left = min(np.min(qrs) for qrs in qr_values) + right = max(np.max(qrs) for qrs in qr_values) + + # Colors can be changed here + colors = ['lightgreen', 'lightblue', 'pink'] + + breathing_rates = ['Seated', 'Light activity', 'Heavy activity'] + activities = ['Breathing', 'Speaking', 'Shouting'] + lines = [mlines.Line2D([], [], color=color, markersize=15, label=label) + for color, label in zip(colors, breathing_rates)] + + for i in range(3): + axs[i].hist(qr_values[3 * i:3 * (i + 1)], bins=2000, histtype='step', + color=colors, cumulative=True, range=(left, right)) + axs[i].set_xlim(left, right) + axs[i].set_yticks([0, samples / 2, samples]) + axs[i].set_yticklabels(['0.0', '0.5', '1.0']) + axs[i].yaxis.set_label_position("right") + axs[i].set_ylabel(activities[i]) + axs[i].legend(handles=lines, loc='upper left') + + plt.xlabel('qR [$q\;h^{-1}$]') + tick_positions = np.arange(int(np.ceil(left)), int(np.ceil(right)), 2) + plt.xticks(ticks=tick_positions, labels=['$\;10^{' + str(i) + '}$' for i in tick_positions]) + + fig.set_figheight(8) + fig.set_figwidth(5) + plt.show() + + def buonanno_exposure_model() -> MCExposureModel: """ :return: An MCExposureModel object corresponding to one of the scenarios outlined in the paper of Buonanno et al. @@ -814,30 +870,31 @@ exposure_models = [MCExposureModel( ) ) for qid in (100, 60)] +generate_cdf_curves_vs_qr(masked=False) -rs = [model.expected_new_cases() for model in large_population_baselines] - -print(f"R0 - original variant:\t{np.mean(rs[0])}") -print(f"R0 - english variant:\t{np.mean(rs[1])}") -print(f"Ratio between R0's:\t\t{np.mean(rs[1]) / np.mean(rs[0])}") - +# rs = [model.expected_new_cases() for model in large_population_baselines] +# +# print(f"R0 - original variant:\t{np.mean(rs[0])}") +# print(f"R0 - english variant:\t{np.mean(rs[1])}") +# print(f"Ratio between R0's:\t\t{np.mean(rs[1]) / np.mean(rs[0])}") +# # compare_infection_probabilities_vs_viral_loads(*exposure_models) - - -present_model(exposure_models[0].concentration_model) -plot_pi_vs_qid(fixed_vl_exposure_models, labels=['Viral load = $10^{' + str(i) + '}$' for i in range(6, 11)], - qid_min=5, qid_max=2000, qid_samples=200) - -plot_pi_vs_qid(fixed_vl_exposure_models, labels=['Viral load = $10^{' + str(i) + '}$' for i in range(6, 11)], - qid_min=100, qid_max=400, qid_samples=100) - - -plot_pi_vs_viral_load(exposure_models, labels=['Without masks', 'With masks']) - -for model in exposure_models: - present_model(model.concentration_model, title=f'Model summary - {"English" if model.concentration_model.infected.qid == 60 else "Original"} variant') - plt.hist(model.infection_probability(), bins=200) - plt.xlabel('Percentage probability of infection') - plt.title(f'Probability of infection in baseline case - {"English" if model.concentration_model.infected.qid == 60 else "Original"} variant') - plt.yticks([], []) - plt.show() +# +# +# present_model(exposure_models[0].concentration_model) +# plot_pi_vs_qid(fixed_vl_exposure_models, labels=['Viral load = $10^{' + str(i) + '}$' for i in range(6, 11)], +# qid_min=5, qid_max=2000, qid_samples=200) +# +# plot_pi_vs_qid(fixed_vl_exposure_models, labels=['Viral load = $10^{' + str(i) + '}$' for i in range(6, 11)], +# qid_min=100, qid_max=400, qid_samples=100) +# +# +# plot_pi_vs_viral_load(exposure_models, labels=['Without masks', 'With masks']) +# +# for model in exposure_models: +# present_model(model.concentration_model, title=f'Model summary - {"English" if model.concentration_model.infected.qid == 60 else "Original"} variant') +# plt.hist(model.infection_probability(), bins=200) +# plt.xlabel('Percentage probability of infection') +# plt.title(f'Probability of infection in baseline case - {"English" if model.concentration_model.infected.qid == 60 else "Original"} variant') +# plt.yticks([], []) +# plt.show()