diff --git a/caimira/apps/expert.py b/caimira/apps/expert.py index c15234c0..5162d8f0 100644 --- a/caimira/apps/expert.py +++ b/caimira/apps/expert.py @@ -126,7 +126,7 @@ class ExposureModelResult(View): self.figure.canvas, ]) - def initialize_axes(self) -> matplotlib.figure.Axes: + def initialize_axes(self) -> typing.Tuple[matplotlib.axes.Axes, matplotlib.axes.Axes]: ax = self.figure.add_subplot(1, 1, 1) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) @@ -230,7 +230,7 @@ class ExposureComparisonResult(View): # unless the widget is wrapped in a container (it is seen on all tabs otherwise!). return widgets.HBox([self.figure.canvas]) - def initialize_axes(self) -> matplotlib.figure.Axes: + def initialize_axes(self) -> typing.Tuple[matplotlib.axes.Axes, matplotlib.axes.Axes]: ax = self.figure.add_subplot(1, 1, 1) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) diff --git a/caimira/apps/expert_co2.py b/caimira/apps/expert_co2.py index 345f7cf7..5c3437bf 100644 --- a/caimira/apps/expert_co2.py +++ b/caimira/apps/expert_co2.py @@ -69,7 +69,7 @@ class ExposureModelResult(View): self.figure.canvas, ]) - def initialize_axes(self) -> matplotlib.figure.Axes: + def initialize_axes(self) -> matplotlib.axes.Axes: ax = self.figure.add_subplot(1, 1, 1) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) @@ -140,7 +140,7 @@ class ExposureComparisonResult(View): # unless the widget is wrapped in a container (it is seen on all tabs otherwise!). return widgets.HBox([self.figure.canvas]) - def initialize_axes(self) -> matplotlib.figure.Axes: + def initialize_axes(self) -> matplotlib.axes.Axes: ax = self.figure.add_subplot(1, 1, 1) ax.spines[['right', 'top']].set_visible(False) diff --git a/caimira/models.py b/caimira/models.py index e8de7667..24ef62e2 100644 --- a/caimira/models.py +++ b/caimira/models.py @@ -1055,8 +1055,15 @@ class _ConcentrationModelBase: """ V = self.room.volume RR = self.removal_rate(time) - - return (self.population.people_present(time) / (RR * V) + + + if isinstance(RR, np.ndarray): + invRR = np.empty(RR.shape, dtype=np.float64) + invRR[RR == 0.] = np.nan + invRR[RR != 0.] = 1. / RR[RR != 0.] + else: + invRR = np.nan if RR == 0. else 1. / RR # type: ignore + + return (self.population.people_present(time) * invRR / V + self.min_background_concentration()/self.normalization_factor()) @method_cache @@ -1132,23 +1139,26 @@ class _ConcentrationModelBase: 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) - # If RR is 0, conc_limit does not play a role but its computation - # would raise an error -> we set it to zero. - try: - conc_limit = self._normed_concentration_limit(next_state_change_time) - except ZeroDivisionError: - conc_limit = 0. t_last_state_change = self.last_state_change(time) conc_at_last_state_change = self._normed_concentration_cached(t_last_state_change) - delta_time = time - t_last_state_change + fac = np.exp(-RR * delta_time) + if isinstance(RR, np.ndarray): + 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.]) + 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) - return conc_limit * (1 - fac) + conc_at_last_state_change * fac - + return curr_conc_state + conc_at_last_state_change * fac + def concentration(self, time: float) -> _VectorisedFloat: """ Total concentration as a function of time. The normalization @@ -1260,9 +1270,7 @@ class CO2ConcentrationModel(_ConcentrationModelBase): return self.CO2_emitters def removal_rate(self, time: float) -> _VectorisedFloat: - # Setting minimum air exchange rate to 1e-6 to avoid divisions by - # zero when computing the CO2 concentration. - return np.maximum(1e-6,self.ventilation.air_exchange(self.room, time)) + return self.ventilation.air_exchange(self.room, time) def min_background_concentration(self) -> _VectorisedFloat: """ diff --git a/caimira/tests/models/test_concentration_model.py b/caimira/tests/models/test_concentration_model.py index f8f42106..27703a2f 100644 --- a/caimira/tests/models/test_concentration_model.py +++ b/caimira/tests/models/test_concentration_model.py @@ -252,8 +252,11 @@ def test_normed_integrated_concentration_vectorisation( "known_min_background_concentration", "expected_concentration"], [ - [0., 240., 240.], - [0., np.array([240., 240.]), np.array([240., 240.])] + [0., 240., 240. + 0.5/75], + [0.0001, 240.0, 240. + 0.5/75], + [1e-6, 240.0, 240 + 0.5/75], + [0., np.array([240., 240.]), np.array([240. + 0.5/75, 240. + 0.5/75])], + [np.array([0.0001, 1e-6]), np.array([240., 240.]), np.array([240. + 0.5/75, 240. + 0.5/75])], ] ) def test_zero_ventilation_rate( @@ -272,4 +275,5 @@ def test_zero_ventilation_rate( known_min_background_concentration = known_min_background_concentration) normed_concentration = known_conc_model.concentration(1) - npt.assert_almost_equal(normed_concentration, expected_concentration) + assert normed_concentration == pytest.approx(expected_concentration, abs=1e-6) + \ No newline at end of file