diff --git a/Short-Range_paper/results_paper.py b/Short-Range_paper/results_paper.py index 4f6a6c2a8815d6bc5545883874ad3c1e187e8a82..685a69fc183774e11257d230aae70e49230be529 100644 --- a/Short-Range_paper/results_paper.py +++ b/Short-Range_paper/results_paper.py @@ -43,6 +43,8 @@ from model_scenarios_paper import * ############# FIGURES ############ +# Plots with viral loads and emission rates + statistical data ############ +# present_vl_er_histograms(physical_activity='Standing', expiratory_activities=['Breathing', 'Speaking', 'Shouting']) # Compare covid_overal_vl_data vs. symptomatic_vl_frequencies distributions # compare_vl_hist() diff --git a/Short-Range_paper/scripts_paper.py b/Short-Range_paper/scripts_paper.py index 628bf6e48bd9cc677f323a2ff2dc4f49035d46a0..e059e5d5b358167556f4f2d75805afc0a58460e3 100644 --- a/Short-Range_paper/scripts_paper.py +++ b/Short-Range_paper/scripts_paper.py @@ -18,6 +18,7 @@ import matplotlib.pyplot as plt from matplotlib.lines import Line2D import matplotlib as mpl from matplotlib.colors import rgb2hex + import typing import dataclasses @@ -865,8 +866,15 @@ def statistics_for_random_variable(physical_activity: bool = False, return None -def short_range_emission_rate_values(physical_activity: str, expiratory_activity: str) -> None: - SAMPLE_SIZE = 250_000 if expiratory_activity == "Breathing" else 1_000_000 +def short_range_emission_rate_values(physical_activity: str, + expiratory_activity: str, + virus_dist: typing.Optional[models.Virus] = None, + ): + + if virus_dist: + SAMPLE_SIZE = len(virus_dist.viral_load_in_sputum) + else: + SAMPLE_SIZE = 250_000 if expiratory_activity == "Breathing" else 1_000_000 # Activity distribution activity_dist: models.Activity = activity_distributions(data_registry)[physical_activity].build_model(SAMPLE_SIZE) @@ -878,7 +886,8 @@ def short_range_emission_rate_values(physical_activity: str, expiratory_activity short_range_distances_dist: models._VectorisedFloat = short_range_distances(data_registry).generate_samples(SAMPLE_SIZE) # Virus distribution - virus_dist: models.Virus = virus_distributions(data_registry)['SARS_CoV_2_OMICRON'].build_model(SAMPLE_SIZE) + if not virus_dist: + virus_dist: models.Virus = virus_distributions(data_registry)['SARS_CoV_2_OMICRON'].build_model(SAMPLE_SIZE) # Infected population infected: models.InfectedPopulation = models.InfectedPopulation( @@ -908,4 +917,77 @@ def short_range_emission_rate_values(physical_activity: str, expiratory_activity print('\n<<<<<<<<<<< ' + expiratory_activity + ' model statistics >>>>>>>>>>>') print_er_info(emission_rate, log10_emission_rate) - return None + return emission_rate, log10_emission_rate + +def present_vl_er_histograms(physical_activity: str, expiratory_activities: typing.List) -> None: + # Virus distribution + virus_dist: models.Virus = virus_distributions(data_registry)['SARS_CoV_2_OMICRON'].build_model(SAMPLE_SIZE) + + emission_rate_dict: typing.Dict = {} + log10_emission_rate_dict: typing.Dict = {} + for exp_act in expiratory_activities: + emission_rate, log10_emission_rate = short_range_emission_rate_values(physical_activity, exp_act, virus_dist) + emission_rate_dict[exp_act] = emission_rate + log10_emission_rate_dict[exp_act] = log10_emission_rate + + # Viral load + viral_loads = models._VectorisedFloat = [np.log10(vl) for vl in virus_dist.viral_load_in_sputum] + + # Plot material + fig, axs = plt.subplots(1, 2, sharex=False, sharey=False) + fig.set_figheight(5) + fig.set_figwidth(10) + plt.tight_layout() + plt.subplots_adjust(wspace=0.2) + plt.subplots_adjust(top=0.9) + plt.subplots_adjust(bottom=0.15) + + axs[0].hist(viral_loads, bins = 300, color='lightgrey') + axs[0].set_xlabel('vl$_{\mathrm{in}}$ (log$_{10}$ RNA copies mL$^{-1}$)') + + # Viral load mean value + mean = np.mean(viral_loads) + axs[0].vlines(x=(mean), ymin=0, ymax=axs[0].get_ylim()[1], colors=('grey'), linestyles=('dashed')) + + # Plot colors + hist_colors = ['lightsteelblue', 'wheat', 'darkseagreen'] + hist_vlines_colors = ['cornflowerblue', 'goldenrod', 'olivedrab'] + + # Plot the histograms + log10_emission_rate_dict_mean = {} + for exp_act_key, hist_color in zip(log10_emission_rate_dict, hist_colors): + exp_act_value = log10_emission_rate_dict[exp_act_key] + axs[1].hist(exp_act_value, bins=300, color=hist_color) + # Save the mean value + log10_emission_rate_dict_mean[exp_act_key] = np.mean(exp_act_value) + + # Max value for y axis + ymax: float = axs[1].get_ylim()[1] + + # Plot the vlines + for exp_act_key, vline_color in zip(log10_emission_rate_dict_mean, hist_vlines_colors): + exp_act_mean = log10_emission_rate_dict_mean[exp_act_key] + axs[1].vlines(x=exp_act_mean, ymin=0, ymax=ymax, colors=vline_color, alpha=0.75, linestyles='dashed') + + axs[1].set_xlabel('vR (log$_{10}$ virions h$^{-1}$)') + + from matplotlib.patches import Patch + import matplotlib.lines as mlines + + hist_colors_legend: typing.List = ['lightgrey'] + hist_colors + hist_str_legend: typing.List = ['Viral Load'] + expiratory_activities + + labels = [Patch(facecolor=[], edgecolor=[], color=color, label=label) + for color, label in zip(hist_colors_legend, hist_str_legend)] + labels.append(mlines.Line2D([], [], color='black', + marker='', linestyle='dashed', label='Mean')) + + for x in (0, 1): + axs[x].set_yticklabels([]) + axs[x].set_yticks([]) + + plt.legend(handles=labels, loc='upper left', bbox_to_anchor=(1, 1)) + plt.tight_layout() + plt.show() + + return None \ No newline at end of file