From e152b296c84ad64fa656cbb748d2dac29fb0f463 Mon Sep 17 00:00:00 2001
From: Luis Aleixo <luis.aleixo@cern.ch>
Date: Mon, 6 Sep 2021 14:44:14 +0200
Subject: [PATCH] breathing, speaking and shouting histograms all in one

---
 cara/plot_output.py   | 11 +++++++++++
 cara/results_paper.py | 23 +++++++++++++++++++++++
 2 files changed, 34 insertions(+)

diff --git a/cara/plot_output.py b/cara/plot_output.py
index 8cdb94e6..973d0e42 100644
--- a/cara/plot_output.py
+++ b/cara/plot_output.py
@@ -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()
diff --git a/cara/results_paper.py b/cara/results_paper.py
index 9f7443a7..a3ff2617 100644
--- a/cara/results_paper.py
+++ b/cara/results_paper.py
@@ -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 #########