Normalization of the computation of the concentration

cara/models.py

@@ -701,11 +701,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
@@ -721,13 +721,6 @@ class InfectedPopulation(Population):
 
         return self.emission_rate_when_present()
 
-    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:
@@ -753,16 +746,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]:
@@ -815,15 +814,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.
@@ -837,25 +837,36 @@ class ConcentrationModel:
             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
@@ -863,17 +874,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)
@@ -890,14 +908,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_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],
     ]