new plot with multiple cn lines

cara/model_scenarios.py

@@ -1,3 +1,4 @@
+from numpy.core.function_base import linspace
 from cara import models
 from cara.monte_carlo.data import activity_distributions
 from tqdm import tqdm
@@ -8,6 +9,7 @@ from scipy.spatial import ConvexHull
 import pandas as pd
 import matplotlib.lines as mlines
 from matplotlib.patches import Rectangle
+import matplotlib as mpl
 
 ######### Plot material #########
 fig = plt.figure()
@@ -128,6 +130,104 @@ def exposure_model_from_vl_talking(viral_loads):
 
     return er_means
 
+def exposure_model_from_vl_talking_cn(viral_loads):
+    
+    n_lines = 5
+    cns = np.linspace(0., .2, n_lines)
+    norm = mpl.colors.Normalize(vmin=cns.min(), vmax=cns.max())
+    cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.jet)
+    cmap.set_array([])
+
+    for cn in tqdm(cns):
+        er_means = []
+        er_medians = []
+        lower_percentiles = []
+        upper_percentiles = []
+        for vl in 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.Expiration.types['Talking'],
+                    ),
+                ),
+                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)
+            # divide by 4 to have in 15min (quarter of an hour)
+            emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(cn)/4
+            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))
+
+        # divide by 4 to have in 15min (quarter of an hour)
+        coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
+
+        
+        ax.plot(viral_loads, er_means, color=cmap.to_rgba(cn))
+        ax.fill_between(viral_loads, lower_percentiles,
+                        upper_percentiles, alpha=0.2, color=cmap.to_rgba(cn))
+    
+    fig.colorbar(cmap, ticks=cns)
+    ax.set_yscale('log')
+
+    ############# Coleman #############
+    plt.scatter(coleman_etal_vl_talking, coleman_etal_er_talking_2, marker='x')
+    x_hull, y_hull = get_enclosure_points(
+        coleman_etal_vl_talking, coleman_etal_er_talking_2)
+    # plot shape
+    plt.fill(x_hull, y_hull, '--', c='orange', alpha=0.2)
+
+    ############# Markers #############
+    markers = [5, 'd', 4]
+
+    ############ Legend ############
+    result_from_model = mlines.Line2D(
+        [], [], color='blue', marker='_', linestyle='None')
+    coleman = mlines.Line2D([], [], color='orange',
+                            marker='x', linestyle='None')
+
+    title_proxy = Rectangle((0, 0), 0, 0, color='w')
+    titles = ["$\\bf{CARA \, \\it{(SARS-CoV-2)}:}$", "$\\bf{Coleman \, et \, al. \, \\it{(SARS-CoV-2)}:}$"]
+    leg = plt.legend([title_proxy, result_from_model, title_proxy, coleman],
+                     [titles[0], "Result from model", titles[1], "Dataset"])
+
+    # Move titles to the left
+    for item, label in zip(leg.legendHandles, leg.texts):
+        if label._text in titles:
+            width = item.get_window_extent(fig.canvas.get_renderer()).width
+            label.set_ha('left')
+            label.set_position((-3*width, 0))
+
+    ############ Plot ############
+    plt.title('Exhaled virions while talking for 15min',
+              fontsize=16, fontweight="bold")
+    plt.ylabel(
+        'Aerosol viral load, $\mathrm{vl_{out}}$\n(RNA copies)', fontsize=14)
+    plt.xticks(ticks=[i for i in range(2, 13)], labels=[
+        '$10^{' + str(i) + '}$' for i in range(2, 13)])
+    plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
+    plt.show()
+
+    return er_means
 
 def exposure_model_from_vl_breathing(viral_loads):
 

cara/models.py

@@ -555,7 +555,7 @@ class Expiration(_ExpirationBase):
     BLO_factors: typing.Tuple[float, float, float]
 
     @cached()
-    def aerosols(self, mask: Mask):
+    def aerosols(self, mask: Mask, cn: float):
         """ Result is in mL.cm^-3 """
         def volume(d):
             return (np.pi * d**3) / 6.
@@ -565,9 +565,9 @@ class Expiration(_ExpirationBase):
             return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) *
                     np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2)))
 
-        def _Lmode(d: float) -> float:
+        def _Lmode(d: float, cn: float) -> float:
             # L-mode (see ref. above).
-            return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.506818)) *
+            return ( (1 / d) * (cn / (np.sqrt(2 * np.pi) * 0.506818)) *
                     np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2)))
 
         def _Omode(d: float) -> float:
@@ -577,7 +577,7 @@ class Expiration(_ExpirationBase):
 
         def integrand(d: float) -> float:
             return (self.BLO_factors[0] * _Bmode(d) +
-                    self.BLO_factors[1] * _Lmode(d) +
+                    self.BLO_factors[1] * _Lmode(d, cn) +
                     self.BLO_factors[2] * _Omode(d)
                     ) * volume(d) * (1 - mask.exhale_efficiency(d))
 
@@ -670,7 +670,7 @@ class InfectedPopulation(Population):
     #: The type of expiration that is being emitted whilst doing the activity.
     expiration: _ExpirationBase
 
-    def emission_rate_when_present(self) -> _VectorisedFloat:
+    def emission_rate_when_present(self, cn: float) -> _VectorisedFloat:
         """
         The emission rate if the infected population is present.
 
@@ -680,7 +680,7 @@ class InfectedPopulation(Population):
         # Emission Rate (virions / h)
         # Note on units: exhalation rate is in m^3/h, aerosols in mL/cm^3
         # and viral load in virus/mL -> 1e6 conversion factor
-        aerosols = self.expiration.aerosols(self.mask)
+        aerosols = self.expiration.aerosols(self.mask, cn)
 
         ER = (self.virus.viral_load_in_sputum *
               self.activity.exhalation_rate *

cara/plot_output.py

@@ -1,11 +1,12 @@
-from cara.model_scenarios import exposure_model_from_vl_breathing, exposure_model_from_vl_talking
+from cara.model_scenarios import exposure_model_from_vl_breathing, exposure_model_from_vl_talking, exposure_model_from_vl_talking_cn
 import numpy as np
 import csv
 
 viral_loads = np.linspace(2, 12, 600)
 
 #er_means = exposure_model_from_vl_talking(viral_loads)
-er_means = exposure_model_from_vl_breathing(viral_loads)
+#er_means = exposure_model_from_vl_breathing(viral_loads)
+er_means = exposure_model_from_vl_talking_cn(viral_loads)
 
 with open('data.csv', 'w', newline='') as csvfile:
     fieldnames = ['viral load', 'emission rate']