- CARA is made available for internal CERN use only. - It is intended for Members of Personnel with roles related to Supervision, Health & Safety or Space Management, in order to simulate the concerned workplaces on CERN sites. - For use outside of this scope, it has to be performed under a pre-defined framework and licence agreement issued by CERN Knowledge Transfer (kt@cern.ch). + At CERN, CARA is intended for Members of Personnel with roles related to Supervision, Health & Safety or Space Management, in order to simulate the concerned workplaces on CERN sites.
{% endblock disclaimer %} diff --git a/cara/apps/templates/layout.html.j2 b/cara/apps/templates/layout.html.j2 index e7658ad2..566c545d 100644 --- a/cara/apps/templates/layout.html.j2 +++ b/cara/apps/templates/layout.html.j2 @@ -193,7 +193,6 @@ policy issues. Any initiative is conducted on a best effort and as-is basis, without liability or warranty. - @@ -231,9 +230,22 @@ ++ CARA is Apache 2.0 licensed open-source + software developed at CERN. + You can find the source code at https://gitlab.cern.ch/cara/cara, + where we welcome contributions, feature requests and issue reports. +
+ @@ -283,7 +295,7 @@ © - 2020 + 2020 - CERN diff --git a/cara/models.py b/cara/models.py index bee4a0ff..d761b0c3 100644 --- a/cara/models.py +++ b/cara/models.py @@ -31,9 +31,11 @@ the same for all parameters of a single model. """ from dataclasses import dataclass -import numpy as np import typing +import numpy as np +from scipy.interpolate import interp1d + if not typing.TYPE_CHECKING: from memoization import cached else: @@ -133,15 +135,18 @@ class PiecewiseConstant: transition_times: typing.Tuple[float, ...] #: values of the function between transitions - values: typing.Tuple[float, ...] + values: typing.Tuple[_VectorisedFloat, ...] def __post_init__(self): if len(self.transition_times) != len(self.values)+1: raise ValueError("transition_times should contain one more element than values") if tuple(sorted(set(self.transition_times))) != self.transition_times: raise ValueError("transition_times should not contain duplicated elements and should be sorted") + shapes = [np.array(v).shape for v in self.values] + if not all(shapes[0] == shape for shape in shapes): + raise ValueError("All values must have the same shape") - def value(self, time) -> float: + def value(self, time) -> _VectorisedFloat: if time <= self.transition_times[0]: return self.values[0] elif time > self.transition_times[-1]: @@ -159,21 +164,21 @@ class PiecewiseConstant: for t1, t2, value in zip(self.transition_times[:-1], self.transition_times[1:], self.values): if value: - present_times.append((t1,t2)) + present_times.append((t1, t2)) return SpecificInterval(present_times=tuple(present_times)) - def refine(self, refine_factor=10): + def refine(self, refine_factor=10) -> "PiecewiseConstant": # build a new PiecewiseConstant object with a refined mesh, # using a linear interpolation in-between the initial mesh points refined_times = np.linspace(self.transition_times[0], self.transition_times[-1], (len(self.transition_times)-1) * refine_factor+1) + interpolator = interp1d( + self.transition_times, + np.concatenate([self.values, self.values[-1:]], axis=0), + axis=0) return PiecewiseConstant( tuple(refined_times), - tuple(np.interp( - refined_times[:-1], - self.transition_times, - self.values + (self.values[-1], ), - )), + tuple(interpolator(refined_times)[:-1]), ) @@ -254,10 +259,10 @@ class WindowOpening(Ventilation): outside_temp: PiecewiseConstant #: The height of the window (m). - window_height: float + window_height: _VectorisedFloat #: The length of the opening-gap when the window is open (m). - opening_length: float + opening_length: _VectorisedFloat #: The number of windows of the given dimensions. number_of_windows: int = 1 @@ -266,7 +271,7 @@ class WindowOpening(Ventilation): min_deltaT: float = 0.1 @property - def discharge_coefficient(self) -> float: + def discharge_coefficient(self) -> _VectorisedFloat: """ Discharge coefficient (or cd_b): what portion effective area is used to exchange air (0 <= discharge_coefficient <= 1). @@ -286,14 +291,14 @@ class WindowOpening(Ventilation): return 0. # Reminder, no dependence on time in the resulting calculation. - inside_temp = self.inside_temp.value(time) - outside_temp = self.outside_temp.value(time) + inside_temp: _VectorisedFloat = self.inside_temp.value(time) + outside_temp: _VectorisedFloat = self.outside_temp.value(time) # The inside_temperature is forced to be always at least min_deltaT degree # warmer than the outside_temperature. Further research needed to # handle the buoyancy driven ventilation when the temperature gradient # is inverted. - inside_temp = max(inside_temp, outside_temp + self.min_deltaT) + inside_temp = np.maximum(inside_temp, outside_temp + self.min_deltaT) # type: ignore temp_gradient = (inside_temp - outside_temp) / outside_temp root = np.sqrt(9.81 * self.window_height * temp_gradient) window_area = self.window_height * self.opening_length * self.number_of_windows @@ -307,7 +312,7 @@ class SlidingWindow(WindowOpening): the horizontal plane). """ @property - def discharge_coefficient(self) -> float: + def discharge_coefficient(self) -> _VectorisedFloat: """ Average measured value of discharge coefficient for sliding or side-hung windows. @@ -322,14 +327,14 @@ class HingedWindow(WindowOpening): horizontal plane). """ #: Window width (m). - window_width: float = 0.0 + window_width: _VectorisedFloat = 0.0 def __post_init__(self): - if not self.window_width > 0: + if self.window_width is 0.0: raise ValueError('window_width must be set') @property - def discharge_coefficient(self) -> float: + def discharge_coefficient(self) -> _VectorisedFloat: """ Simple model to compute discharge coefficient for top or bottom hung hinged windows, in the absence of empirical test results @@ -339,19 +344,15 @@ class HingedWindow(WindowOpening): see Section 8.3 of BB101 and Section 11.3 of ESFA Output Specification Annex 2F on Ventilation opening areas. """ - window_ratio = self.window_width / self.window_height - if window_ratio < 0.5: - M = 0.06 - cd_max = 0.612 - elif window_ratio < 1: - M = 0.048 - cd_max = 0.589 - elif window_ratio < 2: - M = 0.04 - cd_max = 0.563 - else: - M = 0.038 - cd_max = 0.548 + window_ratio = np.array(self.window_width / self.window_height) + coefs = np.empty(window_ratio.shape + (2, ), dtype=np.float64) + + coefs[window_ratio < 0.5] = (0.06, 0.612) + coefs[np.bitwise_and(0.5 <= window_ratio, window_ratio < 1)] = (0.048, 0.589) + coefs[np.bitwise_and(1 <= window_ratio, window_ratio < 2)] = (0.04, 0.563) + coefs[window_ratio >= 2] = (0.038, 0.548) + M, cd_max = coefs.T + window_angle = 2.*np.rad2deg(np.arcsin(self.opening_length/(2.*self.window_height))) return cd_max*(1-np.exp(-M*window_angle)) diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index 3d4f0fdc..9be1cbbf 100644 --- a/cara/tests/models/test_concentration_model.py +++ b/cara/tests/models/test_concentration_model.py @@ -13,7 +13,6 @@ import cara.models {'coefficient_of_infectivity': np.array([0.02, 0.05])}, {'η_exhale': np.array([0.92, 0.95])}, {'η_leaks': np.array([0.15, 0.20])}, - ] ) def test_concentration_model_vectorisation(override_params): diff --git a/cara/tests/models/test_piecewiseconstant.py b/cara/tests/models/test_piecewiseconstant.py new file mode 100644 index 00000000..a0ba14f0 --- /dev/null +++ b/cara/tests/models/test_piecewiseconstant.py @@ -0,0 +1,85 @@ +import numpy as np +import pytest + +from cara import models +from cara import data + + +def test_piecewiseconstantfunction_wrongarguments(): + # number of values should be 1+number of transition times + pytest.raises(ValueError, models.PiecewiseConstant, (0, 1), (0, 0)) + pytest.raises(ValueError, models.PiecewiseConstant, (0,), (0, 0)) + # two transition times cannot be equal + pytest.raises(ValueError, models.PiecewiseConstant, (0, 2, 2), (0, 0)) + # unsorted transition times are not allowed + pytest.raises(ValueError, models.PiecewiseConstant, (2, 0), (0, 0)) + + # If vectors, must all be same length. + with pytest.raises(ValueError, match="All values must have the same shape"): + models.PiecewiseConstant( + (0, 8, 16), (np.array([5, 7]), np.array([8, 9, 10])), + ) + + + +@pytest.mark.parametrize( + "time, expected_value", + [ + [10, 5], + [20.5, 8], + [8, 2], + [0, 2], + [24, 8], + [-1, 2], + [25, 8], + ], +) +def test_piecewiseconstant(time, expected_value): + transition_times = (0, 8, 16, 24) + values = (2, 5, 8) + fun = models.PiecewiseConstant(transition_times, values) + assert fun.value(time) == expected_value + + +def test_piecewiseconstant_interp(): + transition_times = (0, 8, 16, 24) + values = (2, 5, 8) + refined_fun = models.PiecewiseConstant(transition_times, values).refine(refine_factor=2) + assert refined_fun.transition_times == (0, 4, 8, 12, 16, 20, 24) + assert refined_fun.values == (2, 3.5, 5, 6.5, 8, 8) + + +def test_piecewiseconstant_interp_vectorised(): + transition_times = (0, 8, 16, 24) + values = (np.array([2, 3]), np.array([5, 7]), np.array([8, 9])) + refined_fun = models.PiecewiseConstant(transition_times, values).refine(refine_factor=2) + assert refined_fun.transition_times == (0, 4, 8, 12, 16, 20, 24) + np.testing.assert_almost_equal( + refined_fun.values, ((2, 3), (3.5, 5), (5, 7), (6.5, 8), (8, 9), (8, 9)), + ) + + +def test_constantfunction(): + transition_times = (0, 24) + values = (20,) + fun = models.PiecewiseConstant(transition_times, values) + for t in [0, 1, 8, 10, 16, 20.1, 24]: + assert (fun.value(t) == 20) + + +@pytest.mark.parametrize( + "time", + [0, 1, 8, 10, 16, 20.1, 24], +) +def test_piecewiseconstant_vs_interval(time): + transition_times = (0, 8, 16, 24) + values = (0, 1, 0) + fun = models.PiecewiseConstant(transition_times, values) + interval = models.SpecificInterval(present_times=[(8,16)]) + assert interval.transition_times() == fun.interval().transition_times() + assert fun.interval().triggered(time) == interval.triggered(time) + + +def test_piecewiseconstant_transition_times(): + outside_temp = data.GenevaTemperatures['Jan'] + assert set(outside_temp.transition_times) == outside_temp.interval().transition_times() diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 042d9db0..cc3a6249 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -190,70 +190,6 @@ def test_expiration_aerosols(): npt.assert_allclose(exp1.aerosols(mask), exp2.aerosols(mask), rtol=1e-5) -def test_piecewiseconstantfunction_wrongarguments(): - # number of values should be 1+number of transition times - pytest.raises(ValueError,models.PiecewiseConstant,(0,1),(0,0)) - pytest.raises(ValueError,models.PiecewiseConstant,(0,),(0,0)) - # two transition times cannot be equal - pytest.raises(ValueError,models.PiecewiseConstant,(0,2,2),(0,0)) - # unsorted transition times are not allowed - pytest.raises(ValueError,models.PiecewiseConstant,(2,0),(0,0)) - - -@pytest.mark.parametrize( - "time, expected_value", - [ - [10, 5], - [20.5, 8], - [8, 2], - [0, 2], - [24, 8], - [-1, 2], - [25, 8], - ], -) -def test_piecewiseconstant(time, expected_value): - transition_times = (0, 8, 16, 24) - values = (2, 5, 8) - fun = models.PiecewiseConstant(transition_times,values) - assert fun.value(time) == expected_value - - -def test_piecewiseconstant_interp(): - transition_times = (0, 8, 16, 24) - values = (2, 5, 8) - fun = models.PiecewiseConstant(transition_times,values) - refined_fun = models.PiecewiseConstant(transition_times,values).refine(refine_factor=2) - assert refined_fun.transition_times == (0, 4, 8, 12, 16, 20 ,24) - assert refined_fun.values == (2, 3.5, 5, 6.5, 8, 8) - - -def test_constantfunction(): - transition_times = (0,24) - values = (20,) - fun = models.PiecewiseConstant(transition_times,values) - for t in [0,1,8,10,16,20.1,24]: - assert (fun.value(t) == 20) - - -@pytest.mark.parametrize( - "time", - [0,1,8,10,16,20.1,24], -) -def test_piecewiseconstant_vs_interval(time): - transition_times = (0,8,16,24) - values = (0,1,0) - fun = models.PiecewiseConstant(transition_times,values) - interval = models.SpecificInterval(present_times=[(8,16)]) - assert interval.transition_times() == fun.interval().transition_times() - assert fun.interval().triggered(time) == interval.triggered(time) - - -def test_piecewiseconstant_transition_times(): - outside_temp=data.GenevaTemperatures['Jan'] - assert set(outside_temp.transition_times) == outside_temp.interval().transition_times() - - @pytest.mark.parametrize( "time, expected_value", [ diff --git a/cara/tests/test_ventilation.py b/cara/tests/test_ventilation.py index d98fdff8..f10a6263 100644 --- a/cara/tests/test_ventilation.py +++ b/cara/tests/test_ventilation.py @@ -46,9 +46,8 @@ def test_number_of_windows(baseline_slidingwindow): [4., 0.308], ], ) -def test_hinged_window(baseline_hingedwindow,window_width, +def test_hinged_window(baseline_hingedwindow, window_width, expected_discharge_coefficient): - room = models.Room(75) hinged_window = dataclasses.replace(baseline_hingedwindow, window_width=window_width) @@ -56,9 +55,40 @@ def test_hinged_window(baseline_hingedwindow,window_width, expected_discharge_coefficient, rtol=1e-2) -def test_sliding_window(baseline_slidingwindow): - room = models.Room(75) +@pytest.mark.parametrize( + "override_params", [ + {'window_height': np.array([0.15, 0.20, 0.25])}, + {'window_width': np.array([0.15, 0.20, 0.25])}, + {'opening_length': np.array([0.15, 0.20, 0.25])}, + {'outside_temp': models.PiecewiseConstant( + (0, 2, 3), (np.array([20, 30, 28]), np.array([25, 30, 27])) + )}, + {'inside_temp': models.PiecewiseConstant( + (0, 20), (np.array([20, 30, 25]), ) + )}, + ] +) +def test_HingedWindow_vectorisation(override_params): + defaults = { + 'window_height': 0.15, + 'window_width': 0.15, + 'opening_length': 0.15, + 'inside_temp': models.PiecewiseConstant((0, 2, 3), (20, 25)), + 'outside_temp': models.PiecewiseConstant((0, 2, 3), (10, 15)), + } + defaults.update(override_params) + room = models.Room(volume=75) + t = 0.5 + window = models.HingedWindow(models.PeriodicInterval(60, 30), **defaults) + if {'window_height', 'opening_length', 'window_width'}.intersection(override_params): + assert isinstance(window.discharge_coefficient, np.ndarray) + else: + assert isinstance(window.discharge_coefficient, float) + assert isinstance(window.air_exchange(room, t), np.ndarray) + + +def test_sliding_window(baseline_slidingwindow): assert baseline_slidingwindow.discharge_coefficient == 0.6 diff --git a/requirements.txt b/requirements.txt index 2c87413a..1bbecfe1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -61,6 +61,7 @@ pyrsistent==0.17.3 python-dateutil==2.8.1 pyzmq==22.0.3 qrcode==6.1 +scipy==1.5.4 Send2Trash==1.5.0 six==1.15.0 sniffio==1.2.0 diff --git a/setup.cfg b/setup.cfg index 7dac8522..f40978b4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -19,3 +19,6 @@ ignore_missing_imports = True [mypy-qrcode.*] ignore_missing_imports = True + +[mypy-scipy.*] +ignore_missing_imports = True diff --git a/setup.py b/setup.py index bf44d292..8868599f 100644 --- a/setup.py +++ b/setup.py @@ -26,6 +26,7 @@ REQUIREMENTS: dict = { 'mistune', 'numpy', 'qrcode[pil]', + 'scipy', 'tornado', 'voila >=0.2.4', ], @@ -58,6 +59,7 @@ setup( classifiers=[ "Programming Language :: Python :: 3", "Operating System :: OS Independent", + "License :: OSI Approved :: Apache Software License", # Apache 2.0 ], install_requires=REQUIREMENTS['core'],