Merge https://gitlab.cern.ch/caimira/caimira into feature/multiple_concentrationmodels

Author: mheie

caimira/src/caimira/calculator/models/models.py

@@ -1116,18 +1116,6 @@ def last_state_change(self, time: float) -> float:
         t_index = max([t_index - 1, 0])
         return times[t_index]
 
-    def _next_state_change(self, time: float) -> float:
-        """
-        Find the nearest future state change.
-        """
-        for change_time in self.state_change_times():
-            if change_time >= time:
-                return change_time
-        raise ValueError(
-            f"The requested time ({time}) is greater than last available "
-            f"state change time ({change_time})"
-        )
-
     @method_cache
     def _normed_concentration_cached(self, time: float) -> _VectorisedFloat:
         """
@@ -1153,8 +1141,7 @@ def _normed_concentration(self, time: float) -> _VectorisedFloat:
         if time <= self._first_presence_time():
             return self.min_background_concentration()/self.normalization_factor()
 
-        next_state_change_time = self._next_state_change(time)
-        RR = self.removal_rate(next_state_change_time)
+        RR = self.removal_rate(time)
 
         t_last_state_change = self.last_state_change(time)
         conc_at_last_state_change = self._normed_concentration_cached(t_last_state_change)
@@ -1165,12 +1152,12 @@ def _normed_concentration(self, time: float) -> _VectorisedFloat:
             curr_conc_state = np.empty(RR.shape, dtype=np.float64)
             curr_conc_state[RR == 0.] = delta_time * self.population.people_present(time) / (
                 self.room.volume[RR == 0.] if isinstance(self.room.volume,np.ndarray) else self.room.volume)
-            curr_conc_state[RR != 0.] = self._normed_concentration_limit(next_state_change_time)[RR != 0.] * (1 - fac[RR != 0.])
+            curr_conc_state[RR != 0.] = self._normed_concentration_limit(time)[RR != 0.] * (1 - fac[RR != 0.])
         else:
             if RR == 0.:
                 curr_conc_state = delta_time * self.population.people_present(time) / self.room.volume
             else:
-                curr_conc_state = self._normed_concentration_limit(next_state_change_time) * (1 - fac)
+                curr_conc_state = self._normed_concentration_limit(time) * (1 - fac)
 
         return curr_conc_state + conc_at_last_state_change * fac
 
@@ -1193,21 +1180,21 @@ def normed_integrated_concentration(self, start: float, stop: float) -> _Vectori
         """
         if stop <= self._first_presence_time():
             return (stop - start)*self.min_background_concentration()/self.normalization_factor()
-        state_change_times = self.state_change_times()
+        change_times = self.state_change_times()
+        if stop > change_times[-1]:
+            change_times.append(stop)
         req_start, req_stop = start, stop
         total_normed_concentration = 0.
-        for interval_start, interval_stop in zip(state_change_times[:-1], state_change_times[1:]):
+        for interval_start, interval_stop in zip(change_times[:-1], change_times[1:]):
             if req_start > interval_stop or req_stop < interval_start:
                 continue
             # Clip the current interval to the requested range.
             start = max([interval_start, req_start])
             stop = min([interval_stop, req_stop])
 
             conc_start = self._normed_concentration_cached(start)
-
-            next_conc_state = self._next_state_change(stop)
-            conc_limit = self._normed_concentration_limit(next_conc_state)
-            RR = self.removal_rate(next_conc_state)
+            conc_limit = self._normed_concentration_limit(stop)
+            RR = self.removal_rate(stop)
             delta_time = stop - start
             total_normed_concentration += (
                 conc_limit * delta_time +

caimira/tests/models/test_co2_concentration_model.py

@@ -20,6 +20,19 @@ def simple_co2_conc_model(data_registry):
         ),
     )
 
+@pytest.fixture
+def simple_co2_conc_model_extended_presence(data_registry):
+    return models.CO2ConcentrationModel(
+        data_registry=data_registry,
+        room=models.Room(200, models.PiecewiseConstant((0., 24.), (293,))),
+        ventilation=models.AirChange(models.PeriodicInterval(period=120, duration=120), 0.25),
+        CO2_emitters=models.SimplePopulation(
+            number=5,
+            presence=models.SpecificInterval((([0., 4.], [20., 20.1]))),
+            activity=models.Activity.types['Seated'],
+        ),
+    )
+
 
 @pytest.mark.parametrize(
     "time, expected_co2_concentration", [
@@ -114,4 +127,30 @@ def root_mean_square_error_percentage(actual, predicted) -> float:
 
     acceptable_rmsep = 10 # Threshold of 10% for the accepted error margin
     assert rmsep <= acceptable_rmsep, f"RMSEP {rmsep} exceeds acceptable threshold {acceptable_rmsep}"
+
+@pytest.mark.parametrize("time", [4.1,10])
+def test_concentration_limit_last_state_change(simple_co2_conc_model, time):
+    npt.assert_almost_equal(simple_co2_conc_model._normed_concentration_limit(time), simple_co2_conc_model.min_background_concentration()/simple_co2_conc_model.normalization_factor())
+
+@pytest.mark.parametrize([
+    "start",
+    "stop"],
+    [
+        [10, 19],
+        [4, 19],
+        [0., 19],
+        [0., 20],
+    ]
+)
+def test_concentration_after_last_state_change(simple_co2_conc_model, simple_co2_conc_model_extended_presence, start, stop):
+    """
+    Test that the concentration results for a model with no more presence of emitters  
+    equals the concentration results of a model where the emitter will reenter at a later point.
+    """
+    time = (start+stop)/2
+    npt.assert_almost_equal(simple_co2_conc_model.removal_rate(time), simple_co2_conc_model_extended_presence.removal_rate(time))
+    npt.assert_almost_equal(simple_co2_conc_model.concentration(time), simple_co2_conc_model_extended_presence.concentration(time))
+    npt.assert_almost_equal(simple_co2_conc_model._normed_concentration(time), simple_co2_conc_model_extended_presence._normed_concentration(time))
+    npt.assert_almost_equal(simple_co2_conc_model.normed_integrated_concentration(start, stop), simple_co2_conc_model_extended_presence.normed_integrated_concentration(start, stop))
+    npt.assert_almost_equal(simple_co2_conc_model.integrated_concentration(start, stop), simple_co2_conc_model_extended_presence.integrated_concentration(start, stop))
 

caimira/tests/models/test_concentration_model.py

@@ -107,6 +107,26 @@ def simple_conc_model(data_registry):
         evaporation_factor=0.3,
     )
 
+@pytest.fixture
+def simple_conc_model_extended_presence(data_registry):
+    ventilation_times = models.SpecificInterval(([0.5, 1.], [1.1, 2], [2., 3.]), )
+    return models.ConcentrationModel(
+        data_registry=data_registry,
+        room = models.Room(75, models.PiecewiseConstant((0., 24.), (293,))),
+        ventilation = models.AirChange(ventilation_times, 100),
+        infected = models.InfectedPopulation(
+            data_registry=data_registry,
+            number=1,
+            presence=models.SpecificInterval(([0.5, 1.], [1.1, 2], [2., 3.], [20., 20.001]), ),
+            mask=models.Mask.types['Type I'],
+            activity=models.Activity.types['Seated'],
+            virus=models.Virus.types['SARS_CoV_2'],
+            expiration=models.Expiration.types['Breathing'],
+            host_immunity=0.,
+        ),
+        evaporation_factor=0.3,
+    )
+
 
 @pytest.fixture
 def dummy_population(simple_conc_model) -> models.Population:
@@ -143,35 +163,6 @@ def test_last_state_change_time(
     assert simple_conc_model.last_state_change(float(time)) == expected_last_state_change
 
 
-@pytest.mark.parametrize(
-    "time, expected_next_state_change", [
-        [0.0, 0.0],
-        [0.5, 0.5],
-        [1, 1],
-        [1.05, 1.1],
-        [1.1, 1.1],
-        [1.11, 2],
-        [2, 2],
-        [2.1, 3],
-        [3, 3],
-    ]
-)
-def test_next_state_change_time(
-        simple_conc_model: models.ConcentrationModel,
-        time,
-        expected_next_state_change,
-):
-    assert simple_conc_model._next_state_change(float(time)) == expected_next_state_change
-
-
-def test_next_state_change_time_out_of_range(simple_conc_model: models.ConcentrationModel):
-    with pytest.raises(
-            ValueError,
-            match=re.escape("The requested time (3.1) is greater than last available state change time (3.0)")
-    ):
-        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
 
@@ -286,3 +277,29 @@ def test_zero_ventilation_rate(
 
     normed_concentration = known_conc_model.concentration(1)
     assert normed_concentration == pytest.approx(expected_concentration, abs=1e-6)
+
+@pytest.mark.parametrize("time", [3.1,10])
+def test_concentration_limit_last_state_change(simple_conc_model, time):
+    npt.assert_almost_equal(simple_conc_model._normed_concentration_limit(time), simple_conc_model.min_background_concentration()/simple_conc_model.normalization_factor())
+
+@pytest.mark.parametrize([
+    "start",
+    "stop"],
+    [
+        [10, 19],
+        [3, 19],
+        [0., 19],
+        [0., 20],
+    ]
+)
+def test_concentration_after_last_state_change(simple_conc_model, simple_conc_model_extended_presence, start, stop):
+    """
+    Test that the concentration results for a model with no more presence of emitters  
+    equals the concentration results of a model where the emitter will reenter at a later point.
+    """
+    time = (start+stop)/2
+    npt.assert_almost_equal(simple_conc_model.removal_rate(time), simple_conc_model_extended_presence.removal_rate(time))
+    npt.assert_almost_equal(simple_conc_model.concentration(time), simple_conc_model_extended_presence.concentration(time))
+    npt.assert_almost_equal(simple_conc_model._normed_concentration(time), simple_conc_model_extended_presence._normed_concentration(time))
+    npt.assert_almost_equal(simple_conc_model.normed_integrated_concentration(start, stop), simple_conc_model_extended_presence.normed_integrated_concentration(start, stop))
+    npt.assert_almost_equal(simple_conc_model.integrated_concentration(start, stop), simple_conc_model_extended_presence.integrated_concentration(start, stop))

caimira/tests/models/test_exposure_model.py

@@ -19,18 +19,19 @@ class KnownNormedconcentration(models.ConcentrationModel):
 
     """
     normed_concentration_function: typing.Callable = lambda x: 0
+    normed_concentration_limit_function: typing.Callable = normed_concentration_function
 
     def removal_rate(self, time: float) -> models._VectorisedFloat:
         # Very large decay constant -> same as constant concentration
         return 1.e50
 
     def _normed_concentration_limit(self, time: float) -> models._VectorisedFloat:
-        return self.normed_concentration_function(time) * self.infected.number
+        return self.normed_concentration_limit_function(time) * self.infected.number
 
     def state_change_times(self):
         return [0., 24.]
 
-    def _next_state_change(self, time: float):
+    def _next_state_change(self, time: float): # Outdated method, see merge_requests/539
         return 24.
 
     def _normed_concentration(self, time: float) -> models._VectorisedFloat:  # noqa
@@ -56,7 +57,7 @@ def _normed_concentration(self, time: float) -> models._VectorisedFloat:  # noqa
     ),
 ]
 
-def known_concentrations(func, data_registry=DataRegistry()):
+def known_concentrations(func, func_lim=None, data_registry=DataRegistry()):
     dummy_room = models.Room(50, 0.5)
     dummy_ventilation = models._VentilationBase()
     dummy_infected_population = models.InfectedPopulation(
@@ -71,8 +72,13 @@ def known_concentrations(func, data_registry=DataRegistry()):
     )
     normed_func = lambda x: (func(x) /
         dummy_infected_population.emission_rate_per_person_when_present())
+    if func_lim:
+        normed_func_lim = lambda x: (func_lim(x) /
+            dummy_infected_population.emission_rate_per_person_when_present())
+    else:
+        normed_func_lim = normed_func
     return KnownNormedconcentration(data_registry, dummy_room, dummy_ventilation,
-                                dummy_infected_population, 0.3, normed_func)
+                                dummy_infected_population, 0.3, normed_func, normed_func_lim)
 
 
 @pytest.mark.parametrize(
@@ -114,8 +120,13 @@ def test_exposure_model_ndarray(data_registry, population, cm,
         [populations[2], np.array([1.36390289, 1.52436206])],
     ])
 def test_exposure_model_ndarray_and_float_mix(data_registry, population, expected_deposited_exposure, sr_model, cases_model):
-    cm = known_concentrations(
-        lambda t: 0. if np.floor(t) % 2 else np.array([0.6, 0.6]))
+    func = lambda t: 0. if np.floor(t) % 2 else np.array([0.6, 0.6])
+
+    # After merge_requests/539 normed_integrated_concentration computes normed_concentration_limit(t).
+    # However, the expected_deposited_exposure values assume normed_integrated_concentration computes
+    # normed_concentration_limit(_next_state_change(t)) = np.array([0.6, 0.6]) for all t
+    func_lim = lambda t: func(24)
+    cm = known_concentrations(func,func_lim)
     model = ExposureModel(data_registry, cm, sr_model, population, cases_model)
 
     np.testing.assert_almost_equal(