From beb2fd737175d146967cd0f2c64d87dca894c8a9 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Thu, 27 May 2021 12:33:12 +0200
Subject: [PATCH 01/16] Adding tests for MultipleExpiration

---
 cara/tests/test_expiration.py | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 cara/tests/test_expiration.py

diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py
new file mode 100644
index 00000000..7638d613
--- /dev/null
+++ b/cara/tests/test_expiration.py
@@ -0,0 +1,29 @@
+import re
+
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from cara import models
+
+
+def test_multiple_wrong_weight_size():
+    weights = (1., 2., 3.)
+    e_base = models.Expiration((0.084, 0.009, 0.003, 0.002))
+    with pytest.raises(
+            ValueError,
+            match=re.escape("expirations and weigths should contain the"
+                            "same number of elements")
+    ):
+        e = models.MultipleExpiration([e_base, e_base], weights)
+
+
+def test_multiple():
+    weights = (1., 2.)
+    e1 = models.Expiration((0.03, 0.02, 0.01, 0.005))
+    e2 = models.Expiration((0.05, 0.04, 0.03, 0.01))
+    e = models.MultipleExpiration([e1, e2], weights)
+    assert e.aerosols(models.Mask.types['No mask']) == (
+              e1.aerosols(models.Mask.types['No mask'])/3. +
+            2*e2.aerosols(models.Mask.types['No mask'])/3.
+    )

From a8081d8b0b74a6ad134752d1b92e2620e2aa9909 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Thu, 27 May 2021 12:34:14 +0200
Subject: [PATCH 02/16] Adapting tests for model_generator

---
 cara/tests/apps/calculator/test_model_generator.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py
index 6c5050e4..eb216727 100644
--- a/cara/tests/apps/calculator/test_model_generator.py
+++ b/cara/tests/apps/calculator/test_model_generator.py
@@ -8,6 +8,7 @@ from cara.apps.calculator.model_generator import minutes_since_midnight
 from cara import models
 from cara import data
 import numpy as np
+import numpy.testing as npt
 
 
 def test_model_from_dict(baseline_form_data):
@@ -24,10 +25,11 @@ def test_model_from_dict_invalid(baseline_form_data):
 def test_blend_expiration():
     blend = {'Breathing': 2, 'Talking': 1}
     r = model_generator.build_expiration(blend)
+    mask = models.Mask.types['Type I']
     expected = models.Expiration(
         (0.13466666666666668, 0.02866666666666667, 0.004333333333333334, 0.005)
     )
-    assert r == expected
+    npt.assert_almost_equal(r.aerosols(mask), expected.aerosols(mask))
 
 
 def test_ventilation_slidingwindow(baseline_form: model_generator.FormData):

From af730ccd7c249042a3b5c42924acf05c74dccf99 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Thu, 27 May 2021 13:40:46 +0200
Subject: [PATCH 03/16] Introducing _ExpirationBase and MultipleExpiration
 classes; adapting tests and model_generator accordingly (removing now
 obsolete expiration_blend function)

---
 cara/apps/calculator/model_generator.py | 32 +++-------------
 cara/models.py                          | 51 ++++++++++++++++++++++---
 2 files changed, 50 insertions(+), 33 deletions(-)

diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py
index 1ec76a77..0d7506ae 100644
--- a/cara/apps/calculator/model_generator.py
+++ b/cara/apps/calculator/model_generator.py
@@ -537,38 +537,16 @@ class FormData:
         )
 
 
-def build_expiration(expiration_definition) -> models.Expiration:
+def build_expiration(expiration_definition) -> models._ExpirationBase:
     if isinstance(expiration_definition, str):
-        return models.Expiration.types[expiration_definition]
+        return models._ExpirationBase.types[expiration_definition]
     elif isinstance(expiration_definition, dict):
-        return expiration_blend({
-            build_expiration(exp): amount
-            for exp, amount in expiration_definition.items()
-            }
+        return models.MultipleExpiration(
+            tuple([build_expiration(exp) for exp in expiration_definition.keys()]),
+            tuple(expiration_definition.values())
         )
 
 
-def expiration_blend(expiration_weights: typing.Dict[models.Expiration, int]) -> models.Expiration:
-    """
-    Combine together multiple types of Expiration, using a weighted mean to
-    compute their ejection factor and particle sizes.
-
-    """
-    ejection_factor = np.zeros(4)
-    particle_sizes = np.zeros(4)
-
-    total_weight = 0
-    for expiration, weight in expiration_weights.items():
-        total_weight += weight
-        ejection_factor += np.array(expiration.ejection_factor) * weight
-        particle_sizes += np.array(expiration.particle_sizes) * weight
-
-    r_ejection_factor: typing.Tuple[float, float, float, float] = tuple(ejection_factor/total_weight)  # type: ignore
-    r_particle_sizes: typing.Tuple[float, float, float, float] = tuple(particle_sizes/total_weight)  # type: ignore
-
-    return models.Expiration(ejection_factor=r_ejection_factor, particle_sizes=r_particle_sizes)
-
-
 def model_from_form(form: FormData) -> models.ExposureModel:
     # Initializes room with volume either given directly or as product of area and height
     if form.volume_type == 'room_volume_explicit':
diff --git a/cara/models.py b/cara/models.py
index 620e6ca4..bb3ece66 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -564,13 +564,28 @@ _MaskBase.types = {
 
 
 @dataclass(frozen=True)
-class Expiration:
+class _ExpirationBase:
+    """
+    Represents the expiration of aerosols by a person.
+    Subclasses of _ExpirationBase represent different models.
+    """
+    #: Pre-populated examples of Masks.
+    types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
+
+    def aerosols(self, mask: _MaskBase):
+        # total volume of aerosols expired (cm^3).
+        raise NotImplementedError("Subclass must implement")
+
+
+@dataclass(frozen=True)
+class Expiration(_ExpirationBase):
+    """
+    Simple model based on four different sizes of particles emitted,
+    with different ejection factors.
+    """
     ejection_factor: typing.Tuple[float, ...]
     particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4)  # In cm.
 
-    #: Pre-populated examples of Expiration.
-    types: typing.ClassVar[typing.Dict[str, "Expiration"]]
-
     def aerosols(self, mask: _MaskBase):
         def volume(diameter):
             return (4 * np.pi * (diameter/2)**3) / 3
@@ -582,7 +597,31 @@ class Expiration:
         return total
 
 
-Expiration.types = {
+@dataclass(frozen=True)
+class MultipleExpiration(_ExpirationBase):
+    """
+    Represents an expiration of aerosols.
+    Group together different modes of expiration, that represent
+    each the main expiration mode for a certain fraction of time (given by
+    the weights).
+
+    """
+    expirations: typing.Tuple[_ExpirationBase, ...]
+    weights: typing.Tuple[float, ...]
+
+    def __post_init__(self):
+        if len(self.expirations) != len(self.weights):
+            raise ValueError("expirations and weigths should contain the"
+                             "same number of elements")
+
+    def aerosols(self, mask: _MaskBase):
+        return np.array([
+            weight * expiration.aerosols(mask) / sum(self.weights)
+            for weight,expiration in zip(self.weights,self.expirations)
+        ]).sum(axis=0)
+
+
+_ExpirationBase.types = {
     'Breathing': Expiration((0.084, 0.009, 0.003, 0.002)),
     'Whispering': Expiration((0.11, 0.014, 0.004, 0.002)),
     'Talking': Expiration((0.236, 0.068, 0.007, 0.011)),
@@ -638,7 +677,7 @@ class InfectedPopulation(Population):
     virus: Virus
 
     #: The type of expiration that is being emitted whilst doing the activity.
-    expiration: Expiration
+    expiration: _ExpirationBase
 
     def emission_rate_when_present(self) -> _VectorisedFloat:
         """

From 22539e2f9fbb3e77cbbb7a845de90b4030dabff8 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Thu, 27 May 2021 18:17:59 +0200
Subject: [PATCH 04/16] Adding test on BLO expiration

---
 cara/tests/test_expiration.py | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py
index 7638d613..5cfe94b9 100644
--- a/cara/tests/test_expiration.py
+++ b/cara/tests/test_expiration.py
@@ -23,7 +23,15 @@ def test_multiple():
     e1 = models.Expiration((0.03, 0.02, 0.01, 0.005))
     e2 = models.Expiration((0.05, 0.04, 0.03, 0.01))
     e = models.MultipleExpiration([e1, e2], weights)
-    assert e.aerosols(models.Mask.types['No mask']) == (
-              e1.aerosols(models.Mask.types['No mask'])/3. +
-            2*e2.aerosols(models.Mask.types['No mask'])/3.
-    )
+    mask = models.Mask.types['No mask']
+    assert e.aerosols(mask) == e1.aerosols(mask)/3. + 2*e2.aerosols(mask)/3.
+
+
+def test_multiple_BLO():
+    weights = (1., 1.)
+    e1 = models.ExpirationBLO((1., 0., 0.))
+    e2 = models.ExpirationBLO((4., 5., 5.))
+    e_expected = models.ExpirationBLO((2.5, 2.5, 2.5))
+    e = models.MultipleExpiration([e1, e2], weights)
+    mask = models.Mask.types['Type I']
+    assert e_expected.aerosols(mask) == e.aerosols(mask)

From 23510c8f1a1a1d56f5d16042bdb7a1adcb2525f9 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Thu, 27 May 2021 18:26:37 +0200
Subject: [PATCH 05/16] Adding BLO model for expiration (ExpirationBLO class)

---
 cara/models.py | 51 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)

diff --git a/cara/models.py b/cara/models.py
index a7c8d334..d491e0da 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -37,6 +37,7 @@ import typing
 
 import numpy as np
 from scipy.interpolate import interp1d
+from scipy.integrate import quad
 
 if not typing.TYPE_CHECKING:
     from memoization import cached
@@ -621,12 +622,62 @@ class MultipleExpiration(_ExpirationBase):
         ]).sum(axis=0)
 
 
+@dataclass(frozen=True)
+class ExpirationBLO(_ExpirationBase):
+    """
+    BLO model for the expiration (G. Johnson et al., Modality of human
+    expired aerosol size distributions, Journal of Aerosol Science,
+    vol. 42, no. 12, pp. 839 – 851, 2011,
+    https://doi.org/10.1016/j.jaerosci.2011.07.009).
+    All diameters are in cm.
+    """
+    #: factors assigned to resp. the B, L and O modes. Depends on the
+    # kind of expiratory activity (e.g. breathing, speaking, singing,
+    # or shouting).
+    BLO_factors: typing.Tuple[float, float, float]
+
+    def aerosols(self, mask: _MaskBase):
+
+        def volume(diameter):
+            return (np.pi * diameter**3) / 6.
+
+        def _Bmode(diameter: float) -> float:
+            # B-mode (see ref. above).
+            d = diameter * 1e4 # microns
+            return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.26)) *
+                    np.exp(-1 * (np.log(d) - 1.0) ** 2 / (2 * 0.26 ** 2)))
+
+        def _Lmode(diameter: float) -> float:
+            # L-mode (see ref. above).
+            d = diameter * 1e4 # microns
+            return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.5)) *
+                    np.exp(-1 * (np.log(d) - 1.4) ** 2 / (2 * 0.5 ** 2)))
+
+        def _Omode(diameter: float) -> float:
+            # O-mode (see ref. above).
+            d = diameter * 1e4 # microns
+            return ( (1 / d) * (0.001 / (np.sqrt(2 * np.pi) * 0.56)) *
+                    np.exp(-1 * (np.log(d) - 4.98) ** 2 / (2 * 0.56 ** 2)))
+
+        def integrand(diameter: float) -> float:
+            return (self.BLO_factors[0] * _Bmode(diameter) +
+                    self.BLO_factors[1] * _Lmode(diameter) +
+                    self.BLO_factors[2] * _Omode(diameter)
+                    ) * volume(diameter) * (1 - mask.exhale_efficiency(diameter))
+
+        return quad(integrand, 0.1, 30)[0]
+
+
 _ExpirationBase.types = {
     'Breathing': Expiration((0.084, 0.009, 0.003, 0.002)),
     'Whispering': Expiration((0.11, 0.014, 0.004, 0.002)),
     'Talking': Expiration((0.236, 0.068, 0.007, 0.011)),
     'Unmodulated Vocalization': Expiration((0.751, 0.139, 0.0139, 0.059)),
     'Superspreading event': Expiration((np.inf, np.inf, np.inf, np.inf)),
+    'Breathing BLO': ExpirationBLO((1., 0., 0.)),
+    'Talking BLO': ExpirationBLO((1., 1., 1.)),
+    'Shouting BLO': ExpirationBLO((5., 5., 5.)),
+    'Singing BLO': ExpirationBLO((5., 5., 5.)),
 }
 
 

From a9cd36df7529213241d0b55d57bf32cc11c99870 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Thu, 27 May 2021 18:35:54 +0200
Subject: [PATCH 06/16] Minor docstring fix in ExpirationBLO

---
 cara/models.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index d491e0da..2bdade1f 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -631,9 +631,9 @@ class ExpirationBLO(_ExpirationBase):
     https://doi.org/10.1016/j.jaerosci.2011.07.009).
     All diameters are in cm.
     """
-    #: factors assigned to resp. the B, L and O modes. Depends on the
-    # kind of expiratory activity (e.g. breathing, speaking, singing,
-    # or shouting).
+    #: factors assigned to resp. the B, L and O modes. They are
+    # charateristics of the kind of expiratory activity (e.g. breathing,
+    # speaking, singing, or shouting).
     BLO_factors: typing.Tuple[float, float, float]
 
     def aerosols(self, mask: _MaskBase):

From 6e49ed4f0213c32a25179593673f863b4ae00bfd Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Fri, 28 May 2021 06:56:13 +0200
Subject: [PATCH 07/16] Improving docstrings in expiration classes

---
 cara/models.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index a7c8d334..a72f1b48 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -573,7 +573,7 @@ class _ExpirationBase:
     types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
 
     def aerosols(self, mask: _MaskBase):
-        # total volume of aerosols expired (cm^3).
+        # total volume of aerosols expired per volume of air (mL/cm^3).
         raise NotImplementedError("Subclass must implement")
 
 
@@ -581,7 +581,10 @@ class _ExpirationBase:
 class Expiration(_ExpirationBase):
     """
     Simple model based on four different sizes of particles emitted,
-    with different ejection factors.
+    with different ejection factors. See Fig. 4 in L. Morawska et al,
+    Size distribution and sites of origin of droplets expelled from the
+    human respiratory tract during expiratory activities,
+    Aerosol Science 40 (2009) pp. 256 - 269.
     """
     ejection_factor: typing.Tuple[float, ...]
     particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4)  # In cm.

From cd3f9057f9cc8d35cd1b0712a8e9afd13e4bce63 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Fri, 28 May 2021 07:08:23 +0200
Subject: [PATCH 08/16] Dealing better with the units in ExpirationBLO.aerosols

---
 cara/models.py | 27 ++++++++++++---------------
 1 file changed, 12 insertions(+), 15 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index c1bee285..496b86e2 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -632,7 +632,7 @@ class ExpirationBLO(_ExpirationBase):
     expired aerosol size distributions, Journal of Aerosol Science,
     vol. 42, no. 12, pp. 839 – 851, 2011,
     https://doi.org/10.1016/j.jaerosci.2011.07.009).
-    All diameters are in cm.
+    Here all diameters (d) are in microns.
     """
     #: factors assigned to resp. the B, L and O modes. They are
     # charateristics of the kind of expiratory activity (e.g. breathing,
@@ -640,33 +640,30 @@ class ExpirationBLO(_ExpirationBase):
     BLO_factors: typing.Tuple[float, float, float]
 
     def aerosols(self, mask: _MaskBase):
+        # result is in mL.cm^-3
+        def volume(d):
+            return (np.pi * d**3) / 6.
 
-        def volume(diameter):
-            return (np.pi * diameter**3) / 6.
-
-        def _Bmode(diameter: float) -> float:
+        def _Bmode(d: float) -> float:
             # B-mode (see ref. above).
-            d = diameter * 1e4 # microns
             return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.26)) *
                     np.exp(-1 * (np.log(d) - 1.0) ** 2 / (2 * 0.26 ** 2)))
 
-        def _Lmode(diameter: float) -> float:
+        def _Lmode(d: float) -> float:
             # L-mode (see ref. above).
-            d = diameter * 1e4 # microns
             return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.5)) *
                     np.exp(-1 * (np.log(d) - 1.4) ** 2 / (2 * 0.5 ** 2)))
 
-        def _Omode(diameter: float) -> float:
+        def _Omode(d: float) -> float:
             # O-mode (see ref. above).
-            d = diameter * 1e4 # microns
             return ( (1 / d) * (0.001 / (np.sqrt(2 * np.pi) * 0.56)) *
                     np.exp(-1 * (np.log(d) - 4.98) ** 2 / (2 * 0.56 ** 2)))
 
-        def integrand(diameter: float) -> float:
-            return (self.BLO_factors[0] * _Bmode(diameter) +
-                    self.BLO_factors[1] * _Lmode(diameter) +
-                    self.BLO_factors[2] * _Omode(diameter)
-                    ) * volume(diameter) * (1 - mask.exhale_efficiency(diameter))
+        def integrand(d: float) -> float:
+            return (self.BLO_factors[0] * _Bmode(d) +
+                    self.BLO_factors[1] * _Lmode(d) +
+                    self.BLO_factors[2] * _Omode(d)
+                    ) * volume(d*1e-4) * (1 - mask.exhale_efficiency(d*1e-4))
 
         return quad(integrand, 0.1, 30)[0]
 

From 8e68bd7a2317b0ccabcbf91298dd639c7c5deda9 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 11:07:53 +0200
Subject: [PATCH 09/16] Removing old Expiration and ref. to ExpirationBLO in
 test_expiration

---
 cara/tests/test_expiration.py | 19 +++++--------------
 1 file changed, 5 insertions(+), 14 deletions(-)

diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py
index 5cfe94b9..ff08798b 100644
--- a/cara/tests/test_expiration.py
+++ b/cara/tests/test_expiration.py
@@ -9,7 +9,7 @@ from cara import models
 
 def test_multiple_wrong_weight_size():
     weights = (1., 2., 3.)
-    e_base = models.Expiration((0.084, 0.009, 0.003, 0.002))
+    e_base = models.Expiration((1., 0., 0.))
     with pytest.raises(
             ValueError,
             match=re.escape("expirations and weigths should contain the"
@@ -19,19 +19,10 @@ def test_multiple_wrong_weight_size():
 
 
 def test_multiple():
-    weights = (1., 2.)
-    e1 = models.Expiration((0.03, 0.02, 0.01, 0.005))
-    e2 = models.Expiration((0.05, 0.04, 0.03, 0.01))
-    e = models.MultipleExpiration([e1, e2], weights)
-    mask = models.Mask.types['No mask']
-    assert e.aerosols(mask) == e1.aerosols(mask)/3. + 2*e2.aerosols(mask)/3.
-
-
-def test_multiple_BLO():
     weights = (1., 1.)
-    e1 = models.ExpirationBLO((1., 0., 0.))
-    e2 = models.ExpirationBLO((4., 5., 5.))
-    e_expected = models.ExpirationBLO((2.5, 2.5, 2.5))
+    e1 = models.Expiration((1., 0., 0.))
+    e2 = models.Expiration((4., 5., 5.))
+    e_expected = models.Expiration((2.5, 2.5, 2.5))
     e = models.MultipleExpiration([e1, e2], weights)
     mask = models.Mask.types['Type I']
-    assert e_expected.aerosols(mask) == e.aerosols(mask)
+    npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask))

From ecf9ce8655a48c6186a8d80e9fc176fd290db756 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 11:08:42 +0200
Subject: [PATCH 10/16] Removing old Expiration class; replaced by previously
 named ExpirationBLO class

---
 cara/models.py | 87 +++++++++++++++++---------------------------------
 1 file changed, 30 insertions(+), 57 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index 284acdd0..f88107c9 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -538,53 +538,6 @@ class _ExpirationBase:
 
 @dataclass(frozen=True)
 class Expiration(_ExpirationBase):
-    """
-    Simple model based on four different sizes of particles emitted,
-    with different ejection factors. See Fig. 4 in L. Morawska et al,
-    Size distribution and sites of origin of droplets expelled from the
-    human respiratory tract during expiratory activities,
-    Aerosol Science 40 (2009) pp. 256 - 269.
-    """
-    ejection_factor: typing.Tuple[float, ...]
-    particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4)  # In cm.
-
-    def aerosols(self, mask: Mask):
-        def volume(diameter):
-            return (4 * np.pi * (diameter/2)**3) / 3
-        total = 0
-        for diameter, factor in zip(self.particle_sizes, self.ejection_factor):
-            contribution = (volume(diameter) * factor *
-                            (1 - mask.exhale_efficiency(diameter*1e4)))
-            total += contribution
-        return total
-
-
-@dataclass(frozen=True)
-class MultipleExpiration(_ExpirationBase):
-    """
-    Represents an expiration of aerosols.
-    Group together different modes of expiration, that represent
-    each the main expiration mode for a certain fraction of time (given by
-    the weights).
-
-    """
-    expirations: typing.Tuple[_ExpirationBase, ...]
-    weights: typing.Tuple[float, ...]
-
-    def __post_init__(self):
-        if len(self.expirations) != len(self.weights):
-            raise ValueError("expirations and weigths should contain the"
-                             "same number of elements")
-
-    def aerosols(self, mask: Mask):
-        return np.array([
-            weight * expiration.aerosols(mask) / sum(self.weights)
-            for weight,expiration in zip(self.weights,self.expirations)
-        ]).sum(axis=0)
-
-
-@dataclass(frozen=True)
-class ExpirationBLO(_ExpirationBase):
     """
     BLO model for the expiration (G. Johnson et al., Modality of human
     expired aerosol size distributions, Journal of Aerosol Science,
@@ -597,7 +550,7 @@ class ExpirationBLO(_ExpirationBase):
     # speaking, singing, or shouting).
     BLO_factors: typing.Tuple[float, float, float]
 
-    def aerosols(self, mask: _MaskBase):
+    def aerosols(self, mask: Mask):
         """ Result is in mL.cm^-3 """
         def volume(d):
             return (np.pi * d**3) / 6.
@@ -626,16 +579,36 @@ class ExpirationBLO(_ExpirationBase):
         return quad(integrand, 0.1, 30)[0]
 
 
+@dataclass(frozen=True)
+class MultipleExpiration(_ExpirationBase):
+    """
+    Represents an expiration of aerosols.
+    Group together different modes of expiration, that represent
+    each the main expiration mode for a certain fraction of time (given by
+    the weights).
+
+    """
+    expirations: typing.Tuple[_ExpirationBase, ...]
+    weights: typing.Tuple[float, ...]
+
+    def __post_init__(self):
+        if len(self.expirations) != len(self.weights):
+            raise ValueError("expirations and weigths should contain the"
+                             "same number of elements")
+
+    def aerosols(self, mask: Mask):
+        return np.array([
+            weight * expiration.aerosols(mask) / sum(self.weights)
+            for weight,expiration in zip(self.weights,self.expirations)
+        ]).sum(axis=0)
+
+
 _ExpirationBase.types = {
-    'Breathing': Expiration((0.084, 0.009, 0.003, 0.002)),
-    'Whispering': Expiration((0.11, 0.014, 0.004, 0.002)),
-    'Talking': Expiration((0.236, 0.068, 0.007, 0.011)),
-    'Unmodulated Vocalization': Expiration((0.751, 0.139, 0.0139, 0.059)),
-    'Superspreading event': Expiration((np.inf, np.inf, np.inf, np.inf)),
-    'Breathing BLO': ExpirationBLO((1., 0., 0.)),
-    'Talking BLO': ExpirationBLO((1., 1., 1.)),
-    'Shouting BLO': ExpirationBLO((5., 5., 5.)),
-    'Singing BLO': ExpirationBLO((5., 5., 5.)),
+    'Breathing': Expiration((1., 0., 0.)),
+    'Talking': Expiration((1., 1., 1.)),
+    'Shouting': Expiration((5., 5., 5.)),
+    'Singing': Expiration((5., 5., 5.)),
+    'Superspreading event': Expiration((np.inf, 0., 0.)),
 }
 
 

From d12e809d8011cb4b4759594df78ce2239ac71784 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 11:29:49 +0200
Subject: [PATCH 11/16] Better import of scipy.integrate

---
 cara/models.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index f88107c9..24e99e4e 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -37,7 +37,7 @@ import typing
 
 import numpy as np
 from scipy.interpolate import interp1d
-from scipy.integrate import quad
+import scipy.integrate
 
 if not typing.TYPE_CHECKING:
     from memoization import cached
@@ -576,7 +576,7 @@ class Expiration(_ExpirationBase):
                     self.BLO_factors[2] * _Omode(d)
                     ) * volume(d*1e-4) * (1 - mask.exhale_efficiency(d*1e-4))
 
-        return quad(integrand, 0.1, 30)[0]
+        return scipy.integrate.quad(integrand, 0.1, 30)[0]
 
 
 @dataclass(frozen=True)

From 612a4c8e08102546bcbe7fd60f22a22c256f67ea Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 16:38:00 +0200
Subject: [PATCH 12/16] New expiration test to make sure expiration rates are
 correct

---
 cara/tests/test_expiration.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py
index ff08798b..1a6fdbb7 100644
--- a/cara/tests/test_expiration.py
+++ b/cara/tests/test_expiration.py
@@ -26,3 +26,18 @@ def test_multiple():
     e = models.MultipleExpiration([e1, e2], weights)
     mask = models.Mask.types['Type I']
     npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask))
+
+
+# expected values obtained from separate model
+@pytest.mark.parametrize(
+    "BLO_weights, expected_aerosols",
+    [
+        [(1.,0.,0.), 1.38924e-12],
+        [(1.,1.,1.), 1.07129e-10],
+        [(1.,5.,5.), 5.30088e-10],
+    ],
+)
+def test_expiration_aerosols(BLO_weights, expected_aerosols):
+    mask = models.Mask.types['No mask']
+    e = models.Expiration(BLO_weights)
+    npt.assert_allclose(e.aerosols(mask), expected_aerosols, rtol=1e-4)

From 3a390435229756cd81bf84ee139d5fbddcb79974 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 16:41:41 +0200
Subject: [PATCH 13/16] More accurate B, L and O modes; correcting units

---
 cara/models.py | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index 24e99e4e..a5c63b2e 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -557,26 +557,27 @@ class Expiration(_ExpirationBase):
 
         def _Bmode(d: float) -> float:
             # B-mode (see ref. above).
-            return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.26)) *
-                    np.exp(-1 * (np.log(d) - 1.0) ** 2 / (2 * 0.26 ** 2)))
+            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:
             # L-mode (see ref. above).
-            return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.5)) *
-                    np.exp(-1 * (np.log(d) - 1.4) ** 2 / (2 * 0.5 ** 2)))
+            return ( (1 / d) * (1.0 / (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:
             # O-mode (see ref. above).
-            return ( (1 / d) * (0.001 / (np.sqrt(2 * np.pi) * 0.56)) *
-                    np.exp(-1 * (np.log(d) - 4.98) ** 2 / (2 * 0.56 ** 2)))
+            return ( (1 / d) * (0.0010008 / (np.sqrt(2 * np.pi) * 0.585005)) *
+                    np.exp(-1 * (np.log(d) - 4.97673) ** 2 / (2 * 0.585005 ** 2)))
 
         def integrand(d: float) -> float:
             return (self.BLO_factors[0] * _Bmode(d) +
                     self.BLO_factors[1] * _Lmode(d) +
                     self.BLO_factors[2] * _Omode(d)
-                    ) * volume(d*1e-4) * (1 - mask.exhale_efficiency(d*1e-4))
+                    ) * volume(d) * (1 - mask.exhale_efficiency(d))
 
-        return scipy.integrate.quad(integrand, 0.1, 30)[0]
+        # final result converted from microns^3/cm3 to mL/cm^3
+        return scipy.integrate.quad(integrand, 0.1, 30.)[0]*1e-12
 
 
 @dataclass(frozen=True)
@@ -606,8 +607,8 @@ class MultipleExpiration(_ExpirationBase):
 _ExpirationBase.types = {
     'Breathing': Expiration((1., 0., 0.)),
     'Talking': Expiration((1., 1., 1.)),
-    'Shouting': Expiration((5., 5., 5.)),
-    'Singing': Expiration((5., 5., 5.)),
+    'Shouting': Expiration((1., 5., 5.)),
+    'Singing': Expiration((1., 5., 5.)),
     'Superspreading event': Expiration((np.inf, 0., 0.)),
 }
 

From 62a4880a5b92227b1b4887337f0c63e523d7ca4a Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 17:12:18 +0200
Subject: [PATCH 14/16] Avoiding vectorisation of factor_exhale, and changing
 test_infected_population accordingly (incompatibility with scipy quad
 integration)

---
 cara/models.py                         | 2 +-
 cara/tests/test_infected_population.py | 8 ++------
 2 files changed, 3 insertions(+), 7 deletions(-)

diff --git a/cara/models.py b/cara/models.py
index a5c63b2e..481074cd 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -478,7 +478,7 @@ class Mask:
     η_inhale: _VectorisedFloat
 
     #: Global factor applied to filtration efficiency of masks when exhaling.
-    factor_exhale: _VectorisedFloat = 1.
+    factor_exhale: float = 1.
 
     #: Pre-populated examples of Masks.
     types: typing.ClassVar[typing.Dict[str, "Mask"]]
diff --git a/cara/tests/test_infected_population.py b/cara/tests/test_infected_population.py
index 8ef94729..7acef7f7 100644
--- a/cara/tests/test_infected_population.py
+++ b/cara/tests/test_infected_population.py
@@ -8,7 +8,6 @@ import cara.models
     "override_params", [
         {'viral_load_in_sputum': np.array([5e8, 1e9])},
         {'quantum_infectious_dose': np.array([50, 20])},
-        {'factor_exhale': np.array([0.92, 0.95])},
         {'exhalation_rate': np.array([0.75, 0.81])},
     ]
 )
@@ -16,7 +15,6 @@ def test_infected_population_vectorisation(override_params):
     defaults = {
         'viral_load_in_sputum': 1e9,
         'quantum_infectious_dose': 50,
-        'factor_exhale': 0.95,
         'exhalation_rate': 0.75,
     }
     defaults.update(override_params)
@@ -26,7 +24,7 @@ def test_infected_population_vectorisation(override_params):
             number=1,
             presence=office_hours,
             mask=cara.models.Mask(
-                factor_exhale=defaults['factor_exhale'],
+                factor_exhale=0.95,
                 η_inhale=0.3,
             ),
             activity=cara.models.Activity(
@@ -37,9 +35,7 @@ def test_infected_population_vectorisation(override_params):
                 viral_load_in_sputum=defaults['viral_load_in_sputum'],
                 quantum_infectious_dose=defaults['quantum_infectious_dose'],
             ),
-            expiration=cara.models.Expiration(
-                ejection_factor=(0.084, 0.009, 0.003, 0.002),
-        ),
+            expiration=cara.models.Expiration((1., 0., 0.)),
     )
     emission_rate = infected.emission_rate(10)
     assert isinstance(emission_rate, np.ndarray)

From 4a211aa7fc453aa6daa8a03c6b8cb106e9bb1ee3 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 23:09:22 +0200
Subject: [PATCH 15/16] Fixing tests with the change of Expiration, in
 particular: removing useless tests on aerosols in test_known_quantities;
 removing reference to Unmodulated vocalization expiration type; avoiding the
 use of factor_exhale as array; changing some integration test values on r0;
 replacing some integration tests on r0 by integration tests on
 quanta_exposure with a simpler model; changing values in some concentration
 tests

---
 cara/tests/conftest.py                        | 12 ++-
 cara/tests/models/test_concentration_model.py |  8 +-
 cara/tests/test_expiration.py                 |  2 +-
 cara/tests/test_known_quantities.py           | 77 ++++++++-----------
 cara/tests/test_monte_carlo.py                |  2 +-
 5 files changed, 39 insertions(+), 62 deletions(-)

diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py
index 24983cf1..44b25423 100644
--- a/cara/tests/conftest.py
+++ b/cara/tests/conftest.py
@@ -7,11 +7,9 @@ import pytest
 def baseline_model():
     model = models.ConcentrationModel(
         room=models.Room(volume=75),
-        ventilation=models.SlidingWindow(
-            active=models.PeriodicInterval(period=120, duration=120),
-            inside_temp=models.PiecewiseConstant((0,24),(293,)),
-            outside_temp=models.PiecewiseConstant((0,24),(283,)),
-            window_height=1.6, opening_length=0.6,
+        ventilation=models.AirChange(
+            active=models.SpecificInterval(((0,24),)),
+            air_exch=30.,
         ),
         infected=models.InfectedPopulation(
             number=1,
@@ -19,7 +17,7 @@ def baseline_model():
             presence=models.SpecificInterval(((0, 4), (5, 8))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -30,7 +28,7 @@ def baseline_exposure_model(baseline_model):
     return models.ExposureModel(
         baseline_model,
         exposed=models.Population(
-            number=10,
+            number=1000,
             presence=baseline_model.infected.presence,
             activity=baseline_model.infected.activity,
             mask=baseline_model.infected.mask,
diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py
index f7230430..26707236 100644
--- a/cara/tests/models/test_concentration_model.py
+++ b/cara/tests/models/test_concentration_model.py
@@ -14,7 +14,6 @@ from cara import models
         {'air_change': np.array([100, 120])},
         {'viral_load_in_sputum': np.array([5e8, 1e9])},
         {'quantum_infectious_dose': np.array([50, 20])},
-        {'factor_exhale': np.array([0.92, 0.95])},
     ]
 )
 def test_concentration_model_vectorisation(override_params):
@@ -24,7 +23,6 @@ def test_concentration_model_vectorisation(override_params):
         'air_change': 100,
         'viral_load_in_sputum': 1e9,
         'quantum_infectious_dose': 50,
-        'factor_exhale': 0.95,
     }
     defaults.update(override_params)
 
@@ -36,7 +34,7 @@ def test_concentration_model_vectorisation(override_params):
             number=1,
             presence=always,
             mask=models.Mask(
-                factor_exhale=defaults['factor_exhale'],
+                factor_exhale=0.95,
                 η_inhale=0.3,
             ),
             activity=models.Activity(
@@ -47,9 +45,7 @@ def test_concentration_model_vectorisation(override_params):
                 viral_load_in_sputum=defaults['viral_load_in_sputum'],
                 quantum_infectious_dose=defaults['quantum_infectious_dose'],
             ),
-            expiration=models.Expiration(
-                ejection_factor=(0.084, 0.009, 0.003, 0.002),
-            ),
+            expiration=models.Expiration((1., 0., 0.)),
         )
     )
     concentrations = c_model.concentration(10)
diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py
index 1a6fdbb7..434caf8c 100644
--- a/cara/tests/test_expiration.py
+++ b/cara/tests/test_expiration.py
@@ -28,7 +28,7 @@ def test_multiple():
     npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask))
 
 
-# expected values obtained from separate model
+# expected values obtained from another code
 @pytest.mark.parametrize(
     "BLO_weights, expected_aerosols",
     [
diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py
index ea8a1ad0..a3b1e5db 100644
--- a/cara/tests/test_known_quantities.py
+++ b/cara/tests/test_known_quantities.py
@@ -6,21 +6,12 @@ import cara.models as models
 import cara.data as data
 
 
-def test_no_mask_aerosols(baseline_model):
-    exp = models.Expiration.types['Unmodulated Vocalization']
-    npt.assert_allclose(
-        exp.aerosols(models.Mask.types['No mask']),
-        6.077541e-12,
-        rtol=1e-5,
-    )
-
-
-def test_no_mask_emission_rate(baseline_model):
-    rate = 151.938514
+def test_no_mask_superspeading_emission_rate(baseline_model):
+    expected_rate = 970.
     npt.assert_allclose(
         [baseline_model.infected.emission_rate(t) for t in [0, 1, 4, 4.5, 5, 8, 9]],
-        [0, rate, rate, 0, 0, rate, 0],
-        rtol=1e-5
+        [0, expected_rate, expected_rate, 0, 0, expected_rate, 0],
+        rtol=1e-12
     )
 
 
@@ -48,23 +39,24 @@ def baseline_periodic_hepa():
 
 
 def test_concentrations(baseline_model):
+    # expected concentrations were computed analytically
     ts = [0, 4, 5, 7, 10]
     concentrations = [baseline_model.concentration(t) for t in ts]
     npt.assert_allclose(
         concentrations,
-        [0.000000e+00, 2.619538e-01, 1.146999e-04, 2.619537e-01, 5.022289e-08],
-        rtol=1e-5
+        [0.000000e+00, 0.4189594, 1.6422648e-14, 0.4189594, 6.4374587e-28],
+        rtol=1e-6
     )
 
 
 def test_smooth_concentrations(baseline_model):
     # We don't care about the actual concentrations in this test, but rather
     # that the curve itself is smooth.
-    dx = 0.1
-    dy_limit = dx * 2  # Anything more than this is a bit steep.
+    dx = 0.002
+    dy_limit = 0.2  # Anything more than this (in relative) is a bit steep.
     ts = np.arange(0, 10, dx)
     concentrations = [baseline_model.concentration(t) for t in ts]
-    assert np.abs(np.diff(concentrations)).max() < dy_limit
+    assert np.abs(np.diff(concentrations)).max()/np.mean(concentrations) < dy_limit
 
 
 def build_model(interval_duration):
@@ -80,7 +72,7 @@ def build_model(interval_duration):
             presence=models.SpecificInterval(((0, 4), (5, 8))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -98,8 +90,8 @@ def test_concentrations_startup(baseline_model):
 def test_r0(baseline_exposure_model):
     # expected r0 was computed with a trapezoidal integration, using
     # a mesh of 100'000 pts per exposed presence interval.
-    p = baseline_exposure_model.infection_probability()
-    npt.assert_allclose(p, 89.001748)
+    r0 = baseline_exposure_model.reproduction_number()
+    npt.assert_allclose(r0, 973.535888)
 
 
 def test_periodic_window(baseline_periodic_window, baseline_room):
@@ -183,15 +175,6 @@ def test_multiple_ventilation_HEPA_HVAC_AirChange(volume, expected_value):
                         expected_value,rtol=1e-5)
 
 
-def test_expiration_aerosols():
-    mask = models.Mask.types['Type I']
-    exp1 = models.Expiration((0.751, 0.139, 0.0139, 0.059),
-                            particle_sizes = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4))
-    exp2 = models.Expiration((0.059, 0.0139, 0.751, 0.139),
-                            particle_sizes = (5.5e-4, 3.5e-4, 0.8e-4, 1.8e-4))
-    npt.assert_allclose(exp1.aerosols(mask), exp2.aerosols(mask), rtol=1e-5)
-
-
 @pytest.mark.parametrize(
     "time, expected_value",
     [
@@ -239,7 +222,7 @@ def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),),
             presence=models.SpecificInterval(intervals_presence_infected),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -260,7 +243,7 @@ def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)):
             presence=models.SpecificInterval(((0, 4), (5, 7.5))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -287,7 +270,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5
             presence=models.SpecificInterval(((0, 4), (5, 7.5))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -370,16 +353,16 @@ def build_exposure_model(concentration_model):
     )
 
 
-# expected r0 were computed with a trapezoidal integration, using
+# expected quanta were computed with a trapezoidal integration, using
 # a mesh of 100'000 pts per exposed presence interval.
 @pytest.mark.parametrize(
-    "month, expected_r0",
+    "month, expected_quanta",
     [
-        ['Jan', 86.258533],
-        ['Jun', 99.972103],
+        ['Jan', 10.136783],
+        ['Jun', 41.800377],
     ],
 )
-def test_r0_hourly_dep(month,expected_r0):
+def test_quanta_hourly_dep(month,expected_quanta):
     m = build_exposure_model(
         build_hourly_dependent_model(
             month,
@@ -387,19 +370,19 @@ def test_r0_hourly_dep(month,expected_r0):
             intervals_presence_infected=((8, 12), (13, 17))
         )
     )
-    p = m.infection_probability()
-    npt.assert_allclose(p, expected_r0)
+    quanta = m.quanta_exposure()
+    npt.assert_allclose(quanta, expected_quanta)
 
-# expected r0 were computed with a trapezoidal integration, using
+# expected quanta were computed with a trapezoidal integration, using
 # a mesh of 100'000 pts per exposed presence interval.
 @pytest.mark.parametrize(
-    "month, expected_r0",
+    "month, expected_quanta",
     [
-        ['Jan', 86.422736],
-        ['Jun', 99.982323],
+        ['Jan', 10.19818],
+        ['Jun', 44.130683],
     ],
 )
-def test_r0_hourly_dep_refined(month,expected_r0):
+def test_quanta_hourly_dep_refined(month,expected_quanta):
     m = build_exposure_model(
         build_hourly_dependent_model(
             month,
@@ -408,5 +391,5 @@ def test_r0_hourly_dep_refined(month,expected_r0):
             temperatures=data.GenevaTemperatures,
         )
     )
-    p = m.infection_probability()
-    npt.assert_allclose(p, expected_r0)
+    quanta = m.quanta_exposure()
+    npt.assert_allclose(quanta, expected_quanta)
diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py
index f3d9fb08..c3bfc2fe 100644
--- a/cara/tests/test_monte_carlo.py
+++ b/cara/tests/test_monte_carlo.py
@@ -53,7 +53,7 @@ def baseline_mc_model() -> cara.monte_carlo.ConcentrationModel:
             presence=cara.models.SpecificInterval(((0, 4), (5, 8))),
             mask=cara.models.Mask.types['No mask'],
             activity=cara.models.Activity.types['Light activity'],
-            expiration=cara.models.Expiration.types['Unmodulated Vocalization'],
+            expiration=cara.models.Expiration.types['Breathing'],
         ),
     )
     return mc_model

From 390d2299445296465075b7981c23dd34ed5955f4 Mon Sep 17 00:00:00 2001
From: Nicolas Mounet <nicolas.mounet@cern.ch>
Date: Mon, 31 May 2021 23:23:42 +0200
Subject: [PATCH 16/16] Caching aerosols computation

---
 cara/models.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/cara/models.py b/cara/models.py
index 481074cd..cf2ffb8f 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -550,6 +550,7 @@ class Expiration(_ExpirationBase):
     # speaking, singing, or shouting).
     BLO_factors: typing.Tuple[float, float, float]
 
+    @cached()
     def aerosols(self, mask: Mask):
         """ Result is in mL.cm^-3 """
         def volume(d):