Merge branch 'master' into paper/deposition_fraction

cara/apps/calculator/__init__.py

@@ -33,7 +33,7 @@ from .user import AuthenticatedUser, AnonymousUser
 # calculator version. If the calculator needs to make breaking changes (e.g. change
 # form attributes) then it can also increase its MAJOR version without needing to
 # increase the overall CARA version (found at ``cara.__version__``).
-__version__ = "3.0.0"
+__version__ = "3.0.1"
 
 
 class BaseRequestHandler(RequestHandler):

cara/apps/calculator/templates/calculator.form.html.j2

@@ -65,7 +65,7 @@ v{{ calculator_version }} <span style="float:right; font-weight:bold">Please sen
           <option value="SARS_CoV_2">SARS-CoV-2 (nominal strain)</option>
           <option value="SARS_CoV_2_B117">SARS-CoV-2 (Alpha VOC)</option>
           <option value="SARS_CoV_2_P1">SARS-CoV-2 (Gamma VOC)</option>
-          <option value="SARS_CoV_2_B16172">SARS-CoV-2 (Delta VOC)</option>
+          <option selected value="SARS_CoV_2_B16172">SARS-CoV-2 (Delta VOC)</option>
         </select>
     </div>
 
@@ -344,13 +344,13 @@ v{{ calculator_version }} <span style="float:right; font-weight:bold">Please sen
   <b>Quick Guide:</b><br>
   This tool simulates the long range airborne spread SARS-CoV-2 virus in a finite volume and estimates the risk of COVID-19 infection. It is based on current scientific data and can be used to compare the effectiveness of different mitigation measures.<br>
   <b>Virus data:</b> <br>
-  SARS-CoV-2 covers typical strains of the virus and three variants of concern (VOC):<br>
+  SARS-CoV-2 covers the original "wild type" strain of the virus and three variants of concern (VOC):<br>
   <ul>
     <li>Alpha (also known as B.1.1.7, first identified in UK, Dec 2020),</li>
     <li>Gamma (also known as P.1, first identified in Brazil/Japan, Jan 2021).</li>
     <li>Delta (also known as B.1.617.2, first identified in India, Oct 2020).</li>
   </ul>
-  Choose variant according to local area prevalence, e.g. for <a href="https://www.covid19.admin.ch/fr/epidemiologic/virus-variants?detGeo=GE">Geneva</a>
+  Modify the default as necessary, according to local area prevalence e.g. for <a href="https://www.covid19.admin.ch/fr/epidemiologic/virus-variants?detGeo=GE">Geneva</a>
   or <a href="https://www.santepubliquefrance.fr/dossiers/coronavirus-covid-19/covid-19-cartographie-des-variants-en-france-donnees-par-region-et-par-departement">Ain (France)</a>.<br>
   <b>Ventilation data:</b> <br>
   <ul>

cara/apps/calculator/templates/userguide.html.j2

@@ -56,12 +56,12 @@ Changing this setting alters the properties of the virus which are used for the
 This has a significant effect on the probability of infection. 
 The choices are:</p>
 <ul>
-<li><code>SARS-CoV-2 (nominal strain)</code>, covering typical strains and varaints which are not of concern from an epidemiologic point of view of the virus;</li>
+<li><code>SARS-CoV-2 (nominal strain)</code>, covering typical strains and variants which are not of concern from an epidemiologic point of view of the virus;</li>
 <li><code>SARS-CoV-2 (Alpha VOC)</code>, first identified in the UK at the end of 2020 which is found to be approximately 1.5x more transmissible compared to the non-VOCs; </li>
 <li><code>SARS-CoV-2 (Gamma VOC)</code>, first identified in Brazil in January 2021 which is found to be approximately 2.2x more transmissible compared to the non-VOCs.</li>
 <li><code>SARS-CoV-2 (Delta VOC)</code>, first identified in India towards the end of 2020 which is found to be approximately 60% more transmissible compared to the ALPHA VOC.</li>
 </ul>
-<p>The user can base their choice according to the prevalence of the different variants in the local area. Access to this information can be found here:</p>
+<p>The user can modify the selected variant from the default, according to the prevalence of the different variants in the local area. Access to this information can be found here:</p>
 <ul>
 <li>Geneva: <a href="https://www.covid19.admin.ch/fr/epidemiologic/virus-variants?detGeo=GE">https://www.covid19.admin.ch/fr/epidemiologic/virus-variants?detGeo=GE</a></li>
 <li>Ain (France): <a href="https://www.santepubliquefrance.fr/dossiers/coronavirus-covid-19/covid-19-cartographie-des-variants-en-france-donnees-par-region-et-par-departement">https://www.santepubliquefrance.fr/dossiers/coronavirus-covid-19/covid-19-cartographie-des-variants-en-france-donnees-par-region-et-par-departement</a></li>

cara/models.py

@@ -710,11 +710,11 @@ class InfectedPopulation(Population):
         elif np.isinf(aerosols):
             ER = 970 * self.virus.infectious_dose
 
-        return ER
+        return ER * self.number
 
-    def individual_emission_rate(self, time) -> _VectorisedFloat:
+    def emission_rate(self, time) -> _VectorisedFloat:
         """
-        The emission rate of a single individual in the population.
+        The emission rate of the population.
 
         """
         # Note: The original model avoids time dependence on the emission rate
@@ -730,13 +730,6 @@ class InfectedPopulation(Population):
 
         return self.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
 
-    def emission_rate(self, time) -> _VectorisedFloat:
-        """
-        The emission rate of the entire population.
-
-        """
-        return self.individual_emission_rate(time) * self.number
-
 
 @dataclass(frozen=True)
 class ConcentrationModel:
@@ -762,16 +755,22 @@ class ConcentrationModel:
         )
 
     @method_cache
-    def _concentration_limit(self, time: float) -> _VectorisedFloat:
+    def _normed_concentration_limit(self, time: float) -> _VectorisedFloat:
         """
         Provides a constant that represents the theoretical asymptotic 
         value reached by the concentration when time goes to infinity,
         if all parameters were to stay time-independent.
+        This is normalized by the emission rate, the latter acting as a
+        multiplicative constant factor for the concentration model that
+        can be put back in front of the concentration after the time
+        dependence has been solved for.
         """
+        if not self.infected.person_present(time):
+            return 0.
         V = self.room.volume
         IVRR = self.infectious_virus_removal_rate(time)
 
-        return (self.infected.emission_rate(time)) / (IVRR * V)
+        return 1. / (IVRR * V)
 
     @method_cache
     def state_change_times(self) -> typing.List[float]:
@@ -785,6 +784,14 @@ class ConcentrationModel:
         state_change_times.update(self.ventilation.transition_times())
         return sorted(state_change_times)
 
+    @method_cache
+    def _first_presence_time(self) -> float:
+        """
+        First presence time. Before that, the concentration is zero.
+
+        """
+        return self.infected.presence.boundaries()[0][0]
+
     def last_state_change(self, time: float) -> float:
         """
         Find the most recent/previous state change.
@@ -816,15 +823,16 @@ class ConcentrationModel:
         )
 
     @method_cache
-    def _concentration_cached(self, time: float) -> _VectorisedFloat:
+    def _normed_concentration_cached(self, time: float) -> _VectorisedFloat:
-        # A cached version of the concentration method. Use this method if you
+        # A cached version of the _normed_concentration method. Use this
-        # expect that there may be multiple concentration calculations for the
+        # method if you expect that there may be multiple concentration
-        # same time (e.g. at state change times).
+        # calculations for the same time (e.g. at state change times).
-        return self.concentration(time)
+        return self._normed_concentration(time)
 
-    def concentration(self, time: float) -> _VectorisedFloat:
+    def _normed_concentration(self, time: float) -> _VectorisedFloat:
         """
-        Virus exposure concentration, as a function of time.
+        Virus exposure concentration, as a function of time, and
+        normalized by the emission rate.
         The formulas used here assume that all parameters (ventilation,
         emission rate) are constant between two state changes - only
         the value of these parameters at the next state change, are used.
@@ -832,27 +840,42 @@ class ConcentrationModel:
         Note that time is not vectorised. You can only pass a single float
         to this method.
         """
-        if time == 0:
+        # The model always starts at t=0, but we avoid running concentration calculations
+        # before the first presence as an optimisation.
+        if time <= self._first_presence_time():
             return 0.0
         next_state_change_time = self._next_state_change(time)
         IVRR = self.infectious_virus_removal_rate(next_state_change_time)
-        concentration_limit = self._concentration_limit(next_state_change_time)
+        conc_limit = self._normed_concentration_limit(next_state_change_time)
 
         t_last_state_change = self.last_state_change(time)
-        concentration_at_last_state_change = self._concentration_cached(t_last_state_change)
+        conc_at_last_state_change = self._normed_concentration_cached(t_last_state_change)
 
         delta_time = time - t_last_state_change
         fac = np.exp(-IVRR * delta_time)
-        return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac
+        return conc_limit * (1 - fac) + conc_at_last_state_change * fac
+
+    def concentration(self, time: float) -> _VectorisedFloat:
+        """
+        Virus exposure concentration, as a function of time.
+
+        Note that time is not vectorised. You can only pass a single float
+        to this method.
+        """
+        return (self._normed_concentration(time) * 
+                self.infected.emission_rate_when_present())
 
     @method_cache
-    def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat:
+    def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat:
         """
-        Get the integrated concentration dose between the times start and stop.
+        Get the integrated concentration dose between the times start and stop,
+        normalized by the emission rate.
         """
+        if stop <= self._first_presence_time():
+            return 0.0
         state_change_times = self.state_change_times()
         req_start, req_stop = start, stop
-        total_concentration = 0.
+        total_normed_concentration = 0.
         for interval_start, interval_stop in zip(state_change_times[:-1], state_change_times[1:]):
             if req_start > interval_stop or req_stop < interval_start:
                 continue
@@ -860,17 +883,24 @@ class ConcentrationModel:
             start = max([interval_start, req_start])
             stop = min([interval_stop, req_stop])
 
-            conc_start = self._concentration_cached(start)
+            conc_start = self._normed_concentration_cached(start)
 
             next_conc_state = self._next_state_change(stop)
-            conc_limit = self._concentration_limit(next_conc_state)
+            conc_limit = self._normed_concentration_limit(next_conc_state)
             IVRR = self.infectious_virus_removal_rate(next_conc_state)
             delta_time = stop - start
-            total_concentration += (
+            total_normed_concentration += (
                 conc_limit * delta_time +
                 (conc_limit - conc_start) * (np.exp(-IVRR*delta_time)-1) / IVRR
             )
-        return total_concentration
+        return total_normed_concentration
+
+    def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat:
+        """
+        Get the integrated concentration dose between the times start and stop.
+        """
+        return (self.normed_integrated_concentration(start, stop) *
+                self.infected.emission_rate_when_present())
 
 
 @dataclass(frozen=True)
@@ -887,14 +917,22 @@ class ExposureModel:
     #: The fraction of viruses actually deposited in the respiratory tract
     fraction_deposited: _VectorisedFloat = 0.6
 
-    def exposure(self) -> _VectorisedFloat:
+    def _normed_exposure(self) -> _VectorisedFloat:
-        """The number of virus per meter^3."""
+        """
-        exposure = 0.0
+        The number of virus per meter^3, normalized by the emission rate
+        of the infected population.
+        """
+        normed_exposure = 0.0
 
         for start, stop in self.exposed.presence.boundaries():
-            exposure += self.concentration_model.integrated_concentration(start, stop)
+            normed_exposure += self.concentration_model.normed_integrated_concentration(start, stop)
 
-        return exposure * self.repeats
+        return normed_exposure * self.repeats
+
+    def exposure(self) -> _VectorisedFloat:
+        """The number of virus per meter^3."""
+        return (self._normed_exposure() *
+                self.concentration_model.infected.emission_rate_when_present())
 
     def infection_probability(self) -> _VectorisedFloat:
         exposure = self.exposure()

cara/tests/models/test_concentration_model.py

@@ -121,6 +121,10 @@ def test_next_state_change_time_out_of_range(simple_conc_model: models.Concentra
         simple_conc_model._next_state_change(3.1)
 
 
+def test_first_presence_time(simple_conc_model):
+    assert simple_conc_model._first_presence_time() == 0.5
+
+
 def test_integrated_concentration(simple_conc_model):
     c1 = simple_conc_model.integrated_concentration(0, 2)
     c2 = simple_conc_model.integrated_concentration(0, 1)

cara/tests/models/test_exposure_model.py

@@ -11,20 +11,20 @@ from cara.dataclass_utils import replace
 
 
 @dataclass(frozen=True)
-class KnownConcentrations(models.ConcentrationModel):
+class KnownNormedconcentration(models.ConcentrationModel):
     """
     A ConcentrationModel which is based on pre-known exposure concentrations and
     which therefore doesn't need other components. Useful for testing.
 
     """
-    concentration_function: typing.Callable
+    normed_concentration_function: typing.Callable
 
     def infectious_virus_removal_rate(self, time: float) -> models._VectorisedFloat:
         # very large decay constant -> same as constant concentration
         return 1.e50
 
-    def _concentration_limit(self, time: float) -> models._VectorisedFloat:
+    def _normed_concentration_limit(self, time: float) -> models._VectorisedFloat:
-        return self.concentration_function(time)
+        return self.normed_concentration_function(time)
 
     def state_change_times(self):
         return [0., 24.]
@@ -32,8 +32,8 @@ class KnownConcentrations(models.ConcentrationModel):
     def _next_state_change(self, time: float):
         return 24.
 
-    def concentration(self, time: float) -> models._VectorisedFloat:  # noqa
+    def _normed_concentration(self, time: float) -> models._VectorisedFloat:  # noqa
-        return self.concentration_function(time)
+        return self.normed_concentration_function(time)
 
 
 halftime = models.PeriodicInterval(120, 60)
@@ -67,7 +67,9 @@ def known_concentrations(func):
         virus=models.Virus.types['SARS_CoV_2_B117'],
         expiration=models.Expiration.types['Talking']
     )
-    return KnownConcentrations(dummy_room, dummy_ventilation, dummy_infected_population, func)
+    normed_func = lambda x: func(x) / dummy_infected_population.emission_rate_when_present()
+    return KnownNormedconcentration(dummy_room, dummy_ventilation,
+                                dummy_infected_population, normed_func)
 
 
 @pytest.mark.parametrize(

cara/tests/test_monte_carlo_full_models.py

@@ -238,7 +238,7 @@ def skagit_chorale_mc():
         ["shared_office_mc", 10.7, 0.32, 57.24, 654],
         ["classroom_mc",     36.1, 6.85, 780.0, 28464],
         ["ski_cabin_mc",     16.3, 0.49, 35.94, 7404],
-        ["gym_mc",           2.25, 0.63, 0.7842, 984],
+        ["gym_mc",           2.25, 0.63, 0.7842, 1968],
         ["waiting_room_mc",  9.72, 1.36, 34.26, 3534],
         ["skagit_chorale_mc",29.9, 17.9, 190.0, 141400],
     ]