pingu_98/cara
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'feature/model_vectorisation_pt2' into 'master'

Browse source 
Vectorise PiecewiseConstant and WindowOpenings

See merge request cara/cara!159
This commit is contained in:
Nicolas Mounet 2021-04-29 04:35:10 +00:00
commit 2d00bc859f
8 changed files with 158 additions and 102 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

67
cara/models.py
View file

@ -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))

1
cara/tests/models/test_concentration_model.py
View file

@ -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):

85
cara/tests/models/test_piecewiseconstant.py Normal file
View file

@ -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()

64
cara/tests/test_known_quantities.py
View file

@ -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",
[

38
cara/tests/test_ventilation.py
View file

@ -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

1
requirements.txt
View file

@ -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

3
setup.cfg
View file

@ -19,3 +19,6 @@ ignore_missing_imports = True
[mypy-qrcode.*]
ignore_missing_imports = True
[mypy-scipy.*]
ignore_missing_imports = True

1
setup.py
View file

@ -26,6 +26,7 @@ REQUIREMENTS: dict = {
'mistune',
'numpy',
'qrcode[pil]',
'scipy',
'tornado',
'voila >=0.2.4',
],