From 81dcbcd39e5fe426ffe68acae0784e97e524580d Mon Sep 17 00:00:00 2001
From: Luis Aleixo <luis.aleixo@cern.ch>
Date: Thu, 5 Aug 2021 15:10:21 +0200
Subject: [PATCH] file used to plot and write data

---
 cara/apps/calculator/report_generator.py |  3 +-
 cara/plot_output.py                      | 78 +++++++++++++++++++++++-
 2 files changed, 79 insertions(+), 2 deletions(-)

diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py
index 1cae1a35..859f8e49 100644
--- a/cara/apps/calculator/report_generator.py
+++ b/cara/apps/calculator/report_generator.py
@@ -130,7 +130,7 @@ def plot(times, concentrations, model: models.ExposureModel):
                 ),
                 infected=mc.InfectedPopulation(
                     number=1,
-                    virus=mc.Virus(
+                    virus=models.Virus(
                         viral_load_in_sputum = 10**vl,
                         infectious_dose = 50.,
                     ),
@@ -151,6 +151,7 @@ def plot(times, concentrations, model: models.ExposureModel):
             ),
         )
         exposure_model = exposure_mc.build_model(size=_DEFAULT_MC_SAMPLE_SIZE)
+        print(exposure_model)
         
         emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present()
         er_means.append(np.mean(emission_rate))
diff --git a/cara/plot_output.py b/cara/plot_output.py
index 9e02d9b8..e74740d3 100644
--- a/cara/plot_output.py
+++ b/cara/plot_output.py
@@ -1,2 +1,78 @@
-from cara import *
+from dataclasses import field
+import numpy as np
+import matplotlib.pyplot as plt
+import csv
 
+import cara.monte_carlo as mc
+from cara import models,data
+from cara.monte_carlo.data import activity_distributions
+from tqdm import tqdm
+
+SAMPLE_SIZE = 50000
+
+fig = plt.figure()
+ax = fig.add_subplot(1, 1, 1)
+
+points = 600
+viral_loads = np.linspace(2, 12, points)
+
+er_means = []
+er_medians = []
+lower_percentiles = []
+upper_percentiles = []
+
+for vl in tqdm(viral_loads):
+    exposure_mc = mc.ExposureModel(
+        concentration_model=mc.ConcentrationModel(
+            room=models.Room(volume=100, humidity=0.5),
+            ventilation=models.AirChange(
+                active=models.SpecificInterval(((0, 24),)),
+                air_exch=0.25,
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                virus=models.Virus(
+                    viral_load_in_sputum = 10**vl,
+                    infectious_dose = 50.,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Seated'],
+                expiration=models.MultipleExpiration(
+                    expirations=(models.Expiration.types['Talking'],
+                                    models.Expiration.types['Breathing']),
+                    weights=(0.3, 0.7)),
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Seated'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+    emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present()
+    er_means.append(np.mean(emission_rate))
+    er_medians.append(np.median(emission_rate))
+    lower_percentiles.append(np.quantile(emission_rate, 0.01))
+    upper_percentiles.append(np.quantile(emission_rate, 0.99))
+
+with open('data.csv', 'w', newline='') as csvfile:
+    fieldnames = ['viral load', 'emission rate']
+    thewriter = csv.DictWriter(csvfile, fieldnames=fieldnames)
+    thewriter.writeheader()
+    for i, vl in enumerate(viral_loads):
+        thewriter.writerow({'viral load' : 10**vl, 'emission rate' : er_means[i]})
+
+
+ax.plot(viral_loads, er_means)
+ax.fill_between(viral_loads, lower_percentiles,
+                upper_percentiles, alpha=0.2)
+ax.set_yscale('log')
+plt.title('Emission rate vs Viral Load')
+plt.ylabel('ER (Virions/h)', fontsize=14)
+plt.xticks(ticks=[i for i in range(2, 13)], labels=[
+            '$10^{' + str(i) + '}$' for i in range(2, 13)])
+plt.xlabel('Viral load (RNA copies mL$^{-1}$)', fontsize=14)
+plt.show()