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

Merge branch 'develop/refine_time_mesh_for_temp' into 'master'

Browse source 
Refine the time mesh for Geneva temperatures

Closes #72

See merge request cara/cara!70
This commit is contained in:
Philip James Elson 2020-11-10 07:33:31 +00:00
commit 050c9fd9ba
5 changed files with 131 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
cara/apps/calculator/model_generator.py
View file

@ -3,6 +3,7 @@ import html
import typing
from cara import models
from cara import data
@dataclass
@ -131,7 +132,7 @@ class FormData:
month = self.recurrent_event_month[:3]
inside_temp = models.PiecewiseConstant((0, 24), (293,))
outside_temp = models.GenevaTemperatures[month]
outside_temp = data.GenevaTemperatures[month]
ventilation = models.WindowOpening(
active=window_interval,

44
cara/data.py Normal file
View file

@ -0,0 +1,44 @@
import numpy as np
from cara import models
# average temperature of each month, hour per hour (from midnight to 11 pm)
Geneva_hourly_temperatures_celsius_per_hour = {
'Jan': [0.2, -0.3, -0.5, -0.9, -1.1, -1.4, -1.5, -1.5, -1.1, 0.1, 1.5,
2.8, 3.8, 4.4, 4.5, 4.4, 4.4, 3.9, 3.1, 2.7, 2.2, 1.7, 1.5, 1.1],
'Feb': [0.9, 0.3, 0.0, -0.5, -0.7, -1.1, -1.2, -1.1, -0.7, 0.8, 2.5,
4.2, 5.4, 6.2, 6.3, 6.2, 6.1, 5.5, 4.5, 4.1, 3.5, 2.8, 2.5, 2.0],
'Mar': [4.2, 3.5, 3.1, 2.5, 2.1, 1.6, 1.5, 1.6, 2.2, 4.0, 6.3, 8.4,
10.0, 11.1, 11.2, 11.1, 11.0, 10.2, 8.9, 8.3, 7.5, 6.7, 6.3, 5.6],
'Apr': [7.4, 6.7, 6.2, 5.5, 5.2, 4.7, 4.5, 4.6, 5.3, 7.2, 9.6, 11.9,
13.7, 14.8, 14.9, 14.8, 14.7, 13.8, 12.4, 11.8, 10.9, 10.1, 9.6, 8.9],
'May': [11.8, 11.1, 10.6, 9.9, 9.5, 8.9, 8.8, 8.9, 9.6, 11.6, 14.2, 16.6,
18.4, 19.6, 19.7, 19.6, 19.4, 18.6, 17.1, 16.5, 15.6, 14.6, 14.2, 13.4],
'Jun': [15.2, 14.4, 13.9, 13.2, 12.7, 12.2, 12.0, 12.1, 12.8, 15.0, 17.7,
20.2, 22.1, 23.3, 23.5, 23.4, 23.2, 22.3, 20.8, 20.1, 19.1, 18.2, 17.7, 16.9],
'Jul': [17.6, 16.7, 16.1, 15.3, 14.9, 14.3, 14.1, 14.2, 15.0, 17.3, 20.2,
23.0, 25.0, 26.3, 26.5, 26.4, 26.2, 25.2, 23.6, 22.8, 21.8, 20.8, 20.2, 19.4],
'Aug': [17.1, 16.2, 15.7, 14.9, 14.5, 13.9, 13.7, 13.8, 14.6, 16.9, 19.7,
22.4, 24.4, 25.6, 25.8, 25.7, 25.5, 24.5, 22.9, 22.2, 21.2, 20.2, 19.7, 18.9],
'Sep': [13.4, 12.7, 12.2, 11.5, 11.2, 10.7, 10.5, 10.6, 11.3, 13.2, 15.6,
17.9, 19.6, 20.8, 20.9, 20.8, 20.7, 19.8, 18.4, 17.8, 16.9, 16.1, 15.6, 14.9],
'Oct': [9.4, 8.8, 8.5, 7.9, 7.6, 7.2, 7.1, 7.2, 7.7, 9.3, 11.2, 13.0,
14.4, 15.3, 15.4, 15.3, 15.2, 14.5, 13.4, 12.9, 12.2, 11.6, 11.2, 10.6],
'Nov': [4.0, 3.6, 3.3, 2.9, 2.6, 2.3, 2.2, 2.2, 2.7, 3.9, 5.5, 6.9, 8.0,
8.7, 8.8, 8.7, 8.7, 8.1, 7.2, 6.8, 6.3, 5.7, 5.5, 5.0],
'Dec': [1.4, 1.0, 0.8, 0.4, 0.2, -0.0, -0.1, -0.1, 0.3, 1.3, 2.6, 3.8,
4.7, 5.2, 5.3, 5.2, 5.2, 4.7, 4.0, 3.7, 3.2, 2.8, 2.6, 2.2]
}
# Geneva hourly temperatures as piecewise constant function (in Kelvin)
GenevaTemperatures_hourly = {
month: models.PiecewiseConstant(tuple(np.arange(25.)),
tuple(273.15+np.array(temperatures)))
for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items()
}
# same temperatures on a finer temperature mesh
GenevaTemperatures = {
month: GenevaTemperatures_hourly[month].refine(refine_factor=10)
for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items()
}

46
cara/models.py
View file

@ -6,34 +6,6 @@ from abc import abstractmethod
from dataclasses import dataclass
# average temperature of each month, hour per hour (from midnight to 11 pm)
Geneva_hourly_temperatures_celsius_per_hour = {
'Jan': [0.2, -0.3, -0.5, -0.9, -1.1, -1.4, -1.5, -1.5, -1.1, 0.1, 1.5,
2.8, 3.8, 4.4, 4.5, 4.4, 4.4, 3.9, 3.1, 2.7, 2.2, 1.7, 1.5, 1.1],
'Feb': [0.9, 0.3, 0.0, -0.5, -0.7, -1.1, -1.2, -1.1, -0.7, 0.8, 2.5,
4.2, 5.4, 6.2, 6.3, 6.2, 6.1, 5.5, 4.5, 4.1, 3.5, 2.8, 2.5, 2.0],
'Mar': [4.2, 3.5, 3.1, 2.5, 2.1, 1.6, 1.5, 1.6, 2.2, 4.0, 6.3, 8.4,
10.0, 11.1, 11.2, 11.1, 11.0, 10.2, 8.9, 8.3, 7.5, 6.7, 6.3, 5.6],
'Apr': [7.4, 6.7, 6.2, 5.5, 5.2, 4.7, 4.5, 4.6, 5.3, 7.2, 9.6, 11.9,
13.7, 14.8, 14.9, 14.8, 14.7, 13.8, 12.4, 11.8, 10.9, 10.1, 9.6, 8.9],
'May': [11.8, 11.1, 10.6, 9.9, 9.5, 8.9, 8.8, 8.9, 9.6, 11.6, 14.2, 16.6,
18.4, 19.6, 19.7, 19.6, 19.4, 18.6, 17.1, 16.5, 15.6, 14.6, 14.2, 13.4],
'Jun': [15.2, 14.4, 13.9, 13.2, 12.7, 12.2, 12.0, 12.1, 12.8, 15.0, 17.7,
20.2, 22.1, 23.3, 23.5, 23.4, 23.2, 22.3, 20.8, 20.1, 19.1, 18.2, 17.7, 16.9],
'Jul': [17.6, 16.7, 16.1, 15.3, 14.9, 14.3, 14.1, 14.2, 15.0, 17.3, 20.2,
23.0, 25.0, 26.3, 26.5, 26.4, 26.2, 25.2, 23.6, 22.8, 21.8, 20.8, 20.2, 19.4],
'Aug': [17.1, 16.2, 15.7, 14.9, 14.5, 13.9, 13.7, 13.8, 14.6, 16.9, 19.7,
22.4, 24.4, 25.6, 25.8, 25.7, 25.5, 24.5, 22.9, 22.2, 21.2, 20.2, 19.7, 18.9],
'Sep': [13.4, 12.7, 12.2, 11.5, 11.2, 10.7, 10.5, 10.6, 11.3, 13.2, 15.6,
17.9, 19.6, 20.8, 20.9, 20.8, 20.7, 19.8, 18.4, 17.8, 16.9, 16.1, 15.6, 14.9],
'Oct': [9.4, 8.8, 8.5, 7.9, 7.6, 7.2, 7.1, 7.2, 7.7, 9.3, 11.2, 13.0,
14.4, 15.3, 15.4, 15.3, 15.2, 14.5, 13.4, 12.9, 12.2, 11.6, 11.2, 10.6],
'Nov': [4.0, 3.6, 3.3, 2.9, 2.6, 2.3, 2.2, 2.2, 2.7, 3.9, 5.5, 6.9, 8.0,
8.7, 8.8, 8.7, 8.7, 8.1, 7.2, 6.8, 6.3, 5.7, 5.5, 5.0],
'Dec': [1.4, 1.0, 0.8, 0.4, 0.2, -0.0, -0.1, -0.1, 0.3, 1.3, 2.6, 3.8,
4.7, 5.2, 5.3, 5.2, 5.2, 4.7, 4.0, 3.7, 3.2, 2.8, 2.6, 2.2]
}
@dataclass(frozen=True)
class Room:
@ -103,6 +75,9 @@ class PeriodicInterval(Interval):
@dataclass(frozen=True)
class PiecewiseConstant:
# TODO: implement rather a periodic version (24-hour period), where
# transition_times and values have the same length.
#: transition times at which the function changes value (hours).
transition_times: typing.Tuple[float, ...]
@ -135,13 +110,14 @@ class PiecewiseConstant:
present_times.append((t1,t2))
return SpecificInterval(present_times=present_times)
# Geneva hourly temperatures as piecewise constant function (in Kelvin)
GenevaTemperatures = {
month: PiecewiseConstant(tuple(np.arange(25.)),
tuple(273.15+np.array(temperatures)))
for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items()
}
def refine(self,refine_factor=10):
# 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)
return PiecewiseConstant(tuple(refined_times),
tuple(np.interp(refined_times[:-1],self.transition_times,
self.values+(self.values[-1],) ) ) )
@dataclass(frozen=True)

7
cara/tests/apps/calculator/test_model_generator.py
View file

@ -2,6 +2,7 @@ import pytest
from cara.apps.calculator import model_generator
from cara import models
from cara import data
import numpy as np
@pytest.fixture
@ -25,7 +26,7 @@ def test_ventilation_window(baseline_form):
window = models.WindowOpening(
active=models.PeriodicInterval(period=120, duration=10),
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
outside_temp=models.GenevaTemperatures['Dec'],
outside_temp=data.GenevaTemperatures['Dec'],
cd_b=0.6, window_height=1.6, opening_length=0.6,
)
baseline_form.ventilation_type = 'natural'
@ -75,7 +76,7 @@ def test_ventilation_window_hepa(baseline_form):
window = models.WindowOpening(
active=models.PeriodicInterval(period=120, duration=10),
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
outside_temp=models.GenevaTemperatures['Dec'],
outside_temp=data.GenevaTemperatures['Dec'],
cd_b=0.6, window_height=1.6, opening_length=0.6,
)
hepa = models.HEPAFilter(
@ -113,4 +114,4 @@ def test_present_intervals(baseline_form):
def test_key_validation(baseline_form_data):
baseline_form_data['activity_type'] = 'invalid key'
with pytest.raises(ValueError):
model_generator.FormData.from_dict(baseline_form_data)
model_generator.FormData.from_dict(baseline_form_data)

77
cara/tests/test_known_quantities.py
View file

@ -3,6 +3,7 @@ import numpy.testing as npt
import pytest
import cara.models as models
import cara.data as data
def test_no_mask_aerosols(baseline_model):
@ -245,6 +246,15 @@ def test_piecewiseconstant(time, expected_value):
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,)
@ -267,7 +277,7 @@ def test_piecewiseconstant_vs_interval(time):
def test_piecewiseconstant_transition_times():
outside_temp=models.GenevaTemperatures['Jan']
outside_temp=data.GenevaTemperatures['Jan']
assert set(outside_temp.transition_times) == outside_temp.interval().transition_times()
@ -289,13 +299,27 @@ def test_windowopening(time, expected_value):
def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),),
intervals_presence_infected=((0, 4), (5, 7.5))):
intervals_presence_infected=((0, 4), (5, 7.5)),
artificial_refinement=False,
temperatures=data.GenevaTemperatures_hourly):
if artificial_refinement:
# 5-fold increase of number of times, WITHOUT interpolation
# (hence transparent for the results)
refine_factor = 2
times_refined = tuple(np.linspace(0.,24,
refine_factor*len(temperatures[month].values)+1))
temperatures_refined = tuple(np.hstack([[v]*refine_factor
for v in temperatures[month].values]))
outside_temp = models.PiecewiseConstant(times_refined,temperatures_refined)
else:
outside_temp = temperatures[month]
model = models.Model(
room=models.Room(volume=75),
ventilation=models.WindowOpening(
active=models.SpecificInterval(intervals_open),
inside_temp=models.PiecewiseConstant((0,24),(293,)),
outside_temp=models.GenevaTemperatures[month],
outside_temp=outside_temp,
cd_b=0.6, window_height=1.6, opening_length=0.6,
),
infected=models.InfectedPerson(
@ -340,7 +364,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5
models.WindowOpening(
active=models.SpecificInterval(intervals_open),
inside_temp=models.PiecewiseConstant((0,24),(293,)),
outside_temp=models.GenevaTemperatures[month],
outside_temp=data.GenevaTemperatures[month],
cd_b=0.6, window_height=1.6, opening_length=0.6,
),
models.HEPAFilter(
@ -366,7 +390,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5
@pytest.mark.parametrize(
"month, temperatures",
models.Geneva_hourly_temperatures_celsius_per_hour.items(),
data.Geneva_hourly_temperatures_celsius_per_hour.items(),
)
@pytest.mark.parametrize(
"time",
@ -381,7 +405,7 @@ def test_concentrations_hourly_dep_temp_vs_constant(month, temperatures, time):
@pytest.mark.parametrize(
"month, temperatures",
models.Geneva_hourly_temperatures_celsius_per_hour.items(),
data.Geneva_hourly_temperatures_celsius_per_hour.items(),
)
@pytest.mark.parametrize(
"time",
@ -402,7 +426,7 @@ def test_concentrations_hourly_dep_multipleventilation():
@pytest.mark.parametrize(
"month_temp_item",
models.Geneva_hourly_temperatures_celsius_per_hour.items(),
data.Geneva_hourly_temperatures_celsius_per_hour.items(),
)
@pytest.mark.parametrize(
"time",
@ -414,3 +438,42 @@ def test_concentrations_hourly_dep_adding_artificial_transitions(month_temp_item
m1 = build_hourly_dependent_model(month,intervals_open=((7.5, 8.5),))
m2 = build_hourly_dependent_model(month,intervals_open=((7.5, 8.5),(8.,8.1)))
npt.assert_allclose(m1.concentration(time), m2.concentration(time), rtol=1e-5)
@pytest.mark.parametrize(
"time",
list(np.random.random_sample(10)*24.)+list(np.arange(0,24.5,0.5)),
)
def test_concentrations_refine_times(time):
month = 'Jan'
m1 = build_hourly_dependent_model(month,intervals_open=((0, 24),))
m2 = build_hourly_dependent_model(month,intervals_open=((0, 24),),
artificial_refinement=True)
npt.assert_allclose(m1.concentration(time), m2.concentration(time), rtol=1e-8)
@pytest.mark.parametrize(
"month, expected_r0",
[
['Jan', 91.06953],
['Jun', 99.46692],
],
)
def test_r0_hourly_dep(month,expected_r0):
m = build_hourly_dependent_model(month,intervals_open=((0,24),),
intervals_presence_infected=((8,12),(13,17)))
p = m.infection_probability()
npt.assert_allclose(p, expected_r0)
@pytest.mark.parametrize(
"month, expected_r0",
[
['Jan', 91.19912],
['Jun', 99.59226],
],
)
def test_r0_hourly_dep_refined(month,expected_r0):
m = build_hourly_dependent_model(month,intervals_open=((0,24),),
intervals_presence_infected=((8,12),(13,17)),
temperatures=data.GenevaTemperatures)
p = m.infection_probability()
npt.assert_allclose(p, expected_r0)