Vectorise PiecewiseConstant and WindowOpenings.

This commit is contained in:
Phil Elson 2021-04-07 10:15:48 +02:00
parent 7152be47e6
commit eccc57702d
8 changed files with 136 additions and 96 deletions

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,7 +135,7 @@ 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:
@ -141,7 +143,7 @@ class PiecewiseConstant:
if tuple(sorted(set(self.transition_times))) != self.transition_times:
raise ValueError("transition_times should not contain duplicated elements and should be sorted")
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 +161,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 +256,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
@ -322,10 +324,10 @@ 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
@ -339,19 +341,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
window_angle = 2.*np.rad2deg(np.arcsin(self.opening_length/(2.*self.window_height)))
return cd_max*(1-np.exp(-M*window_angle))

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

View file

@ -0,0 +1,78 @@
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))
@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()

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

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,34 @@ 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])},
{'window_width': np.array([0.15, 0.20])},
{'opening_length': np.array([0.15, 0.20])},
{'outside_temp': models.PiecewiseConstant((1, 2, 3), (np.array([20, 30]), 25))},
{'outside_temp': np.array([20, 30])},
{'inside_temp': np.array([20, 30])},
]
)
def test_HingedWindow_vectorisation(override_params):
defaults = {
'window_height': 0.15,
'window_width': 0.15,
'opening_length': 0.15,
'inside_temp': 20,
'outside_temp': 10,
}
defaults.update(override_params)
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)
def test_sliding_window(baseline_slidingwindow):
assert baseline_slidingwindow.discharge_coefficient == 0.6

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

View file

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

View file

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