breathing, speaking and shouting histograms all in one

This commit is contained in:
Luis Aleixo 2021-09-06 14:44:14 +02:00
parent cca80d86ea
commit e152b296c8
2 changed files with 34 additions and 0 deletions

View file

@ -7,6 +7,7 @@ Availability: <where it's located> """
from cara import model_scenarios_paper
from cara.results_paper import *
from cara.test_plots import *
from cara.monte_carlo.data import symptomatic_vl_frequencies
# Exhaled virions while talking, seated #
print('\n<<<<<<<<<<< Vlout for Talking, seated >>>>>>>>>>>')
@ -24,12 +25,14 @@ print('\n<<<<<<<<<<< Vlout for Talking, seated with chosen Cn,L >>>>>>>>>>>')
print('\n<<<<<<<<<<< Vlout for Breathing, seated with chosen Cn,B >>>>>>>>>>>')
#exposure_model_from_vl_breathing_cn()
print('\n')
############ Statistical Data ############
############ Breathing model ############
exposure_mc = model_scenarios_paper.breathing_exposure()
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
breathing_er = [np.log10(er) for er in emission_rate]
print('\n<<<<<<<<<<< Breathing model statistics >>>>>>>>>>>')
print_er_info(emission_rate)
@ -37,6 +40,7 @@ print_er_info(emission_rate)
exposure_mc = model_scenarios_paper.speaking_exposure()
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
speaking_er = [np.log10(er) for er in emission_rate]
print('\n<<<<<<<<<<< Speaking model statistics >>>>>>>>>>>')
print_er_info(emission_rate)
@ -44,8 +48,15 @@ print_er_info(emission_rate)
exposure_mc = model_scenarios_paper.shouting_exposure()
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
shouting_er = [np.log10(er) for er in emission_rate]
print('\n<<<<<<<<<<< Shouting model statistics >>>>>>>>>>>')
print_er_info(emission_rate)
############ Plots with Viral loads and emission rates
viral_load_in_sputum = exposure_model.concentration_model.infected.virus.viral_load_in_sputum
present_vl_er_histograms(viral_load_in_sputum, breathing_er, speaking_er, shouting_er)
############ Used for testing ############
#exposure_model_from_vl_talking_new_points()

View file

@ -252,6 +252,29 @@ def exposure_model_from_vl_talking_cn():
plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
plt.show()
def present_vl_er_histograms(viral_load_in_sputum, breathing_er, speaking_er, shouting_er):
fig, axs = plt.subplots(1, 2, sharex=False, sharey=False)
plt.tight_layout()
viral_loads = [np.log10(vl) for vl in viral_load_in_sputum]
axs[0].hist(viral_loads, bins = 200)
axs[0].title.set_text('Viral load')
axs[0].set_xlabel('vl (log$_{10}$(RNA copies mL$^{-1}$))')
axs[1].title.set_text('Viral emission rate')
axs[1].hist(breathing_er, bins = 200, label='Breathing vR', alpha=0.5)
axs[1].hist(speaking_er, bins = 200, label='Speaking vR', alpha=0.5)
axs[1].hist(shouting_er, bins = 200, label='Shouting vR', alpha=0.5)
axs[1].set_xlabel('vR (log$_{10}$)')
for x in (0, 1):
axs[x].set_yticklabels([])
axs[x].set_yticks([])
plt.legend(loc='upper right')
plt.show()
######### Auxiliar functions #########