Merge with master
This commit is contained in:
Luis Aleixo
commit
4704a6e022
12 changed files with 311 additions and 236 deletions
|
|
@ -33,7 +33,7 @@ from .user import AuthenticatedUser, AnonymousUser
|
|||
# calculator version. If the calculator needs to make breaking changes (e.g. change
|
||||
# form attributes) then it can also increase its MAJOR version without needing to
|
||||
# increase the overall CARA version (found at ``cara.__version__``).
|
||||
__version__ = "3.3.0"
|
||||
__version__ = "4.0.0"
|
||||
|
||||
|
||||
class BaseRequestHandler(RequestHandler):
|
||||
|
|
|
|||
|
|
@ -138,7 +138,7 @@ def calculate_report_data(model: models.ExposureModel):
|
|||
exposed_occupants = model.exposed.number
|
||||
expected_new_cases = np.array(model.expected_new_cases()).mean()
|
||||
cumulative_doses = np.cumsum([
|
||||
np.array(model.exposure_between_bounds(float(time1), float(time2))).mean()
|
||||
np.array(model.deposited_exposure_between_bounds(float(time1), float(time2))).mean()
|
||||
for time1, time2 in zip(times[:-1], times[1:])
|
||||
])
|
||||
|
||||
|
|
@ -337,11 +337,13 @@ class ReportGenerator:
|
|||
)
|
||||
context['permalink'] = generate_permalink(base_url, self.calculator_prefix, form)
|
||||
context['calculator_prefix'] = self.calculator_prefix
|
||||
|
||||
# For further information about these values visit https://gitlab.cern.ch/cara/cara/-/merge_requests/321.
|
||||
context['scale_warning'] = {
|
||||
'level': 'red-4',
|
||||
'incidence_rate': 'higher or equal to 100 new cases per 100 000 inhabitants',
|
||||
'onsite_access': 'lower than 4000',
|
||||
'threshold': '5%'
|
||||
'level': 'red-4', # 'red-4' - 'orange-3' - 'yellow-2' - 'green-1'
|
||||
'risk': 'strong', # 'strong' - 'medium' - 'reduced' - ''
|
||||
'onsite_access': '4’000', # '4’000' - '5’000' - '6’500' - '8’000'
|
||||
'threshold': '2%' # '2%' - '' - '' - ''
|
||||
}
|
||||
return context
|
||||
|
||||
|
|
|
|||
|
|
@ -7,9 +7,9 @@
|
|||
{% endblock report_preamble_navtab %}
|
||||
|
||||
{% block warning_animation %}
|
||||
{% if ((prob_inf > 15) or (expected_new_cases >= 1)) %} {% set warning_color= 'bg-danger' %}
|
||||
{% elif (5 <= prob_inf <= 15) %} {% set warning_color = 'bg-warning' %}
|
||||
{% elif (prob_inf < 5) %} {% set warning_color = 'bg-success' %}
|
||||
{% if ((prob_inf > 10) or (expected_new_cases >= 1)) %} {% set warning_color= 'bg-danger' %}
|
||||
{% elif (2 <= prob_inf <= 10) %} {% set warning_color = 'bg-warning' %}
|
||||
{% elif (prob_inf < 2) %} {% set warning_color = 'bg-success' %}
|
||||
{% endif %}
|
||||
|
||||
<div class="intro-banner-vdo-play-btn {{warning_color}} m-auto d-flex align-items-center justify-content-center">
|
||||
|
|
@ -23,41 +23,41 @@
|
|||
|
||||
{% block report_summary %}
|
||||
<div class="flex-row align-self-center">
|
||||
{% if ((prob_inf > 15) or (expected_new_cases >= 1)) %}
|
||||
{% if ((prob_inf > 10) or (expected_new_cases >= 1)) %}
|
||||
<div class="alert alert-danger mb-0" role="alert">
|
||||
<strong>Not Acceptable:</strong>
|
||||
Taking into account the uncertainties tied to the model variables, in this scenario, the <b>probability of one exposed occupant getting infected is {{ prob_inf | non_zero_percentage }}</b> and the <b>expected number of new cases is {{ expected_new_cases | float_format }}</b>*.
|
||||
|
||||
</div>
|
||||
{% elif 5 <= prob_inf <= 15 %}
|
||||
{% elif 2 <= prob_inf <= 10 %}
|
||||
<div class="alert alert-warning mb-0" role="alert">
|
||||
<strong>Attention:</strong>
|
||||
Taking into account the uncertainties tied to the model variables, in this scenario, the <b>probability of one exposed occupant getting infected is {{ prob_inf | non_zero_percentage }}</b> and the <b>expected number of new cases is {{ expected_new_cases | float_format }}</b>*.
|
||||
</div>
|
||||
{% elif prob_inf < 5 %}
|
||||
{% elif prob_inf < 2 %}
|
||||
<div class="alert alert-success mb-0" role="alert">
|
||||
<strong>Acceptable:</strong>
|
||||
Taking into account the uncertainties tied to the model variables, in this scenario, the <b>probability of one exposed occupant getting infected is {{ prob_inf | non_zero_percentage }}</b> and the <b>expected number of new cases is {{ expected_new_cases | float_format }}</b>*.
|
||||
</div>
|
||||
{% endif %}
|
||||
|
||||
{% if (prob_inf > 5) %}
|
||||
{% if (prob_inf > 2) %}
|
||||
<br>
|
||||
{% if scale_warning.level == "green-1" %}
|
||||
<div class="alert alert-dark mb-0" role="alert" style="height:fit-content">
|
||||
Note: the current CERN COVID Scale is <b>Green - 1</b>, which means the incidence rate in the local community is <b>{{scale_warning.incidence_rate}}</b>. Align your risk assessment with the guidance and instructions provided by the HSE Unit.
|
||||
Note: the current CERN COVID Scale is <b>Green – 1</b>. The risk of circulation of asymptomatic or pre-symptomatic infected individuals within the CERN site is considered negligible. There may be more than <b>{{scale_warning.onsite_access}} daily on-site accesses</b>. Align your risk assessment with the guidance and instructions provided by the HSE Unit.
|
||||
</div>
|
||||
{% elif scale_warning.level == "yellow-2" %}
|
||||
<div class="alert alert-dark mb-0" role="alert" style="height:fit-content">
|
||||
Note: the current CERN COVID Scale is <b>Yellow - 2</b>, which means the incidence rate in the local community is <b>{{scale_warning.incidence_rate}}</b>. There is a reduced chance that asymptomatic or pre-symptomatic infected individuals circulate within the CERN site which, during this stage, corresponds to an average daily on-site access {{scale_warning.onsite_access}}. See with your supervisor if this scenario is acceptable.
|
||||
Note: the current CERN COVID Scale is <b>Yellow - 2</b>. There is a {{scale_warning.risk}} risk that asymptomatic or pre-symptomatic infected individuals circulate within the CERN site. There may be around <b>{{scale_warning.onsite_access}} daily on-site accesses</b> during this stage. See with your supervisor, DSO/LEXGLIMOS and space manager if this scenario is acceptable and if any additional measures can be applied (ALARA).
|
||||
</div>
|
||||
{% elif scale_warning.level == "orange-3" %}
|
||||
<div class="alert alert-dark mb-0" role="alert" style="height:fit-content">
|
||||
Warning: the current CERN COVID Scale is <b>Orange - 3</b>, which means the incidence rate in the local community is <b>{{scale_warning.incidence_rate}}</b>. There is a slight chance that asymptomatic or pre-symptomatic infected individuals circulate within the CERN site which, during this stage, corresponds to an average daily on-site access {{scale_warning.onsite_access}}. See with your supervisor if any additional measures can be applied (ALARA).
|
||||
Warning: the current CERN COVID Scale is <b>Orange - 3</b>. There is a {{scale_warning.risk}} risk that asymptomatic or pre-symptomatic infected individuals circulate within the CERN site. There may be around <b>{{scale_warning.onsite_access}} daily on-site accesses</b> during this stage. See with your supervisor, DSO/LEXGLIMOS and space manager if this scenario is acceptable and if any additional measures can be applied (ALARA).
|
||||
</div>
|
||||
{% elif scale_warning.level == "red-4" %}
|
||||
<div class="alert alert-dark mb-0" role="alert" style="height:fit-content">
|
||||
Warning: the current CERN COVID Scale is <b>Red - 4</b>, which means the incidence rate in the local community is <b>{{scale_warning.incidence_rate}}</b>. There is a strong chance that asymptomatic or pre-symptomatic infected individuals circulate within the CERN site which, during this stage, corresponds to an average daily on-site access {{scale_warning.onsite_access}}. Please reduce the value below the threshold of <b>{{scale_warning.threshold}}</b>.
|
||||
Warning: the current CERN COVID Scale is <b>Red - 4</b>. There is a {{scale_warning.risk}} risk that asymptomatic or pre-symptomatic infected individuals circulate within the CERN site. There may be around <b>{{scale_warning.onsite_access}} daily on-site accesses</b> during this stage. Please reduce the value below the threshold of {{scale_warning.threshold}}. See with your supervisor, DSO/LEXGLIMOS and space manager if this scenario is acceptable and if any additional measures are required.
|
||||
</div>
|
||||
{% else %}
|
||||
<p><b>Note:</b> The CERN COVID Level is not specified.</p>
|
||||
|
|
@ -67,13 +67,13 @@
|
|||
{% endblock report_summary %}
|
||||
|
||||
{% block report_summary_footnote %}
|
||||
{% if ((prob_inf > 15) or (expected_new_cases >= 1)) %}
|
||||
{% if ((prob_inf > 10) or (expected_new_cases >= 1)) %}
|
||||
This exceeds the authorised risk threshold or number of expected new cases.
|
||||
The risk level must be reduced before this activity can be undertaken.
|
||||
{% elif (5 <= prob_inf <= 15) %}
|
||||
{% elif (2 <= prob_inf <= 10) %}
|
||||
This activity has an elevated level of risk, ALARA principles must be applied to minimise the level of risk before undertaking the activity.
|
||||
See the footnotes for more details on the ALARA principles.
|
||||
{% elif (prob_inf < 5) %}
|
||||
{% elif (prob_inf < 2) %}
|
||||
This level of risk is within acceptable parameters, no further actions are required.
|
||||
{% endif %}
|
||||
{% endblock report_summary_footnote %}
|
||||
|
|
@ -124,9 +124,9 @@
|
|||
<div class="d-flex">
|
||||
<div class="split">
|
||||
<div class="d-flex flex-column col-xl-8 p-0">
|
||||
<div class="alert alert-success" role="alert">Events with a <strong>P(i) less than 5%</strong> may go ahead without further mitigation measures.</div>
|
||||
<div class="alert alert-warning" role="alert">Events with a <strong>P(i) between 5% and 15%</strong> shall be subject to ALARA principles (see footnote) to minimise the risk before proceeding.</div>
|
||||
<div class="alert alert-danger mb-0" role="alert">Events with a <strong>P(i) exceeding 15% or a number of expected new cases that exceeds 1</strong> may not take place until additional measures are in place and a risk reduction has been performed.</div>
|
||||
<div class="alert alert-success" role="alert">Events with a <strong>P(i) less than 2%</strong> may go ahead without further mitigation measures.</div>
|
||||
<div class="alert alert-warning" role="alert">Events with a <strong>P(i) between 2% and 10%</strong> shall be subject to ALARA principles (see footnote) to minimise the risk before proceeding.</div>
|
||||
<div class="alert alert-danger mb-0" role="alert">Events with a <strong>P(i) exceeding 10% or a number of expected new cases that exceeds 1</strong> may not take place until additional measures are in place and a risk reduction has been performed.</div>
|
||||
</div>
|
||||
<div class="col-xl-3 align-self-center text-center"><img id="cern_level" class="rounded" src="{{ calculator_prefix }}/static/images/warning_scale/{{scale_warning.level}}.png"></div>
|
||||
</div>
|
||||
|
|
|
|||
|
|
@ -22,10 +22,6 @@ geneva_coordinates = (46.204391, 6.143158)
|
|||
local_hourly_temperatures_celsius_per_hour = get_hourly_temperatures_celsius_per_hour(
|
||||
geneva_coordinates)
|
||||
|
||||
# Load the weather data (temperature in kelvin) for Toronto.
|
||||
toronto_coordinates = (43.667, 79.400)
|
||||
toronto_hourly_temperatures_celsius_per_hour = get_hourly_temperatures_celsius_per_hour(
|
||||
toronto_coordinates)
|
||||
|
||||
# Geneva hourly temperatures as piecewise constant function (in Kelvin).
|
||||
GenevaTemperatures_hourly = {
|
||||
|
|
@ -38,25 +34,9 @@ GenevaTemperatures_hourly = {
|
|||
for month, temperatures in local_hourly_temperatures_celsius_per_hour.items()
|
||||
}
|
||||
|
||||
# Toronto hourly temperatures as piecewise constant function (in Kelvin).
|
||||
TorontoTemperatures_hourly = {
|
||||
month: models.PiecewiseConstant(
|
||||
# NOTE: It is important that the time type is float, not np.float, in
|
||||
# order to allow hashability (for caching).
|
||||
tuple(float(time) for time in range(25)),
|
||||
tuple(273.15 + np.array(temperatures)),
|
||||
)
|
||||
for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items()
|
||||
}
|
||||
|
||||
# Same Geneva temperatures on a finer temperature mesh (every 6 minutes).
|
||||
GenevaTemperatures = {
|
||||
month: GenevaTemperatures_hourly[month].refine(refine_factor=10)
|
||||
for month, temperatures in local_hourly_temperatures_celsius_per_hour.items()
|
||||
}
|
||||
|
||||
# Same Toronto temperatures on a finer temperature mesh (every 6 minutes).
|
||||
TorontoTemperatures = {
|
||||
month: TorontoTemperatures_hourly[month].refine(refine_factor=10)
|
||||
for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items()
|
||||
}
|
||||
|
|
|
|||
306
cara/models.py
306
cara/models.py
|
|
@ -51,7 +51,6 @@ from .utils import method_cache
|
|||
from .dataclass_utils import nested_replace
|
||||
|
||||
oneoverln2 = 1 / np.log(2)
|
||||
|
||||
# Define types for items supporting vectorisation. In the future this may be replaced
|
||||
# by ``np.ndarray[<type>]`` once/if that syntax is supported. Note that vectorization
|
||||
# implies 1d arrays: multi-dimensional arrays are not supported.
|
||||
|
|
@ -152,9 +151,9 @@ class PiecewiseConstant:
|
|||
|
||||
def __post_init__(self):
|
||||
if len(self.transition_times) != len(self.values)+1:
|
||||
raise ValueError("transition_times should contain one more element than values")
|
||||
raise ValueError("transition_times must 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")
|
||||
raise ValueError("transition_times must not contain duplicated elements and must 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")
|
||||
|
|
@ -460,8 +459,8 @@ class SARSCoV2(Virus):
|
|||
piecewise constant model (for more details see A. Henriques et al,
|
||||
CERN-OPEN-2021-004, DOI: 10.17181/CERN.1GDQ.5Y75)
|
||||
"""
|
||||
# Taken from Morris et al (https://doi.org/10.7554/eLife.65902) data at T = 22°C and RH = 40 % and
|
||||
# from Doremalen et al (https://www.nejm.org/doi/10.1056/NEJMc2004973).
|
||||
# Taken from Morris et al (https://doi.org/10.7554/eLife.65902) data at T = 22°C and RH = 40 %,
|
||||
# and from Doremalen et al (https://www.nejm.org/doi/10.1056/NEJMc2004973).
|
||||
return np.piecewise(humidity, [humidity <= 0.4, humidity > 0.4], [6.43, 1.1])
|
||||
|
||||
|
||||
|
|
@ -561,6 +560,58 @@ Mask.types = {
|
|||
}
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Particle:
|
||||
"""
|
||||
Represents an aerosol particle.
|
||||
"""
|
||||
|
||||
#: diameter of the aerosol in microns
|
||||
diameter: typing.Union[None,_VectorisedFloat] = None
|
||||
|
||||
def settling_velocity(self, evaporation_factor: float=0.3) -> _VectorisedFloat:
|
||||
"""
|
||||
Settling velocity (i.e. speed of deposition on the floor due
|
||||
to gravity), for aerosols, in m/s. Diameter-dependent expression
|
||||
from https://doi.org/10.1101/2021.10.14.21264988
|
||||
When an aerosol-diameter is not given, returns
|
||||
the default value of 1.88e-4 m/s (corresponds to diameter of
|
||||
2.5 microns, i.e. geometric average of the breathing
|
||||
expiration distribution, taking evaporation into account, see
|
||||
https://doi.org/10.1101/2021.10.14.21264988)
|
||||
evaporation_factor represents the factor applied to the diameter,
|
||||
due to instantaneous evaporation of the particle in the air.
|
||||
"""
|
||||
if self.diameter is None:
|
||||
vg = 1.88e-4
|
||||
else:
|
||||
vg = 1.88e-4 * (self.diameter*evaporation_factor / 2.5)**2
|
||||
return vg
|
||||
|
||||
def fraction_deposited(self, evaporation_factor: float=0.3) -> _VectorisedFloat:
|
||||
"""
|
||||
The fraction of particles actually deposited in the respiratory
|
||||
tract (over the total number of particles). It depends on the
|
||||
particle diameter.
|
||||
From W. C. Hinds, New York, Wiley, 1999 (pp. 233 – 259).
|
||||
evaporation_factor represents the factor applied to the diameter,
|
||||
due to instantaneous evaporation of the particle in the air.
|
||||
"""
|
||||
if self.diameter is None:
|
||||
# model is not evaluated for specific values of aerosol
|
||||
# diameters - we choose a single "average" deposition factor,
|
||||
# as in https://doi.org/10.1101/2021.10.14.21264988.
|
||||
fdep = 0.6
|
||||
else:
|
||||
# deposition fraction depends on aerosol particle diameter.
|
||||
d = (self.diameter * evaporation_factor)
|
||||
IFrac = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8)))))
|
||||
fdep = IFrac * (0.0587
|
||||
+ (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d))))
|
||||
+ (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d)))))
|
||||
return fdep
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class _ExpirationBase:
|
||||
"""
|
||||
|
|
@ -570,6 +621,13 @@ class _ExpirationBase:
|
|||
#: Pre-populated examples of Expirations.
|
||||
types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
"""
|
||||
the Particle object representing the aerosol - here the default one
|
||||
"""
|
||||
return Particle()
|
||||
|
||||
def aerosols(self, mask: Mask):
|
||||
"""
|
||||
total volume of aerosols expired per volume of air (mL/cm^3).
|
||||
|
|
@ -587,9 +645,18 @@ class Expiration(_ExpirationBase):
|
|||
#: diameter of the aerosol in microns
|
||||
diameter: _VectorisedFloat
|
||||
|
||||
#: total concentration of aerosols (cm^-3)
|
||||
#: total concentration of aerosols per unit volume of expired air
|
||||
# (in cm^-3), integrated over all aerosol diameters (corresponding
|
||||
# to c_n,i in Eq. (4) of https://doi.org/10.1101/2021.10.14.21264988)
|
||||
cn: float = 1.
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
"""
|
||||
the Particle object representing the aerosol
|
||||
"""
|
||||
return Particle(diameter=self.diameter)
|
||||
|
||||
@cached()
|
||||
def aerosols(self, mask: Mask):
|
||||
""" Result is in mL.cm^-3 """
|
||||
|
|
@ -608,18 +675,18 @@ class MultipleExpiration(_ExpirationBase):
|
|||
Group together different modes of expiration, that represent
|
||||
each the main expiration mode for a certain fraction of time (given by
|
||||
the weights).
|
||||
This class can only be used with single diameters (it cannot
|
||||
be used with diameter distributions).
|
||||
This class can only be used with single diameters defined in each
|
||||
expiration (it cannot be used with diameter distributions).
|
||||
"""
|
||||
expirations: typing.Tuple[_ExpirationBase, ...]
|
||||
weights: typing.Tuple[float, ...]
|
||||
|
||||
def __post_init__(self):
|
||||
if len(self.expirations) != len(self.weights):
|
||||
raise ValueError("expirations and weigths should contain the"
|
||||
raise ValueError("expirations and weigths must contain the"
|
||||
"same number of elements")
|
||||
if not all(np.isscalar(e.diameter) for e in self.expirations):
|
||||
raise ValueError("diameters should all be scalars")
|
||||
raise ValueError("diameters must all be scalars")
|
||||
|
||||
def aerosols(self, mask: Mask):
|
||||
return np.array([
|
||||
|
|
@ -709,13 +776,28 @@ class _PopulationWithVirus(Population):
|
|||
"""
|
||||
return 1.
|
||||
|
||||
def aerosols(self):
|
||||
"""
|
||||
total volume of aerosols expired per volume of air (mL/cm^3).
|
||||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
|
||||
def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate of virions per fraction of aerosol volume in
|
||||
the expired air, if the infected population is present, in cm^3/(mL.h).
|
||||
It should not be a function of time.
|
||||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
|
||||
@method_cache
|
||||
def emission_rate_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate if the infected population is present
|
||||
(in virions / h). It should not be a function of time.
|
||||
(in virions / h).
|
||||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
return (self.emission_rate_per_aerosol_when_present() *
|
||||
self.aerosols())
|
||||
|
||||
def emission_rate(self, time) -> _VectorisedFloat:
|
||||
"""
|
||||
|
|
@ -734,16 +816,33 @@ class _PopulationWithVirus(Population):
|
|||
|
||||
return self.emission_rate_when_present()
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
"""
|
||||
the Particle object representing the aerosol expired by the
|
||||
population - here we take the default Particle object
|
||||
"""
|
||||
return Particle()
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class EmittingPopulation(_PopulationWithVirus):
|
||||
#: The emission rate of a single individual, in virions / h.
|
||||
known_individual_emission_rate: float
|
||||
|
||||
@method_cache
|
||||
def emission_rate_when_present(self) -> _VectorisedFloat:
|
||||
def aerosols(self):
|
||||
"""
|
||||
The emission rate if the infected population is present.
|
||||
total volume of aerosols expired per volume of air (mL/cm^3).
|
||||
Here arbitrarily set to 1 as the full emission rate is known.
|
||||
"""
|
||||
return 1.
|
||||
|
||||
@method_cache
|
||||
def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate of virions per fraction of aerosol volume in
|
||||
the expired air, if the infected population is present, in cm^3/(mL.h).
|
||||
It should not be a function of time.
|
||||
"""
|
||||
return self.known_individual_emission_rate * self.number
|
||||
|
||||
|
|
@ -757,29 +856,37 @@ class InfectedPopulation(_PopulationWithVirus):
|
|||
def fraction_of_infectious_virus(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The fraction of infectious virus.
|
||||
|
||||
"""
|
||||
return self.virus.viable_to_RNA_ratio * (1 - self.host_immunity)
|
||||
|
||||
@method_cache
|
||||
def emission_rate_when_present(self) -> _VectorisedFloat:
|
||||
def aerosols(self):
|
||||
"""
|
||||
The emission rate if the infected population is present.
|
||||
Note that the rate is not currently time-dependent.
|
||||
total volume of aerosols expired per volume of air (mL/cm^3).
|
||||
"""
|
||||
# Emission Rate (virions / h)
|
||||
# Note on units: exhalation rate is in m^3/h, aerosols in mL/cm^3
|
||||
# and viral load in virus/mL -> 1e6 conversion factor
|
||||
return self.expiration.aerosols(self.mask)
|
||||
|
||||
aerosols = self.expiration.aerosols(self.mask)
|
||||
@method_cache
|
||||
def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate of virions per fraction of aerosol volume in
|
||||
the expired air, if the infected population is present, in cm^3/(mL.h).
|
||||
It should not be a function of time.
|
||||
"""
|
||||
# Note on units: exhalation rate is in m^3/h -> 1e6 conversion factor
|
||||
|
||||
ER = (self.virus.viral_load_in_sputum *
|
||||
self.activity.exhalation_rate *
|
||||
10 ** 6 *
|
||||
aerosols)
|
||||
10 ** 6)
|
||||
|
||||
return ER * self.number
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
"""
|
||||
the Particle object representing the aerosol - here the default one
|
||||
"""
|
||||
return self.expiration.particle
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class ConcentrationModel:
|
||||
|
|
@ -798,20 +905,7 @@ class ConcentrationModel:
|
|||
|
||||
def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat:
|
||||
# Equilibrium velocity of particle motion toward the floor
|
||||
if (isinstance(self.infected, InfectedPopulation)
|
||||
and isinstance(self.infected.expiration, Expiration)):
|
||||
d = self.infected.expiration.diameter * self.evaporation_factor
|
||||
vg = 1.88e-4 * (d / 2.5)**2
|
||||
# see https://doi.org/10.1101/2021.10.14.21264988
|
||||
# (velocity of 1.88e-4 corresponds to diameter of 2.5 microns)
|
||||
else:
|
||||
# model is not evaluated for specific values of aerosol
|
||||
# diameters - we choose a single velocity value
|
||||
# corresponding to that obtained with a diameter of 2.5 microns
|
||||
# (geometric average of the breathing expiration distribution,
|
||||
# taking evaporation into account, see
|
||||
# https://doi.org/10.1101/2021.10.14.21264988)
|
||||
vg = 1.88e-4
|
||||
vg = self.infected.particle.settling_velocity(self.evaporation_factor)
|
||||
# Height of the emission source to the floor - i.e. mouth/nose (m)
|
||||
h = 1.5
|
||||
# Deposition rate (h^-1)
|
||||
|
|
@ -989,26 +1083,14 @@ class ExposureModel:
|
|||
#: The number of times the exposure event is repeated (default 1).
|
||||
repeats: int = 1
|
||||
|
||||
def fraction_deposited(self):
|
||||
def fraction_deposited(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The fraction of viruses actually deposited in the respiratory tract.
|
||||
From W. C. Hinds, New York, Wiley, 1999 (pp. 233 – 259).
|
||||
The fraction of particles actually deposited in the respiratory
|
||||
tract (over the total number of particles). It depends on the
|
||||
particle diameter.
|
||||
"""
|
||||
if not (isinstance(self.concentration_model.infected,InfectedPopulation)
|
||||
and isinstance(self.concentration_model.infected.expiration,Expiration)):
|
||||
# model is not evaluated for specific values of aerosol
|
||||
# diameters - we choose a single "average" deposition factor,
|
||||
# as in https://doi.org/10.1101/2021.10.14.21264988.
|
||||
fdep = 0.6
|
||||
else:
|
||||
# deposition factor depends on aerosol particle diameter.
|
||||
d = (self.concentration_model.infected.expiration.diameter
|
||||
* self.concentration_model.evaporation_factor)
|
||||
IFrac = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8)))))
|
||||
fdep = IFrac * (0.0587
|
||||
+ (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d))))
|
||||
+ (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d)))))
|
||||
return fdep
|
||||
return self.concentration_model.infected.particle.fraction_deposited(
|
||||
self.concentration_model.evaporation_factor)
|
||||
|
||||
def _normed_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
|
||||
"""The number of virions per meter^3 between any two times, normalized
|
||||
|
|
@ -1028,11 +1110,6 @@ class ExposureModel:
|
|||
elif time1 <= start and stop < time2:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(start, stop)
|
||||
return exposure
|
||||
|
||||
def exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
|
||||
"""The number of virions per meter^3 between any two times."""
|
||||
return (self._normed_exposure_between_bounds(time1, time2) *
|
||||
self.concentration_model.infected.emission_rate_when_present())
|
||||
|
||||
def _normed_exposure(self) -> _VectorisedFloat:
|
||||
"""
|
||||
|
|
@ -1046,76 +1123,57 @@ class ExposureModel:
|
|||
|
||||
return normed_exposure * self.repeats
|
||||
|
||||
def exposure(self) -> _VectorisedFloat:
|
||||
def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
|
||||
"""
|
||||
The number of virus per meter^3. With sampled diameters, the
|
||||
aerosol volume has to be put back in the exposure before taking
|
||||
the mean, to obtain the proper result for the exposure (which
|
||||
corresponds to an integration on diameters).
|
||||
The number of virus per m^3 deposited on the respiratory tract
|
||||
between any two times.
|
||||
"""
|
||||
emission_rate = self.concentration_model.infected.emission_rate_when_present()
|
||||
if (not isinstance(self.concentration_model.infected,InfectedPopulation)
|
||||
or not isinstance(self.concentration_model.infected.expiration,Expiration)
|
||||
or np.isscalar(self.concentration_model.infected.expiration.diameter)
|
||||
):
|
||||
# in all these cases, there is no distribution of
|
||||
# diameters that need to be integrated over
|
||||
return self._normed_exposure() * emission_rate
|
||||
else:
|
||||
# the mean of the diameter-dependent exposure (including
|
||||
# aerosols volume, but NOT the other factors) has to be
|
||||
# taken first (this is equivalent to integrating over the
|
||||
# diameters)
|
||||
mask = self.concentration_model.infected.mask
|
||||
aerosols = self.concentration_model.infected.expiration.aerosols(mask)
|
||||
return (np.array(self._normed_exposure()*aerosols).mean() *
|
||||
emission_rate/aerosols)
|
||||
|
||||
def _deposited_exposure(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The number of virus per meter^3 deposited on the respiratory
|
||||
tract. As in the exposure method, with sampled diameters the
|
||||
aerosol volume and the deposited fraction, have to be put
|
||||
in the deposited exposure before taking the mean, to obtain the
|
||||
proper result (which corresponds to an integration on diameters).
|
||||
"""
|
||||
emission_rate = self.concentration_model.infected.emission_rate_when_present()
|
||||
if (not isinstance(self.concentration_model.infected,InfectedPopulation)
|
||||
or not isinstance(self.concentration_model.infected.expiration,Expiration)
|
||||
or np.isscalar(self.concentration_model.infected.expiration.diameter)
|
||||
):
|
||||
# in all these cases, there is no distribution of
|
||||
# diameters that need to be integrated over
|
||||
return (self._normed_exposure() *
|
||||
self.fraction_deposited() *
|
||||
emission_rate)
|
||||
else:
|
||||
# the mean of the diameter-dependent exposure (including
|
||||
# aerosols volume, but NOT the other factors) has to be
|
||||
# taken first (this is equivalent to integrating over the
|
||||
# diameters)
|
||||
mask = self.concentration_model.infected.mask
|
||||
aerosols = self.concentration_model.infected.expiration.aerosols(mask)
|
||||
return (np.array(self._normed_exposure() * aerosols *
|
||||
self.fraction_deposited()).mean() *
|
||||
emission_rate/aerosols)
|
||||
|
||||
def infection_probability(self) -> _VectorisedFloat:
|
||||
deposited_exposure = self._deposited_exposure()
|
||||
|
||||
emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present()
|
||||
aerosols = self.concentration_model.infected.aerosols()
|
||||
fdep = self.fraction_deposited()
|
||||
f_inf = self.concentration_model.infected.fraction_of_infectious_virus()
|
||||
|
||||
inf_aero = (
|
||||
self.exposed.activity.inhalation_rate *
|
||||
(1 - self.exposed.mask.inhale_efficiency()) *
|
||||
deposited_exposure * f_inf
|
||||
)
|
||||
diameter = self.concentration_model.infected.particle.diameter
|
||||
|
||||
if not np.isscalar(diameter) and diameter is not None:
|
||||
# we compute first the mean of all diameter-dependent quantities
|
||||
# to perform properly the Monte-Carlo integration over
|
||||
# particle diameters (doing things in another order would
|
||||
# lead to wrong results).
|
||||
dep_exposure_integrated = np.array(self._normed_exposure_between_bounds(time1, time2) *
|
||||
aerosols *
|
||||
fdep).mean()
|
||||
else:
|
||||
# in the case of a single diameter or no diameter defined,
|
||||
# one should not take any mean at this stage.
|
||||
dep_exposure_integrated = self._normed_exposure_between_bounds(time1, time2)*aerosols*fdep
|
||||
|
||||
# then we multiply by the diameter-independent quantity emission_rate_per_aerosol,
|
||||
# and parameters of the vD equation (i.e. f_inf, BR_k and n_in).
|
||||
return (dep_exposure_integrated * emission_rate_per_aerosol *
|
||||
f_inf * self.exposed.activity.inhalation_rate *
|
||||
(1 - self.exposed.mask.inhale_efficiency()))
|
||||
|
||||
def deposited_exposure(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The number of virus per m^3 deposited on the respiratory tract.
|
||||
"""
|
||||
deposited_exposure = 0.0
|
||||
|
||||
for start, stop in self.exposed.presence.boundaries():
|
||||
deposited_exposure += self.deposited_exposure_between_bounds(start, stop)
|
||||
|
||||
return deposited_exposure * self.repeats
|
||||
|
||||
def infection_probability(self) -> _VectorisedFloat:
|
||||
# viral dose (vD)
|
||||
vD = self.deposited_exposure()
|
||||
|
||||
# oneoverln2 multiplied by ID_50 corresponds to ID_63.
|
||||
infectious_dose = oneoverln2 * self.concentration_model.virus.infectious_dose
|
||||
|
||||
# Probability of infection.
|
||||
return (1 - np.exp(-((inf_aero * (1 - self.exposed.host_immunity))/(infectious_dose *
|
||||
return (1 - np.exp(-((vD * (1 - self.exposed.host_immunity))/(infectious_dose *
|
||||
self.concentration_model.virus.transmissibility_factor)))) * 100
|
||||
|
||||
def expected_new_cases(self) -> _VectorisedFloat:
|
||||
|
|
|
|||
|
|
@ -165,15 +165,15 @@ def expiration_distribution(BLO_factors):
|
|||
"""
|
||||
Returns an Expiration with an aerosol diameter distribution, defined
|
||||
by the BLO factors (a length-3 tuple).
|
||||
The total concentration of aerosols is computed by integrating
|
||||
the distribution between 0.1 and D microns - these boundaries are
|
||||
The total concentration of aerosols, cn, is computed by integrating
|
||||
the distribution between 0.1 and 30 microns - these boundaries are
|
||||
an historical choice based on previous implementations of the model
|
||||
(it limits the influence of the O-mode).
|
||||
"""
|
||||
dscan = np.linspace(0.1, 30. ,3000)
|
||||
return mc.Expiration(CustomKernel(dscan,
|
||||
BLOmodel(BLO_factors).distribution(dscan),kernel_bandwidth=0.1),
|
||||
BLOmodel(BLO_factors).integrate(0.1, 30.))
|
||||
cn=BLOmodel(BLO_factors).integrate(0.1, 30.))
|
||||
|
||||
|
||||
def dilution_factor(distance, D=0.02):
|
||||
|
|
|
|||
|
|
@ -13,8 +13,6 @@ from cara.monte_carlo.data import expiration_distributions
|
|||
|
||||
# TODO: seed better the random number generators
|
||||
np.random.seed(2000)
|
||||
SAMPLE_SIZE = 250000
|
||||
TOLERANCE = 0.02
|
||||
|
||||
def test_model_from_dict(baseline_form_data):
|
||||
form = model_generator.FormData.from_dict(baseline_form_data)
|
||||
|
|
@ -35,6 +33,8 @@ def test_model_from_dict_invalid(baseline_form_data):
|
|||
]
|
||||
)
|
||||
def test_blend_expiration(mask_type):
|
||||
SAMPLE_SIZE = 250000
|
||||
TOLERANCE = 0.02
|
||||
blend = {'Breathing': 2, 'Speaking': 1}
|
||||
r = model_generator.build_expiration(blend).build_model(SAMPLE_SIZE)
|
||||
mask = models.Mask.types[mask_type]
|
||||
|
|
|
|||
|
|
@ -77,25 +77,25 @@ def known_concentrations(func):
|
|||
@pytest.mark.parametrize(
|
||||
"population, cm, expected_exposure, expected_probability", [
|
||||
[populations[1], known_concentrations(lambda t: 36.),
|
||||
np.array([432, 432]), np.array([67.9503762594, 65.2366759251])],
|
||||
np.array([64.02320633, 59.45012016]), np.array([67.9503762594, 65.2366759251])],
|
||||
|
||||
[populations[2], known_concentrations(lambda t: 36.),
|
||||
np.array([432, 432]), np.array([51.6749232285, 55.6374622042])],
|
||||
np.array([40.91708675, 45.73086166]), np.array([51.6749232285, 55.6374622042])],
|
||||
|
||||
[populations[0], known_concentrations(lambda t: np.array([36., 72.])),
|
||||
np.array([432, 864]), np.array([55.6374622042, 80.3196524031])],
|
||||
np.array([45.73086166, 91.46172332]), np.array([55.6374622042, 80.3196524031])],
|
||||
|
||||
[populations[1], known_concentrations(lambda t: np.array([36., 72.])),
|
||||
np.array([432, 864]), np.array([67.9503762594, 87.9151129926])],
|
||||
np.array([64.02320633, 118.90024032]), np.array([67.9503762594, 87.9151129926])],
|
||||
|
||||
[populations[2], known_concentrations(lambda t: np.array([36., 72.])),
|
||||
np.array([432, 864]), np.array([51.6749232285, 80.3196524031])],
|
||||
np.array([40.91708675, 91.46172332]), np.array([51.6749232285, 80.3196524031])],
|
||||
])
|
||||
def test_exposure_model_ndarray(population, cm,
|
||||
expected_exposure, expected_probability):
|
||||
model = ExposureModel(cm, population)
|
||||
np.testing.assert_almost_equal(
|
||||
model.exposure(), expected_exposure
|
||||
model.deposited_exposure(), expected_exposure
|
||||
)
|
||||
np.testing.assert_almost_equal(
|
||||
model.infection_probability(), expected_probability, decimal=10
|
||||
|
|
@ -107,33 +107,43 @@ def test_exposure_model_ndarray(population, cm,
|
|||
assert model.expected_new_cases().shape == (2,)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("population", populations)
|
||||
def test_exposure_model_ndarray_and_float_mix(population):
|
||||
@pytest.mark.parametrize("population, expected_deposited_exposure", [
|
||||
[populations[0], np.array([1.52436206, 1.52436206])],
|
||||
[populations[1], np.array([2.13410688, 1.98167067])],
|
||||
[populations[2], np.array([1.36390289, 1.52436206])],
|
||||
])
|
||||
def test_exposure_model_ndarray_and_float_mix(population, expected_deposited_exposure):
|
||||
cm = known_concentrations(
|
||||
lambda t: 0. if np.floor(t) % 2 else np.array([1.2, 1.2]))
|
||||
model = ExposureModel(cm, population)
|
||||
|
||||
expected_exposure = np.array([14.4, 14.4])
|
||||
np.testing.assert_almost_equal(
|
||||
model.exposure(), expected_exposure
|
||||
model.deposited_exposure(), expected_deposited_exposure
|
||||
)
|
||||
|
||||
assert isinstance(model.infection_probability(), np.ndarray)
|
||||
assert isinstance(model.expected_new_cases(), np.ndarray)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("population", populations)
|
||||
def test_exposure_model_compare_scalar_vector(population):
|
||||
cm_scalar = known_concentrations(lambda t: 1.2)
|
||||
@pytest.mark.parametrize("population, expected_deposited_exposure", [
|
||||
[populations[0], np.array([1.52436206, 1.52436206])],
|
||||
[populations[1], np.array([2.13410688, 1.98167067])],
|
||||
[populations[2], np.array([1.36390289, 1.52436206])],
|
||||
])
|
||||
def test_exposure_model_vector(population, expected_deposited_exposure):
|
||||
cm_array = known_concentrations(lambda t: np.array([1.2, 1.2]))
|
||||
model_scalar = ExposureModel(cm_scalar, population)
|
||||
model_array = ExposureModel(cm_array, population)
|
||||
expected_exposure = 14.4
|
||||
np.testing.assert_almost_equal(
|
||||
model_scalar.exposure(), expected_exposure
|
||||
model_array.deposited_exposure(), np.array(expected_deposited_exposure)
|
||||
)
|
||||
|
||||
|
||||
def test_exposure_model_scalar():
|
||||
cm_scalar = known_concentrations(lambda t: 1.2)
|
||||
model_scalar = ExposureModel(cm_scalar, populations[0])
|
||||
expected_deposited_exposure = 1.52436206
|
||||
np.testing.assert_almost_equal(
|
||||
model_array.exposure(), np.array([expected_exposure]*2)
|
||||
model_scalar.deposited_exposure(), expected_deposited_exposure
|
||||
)
|
||||
|
||||
|
||||
|
|
@ -163,28 +173,28 @@ def conc_model():
|
|||
)
|
||||
|
||||
|
||||
# Expected exposure were computed with a trapezoidal integration, using
|
||||
# Expected deposited exposure were computed with a trapezoidal integration, using
|
||||
# a mesh of 10'000 pts per exposed presence interval.
|
||||
@pytest.mark.parametrize(
|
||||
["exposed_time_interval", "expected_exposure"],
|
||||
["exposed_time_interval", "expected_deposited_exposure"],
|
||||
[
|
||||
[(0., 1.), 266.67176],
|
||||
[(1., 1.01), 3.0879539],
|
||||
[(1.01, 1.02), 3.00082435],
|
||||
[(12., 12.01), 0.095063235],
|
||||
[(12., 24.), 3775.65025],
|
||||
[(0., 24.), 4097.8494],
|
||||
[(0., 1.), 45.6008710],
|
||||
[(1., 1.01), 0.5280401],
|
||||
[(1.01, 1.02), 0.51314096385],
|
||||
[(12., 12.01), 0.016255813185],
|
||||
[(12., 24.), 645.63619275],
|
||||
[(0., 24.), 700.7322474],
|
||||
]
|
||||
)
|
||||
def test_exposure_model_integral_accuracy(exposed_time_interval,
|
||||
expected_exposure, conc_model):
|
||||
expected_deposited_exposure, conc_model):
|
||||
presence_interval = models.SpecificInterval((exposed_time_interval,))
|
||||
population = models.Population(
|
||||
10, presence_interval, models.Mask.types['Type I'],
|
||||
models.Activity.types['Standing'], 0.,
|
||||
)
|
||||
model = ExposureModel(conc_model, population)
|
||||
np.testing.assert_allclose(model.exposure(), expected_exposure)
|
||||
np.testing.assert_allclose(model.deposited_exposure(), expected_deposited_exposure)
|
||||
|
||||
|
||||
def test_infectious_dose_vectorisation():
|
||||
|
|
@ -212,7 +222,7 @@ def test_infectious_dose_vectorisation():
|
|||
10, presence_interval, models.Mask.types['Type I'],
|
||||
models.Activity.types['Standing'], 0.,
|
||||
)
|
||||
model = ExposureModel(cm, population) #, fraction_deposited=1.0
|
||||
model = ExposureModel(cm, population)
|
||||
inf_probability = model.infection_probability()
|
||||
assert isinstance(inf_probability, np.ndarray)
|
||||
assert inf_probability.shape == (3, )
|
||||
|
|
|
|||
|
|
@ -13,7 +13,7 @@ def test_multiple_wrong_weight_size():
|
|||
e_base = models.Expiration(2.5)
|
||||
with pytest.raises(
|
||||
ValueError,
|
||||
match=re.escape("expirations and weigths should contain the"
|
||||
match=re.escape("expirations and weigths must contain the"
|
||||
"same number of elements")
|
||||
):
|
||||
e = models.MultipleExpiration([e_base, e_base], weights)
|
||||
|
|
@ -26,7 +26,7 @@ def test_multiple_wrong_diameters():
|
|||
e3 = models.Expiration(2.)
|
||||
with pytest.raises(
|
||||
ValueError,
|
||||
match=re.escape("diameters should all be scalars")
|
||||
match=re.escape("diameters must all be scalars")
|
||||
):
|
||||
e = models.MultipleExpiration([e1, e2, e3], weights)
|
||||
|
||||
|
|
|
|||
|
|
@ -97,7 +97,7 @@ def test_r0(baseline_exposure_model):
|
|||
# expected r0 was computed with a trapezoidal integration, using
|
||||
# a mesh of 100'000 pts per exposed presence interval.
|
||||
r0 = baseline_exposure_model.reproduction_number()
|
||||
npt.assert_allclose(r0, 776.9419902161412)
|
||||
npt.assert_allclose(r0, 776.941990)
|
||||
|
||||
|
||||
def test_periodic_window(baseline_periodic_window, baseline_room):
|
||||
|
|
@ -389,16 +389,16 @@ def build_exposure_model(concentration_model):
|
|||
)
|
||||
|
||||
|
||||
# expected exposure were computed with a trapezoidal integration, using
|
||||
# expected deposited exposure were computed with a trapezoidal integration, using
|
||||
# a mesh of 100'000 pts per exposed presence interval.
|
||||
@pytest.mark.parametrize(
|
||||
"month, expected_exposure",
|
||||
"month, expected_deposited_exposure",
|
||||
[
|
||||
['Jan', 503.254087759],
|
||||
['Jun', 2294.71115639],
|
||||
['Jan', 377.440565819],
|
||||
['Jun', 1721.03336729],
|
||||
],
|
||||
)
|
||||
def test_exposure_hourly_dep(month,expected_exposure):
|
||||
def test_exposure_hourly_dep(month,expected_deposited_exposure):
|
||||
m = build_exposure_model(
|
||||
build_hourly_dependent_model(
|
||||
month,
|
||||
|
|
@ -406,20 +406,20 @@ def test_exposure_hourly_dep(month,expected_exposure):
|
|||
intervals_presence_infected=((8., 12.), (13., 17.))
|
||||
)
|
||||
)
|
||||
exposure = m.exposure()
|
||||
npt.assert_allclose(exposure, expected_exposure)
|
||||
deposited_exposure = m.deposited_exposure()
|
||||
npt.assert_allclose(deposited_exposure, expected_deposited_exposure)
|
||||
|
||||
# expected exposure were computed with a trapezoidal integration, using
|
||||
# expected deposited exposure were computed with a trapezoidal integration, using
|
||||
# a mesh of 100'000 pts per exposed presence interval and 25 pts per hour
|
||||
# for the temperature discretization.
|
||||
@pytest.mark.parametrize(
|
||||
"month, expected_exposure",
|
||||
"month, expected_deposited_exposure",
|
||||
[
|
||||
['Jan', 511.118941481],
|
||||
['Jun', 2398.90129579],
|
||||
['Jan', 383.339206111],
|
||||
['Jun', 1799.17597184],
|
||||
],
|
||||
)
|
||||
def test_exposure_hourly_dep_refined(month,expected_exposure):
|
||||
def test_exposure_hourly_dep_refined(month,expected_deposited_exposure):
|
||||
m = build_exposure_model(
|
||||
build_hourly_dependent_model(
|
||||
month,
|
||||
|
|
@ -428,5 +428,5 @@ def test_exposure_hourly_dep_refined(month,expected_exposure):
|
|||
temperatures=data.GenevaTemperatures,
|
||||
)
|
||||
)
|
||||
exposure = m.exposure()
|
||||
npt.assert_allclose(exposure, expected_exposure, rtol=0.02)
|
||||
deposited_exposure = m.deposited_exposure()
|
||||
npt.assert_allclose(deposited_exposure, expected_deposited_exposure, rtol=0.02)
|
||||
|
|
|
|||
|
|
@ -89,6 +89,6 @@ def test_build_concentration_model(baseline_mc_model: cara.monte_carlo.Concentra
|
|||
def test_build_exposure_model(baseline_mc_exposure_model: cara.monte_carlo.ExposureModel):
|
||||
model = baseline_mc_exposure_model.build_model(7)
|
||||
assert isinstance(model, cara.models.ExposureModel)
|
||||
prob = model.exposure()
|
||||
prob = model.deposited_exposure()
|
||||
assert isinstance(prob, np.ndarray)
|
||||
assert prob.shape == (7, )
|
||||
|
|
|
|||
|
|
@ -12,6 +12,31 @@ np.random.seed(2000)
|
|||
SAMPLE_SIZE = 250000
|
||||
TOLERANCE = 0.05
|
||||
|
||||
# Load the weather data (temperature in kelvin) for Toronto.
|
||||
toronto_coordinates = (43.667, 79.400)
|
||||
toronto_hourly_temperatures_celsius_per_hour = data.get_hourly_temperatures_celsius_per_hour(
|
||||
toronto_coordinates)
|
||||
|
||||
|
||||
# Toronto hourly temperatures as piecewise constant function (in Kelvin).
|
||||
TorontoTemperatures_hourly = {
|
||||
month: models.PiecewiseConstant(
|
||||
# NOTE: It is important that the time type is float, not np.float, in
|
||||
# order to allow hashability (for caching).
|
||||
tuple(float(time) for time in range(25)),
|
||||
tuple(273.15 + np.array(temperatures)),
|
||||
)
|
||||
for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items()
|
||||
}
|
||||
|
||||
|
||||
# Same Toronto temperatures on a finer temperature mesh (every 6 minutes).
|
||||
TorontoTemperatures = {
|
||||
month: TorontoTemperatures_hourly[month].refine(refine_factor=10)
|
||||
for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items()
|
||||
}
|
||||
|
||||
|
||||
# references values for infection_probability and expected new cases
|
||||
# in the following tests, were obtained from the feature/mc branch
|
||||
|
||||
|
|
@ -22,7 +47,7 @@ def shared_office_mc():
|
|||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=50, humidity=0.5),
|
||||
ventilation = models.MultipleVentilation(
|
||||
ventilation=models.MultipleVentilation(
|
||||
ventilations=(
|
||||
models.SlidingWindow(
|
||||
active=models.PeriodicInterval(period=120, duration=120),
|
||||
|
|
@ -36,7 +61,7 @@ def shared_office_mc():
|
|||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
presence=mc.SpecificInterval(present_times = ((0, 3.5), (4.5, 9))),
|
||||
presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))),
|
||||
virus=virus_distributions['SARS_CoV_2_DELTA'],
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=activity_distributions['Seated'],
|
||||
|
|
@ -51,7 +76,7 @@ def shared_office_mc():
|
|||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=3,
|
||||
presence=mc.SpecificInterval(present_times = ((0, 3.5), (4.5, 9))),
|
||||
presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))),
|
||||
activity=activity_distributions['Seated'],
|
||||
mask=models.Mask.types['No mask'],
|
||||
host_immunity=0.,
|
||||
|
|
@ -66,12 +91,12 @@ def classroom_mc():
|
|||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=160, humidity=0.3),
|
||||
ventilation = models.MultipleVentilation(
|
||||
ventilation=models.MultipleVentilation(
|
||||
ventilations=(
|
||||
models.SlidingWindow(
|
||||
active=models.PeriodicInterval(period=120, duration=120),
|
||||
inside_temp=models.PiecewiseConstant((0., 24.), (293,)),
|
||||
outside_temp=data.TorontoTemperatures['Dec'],
|
||||
outside_temp=TorontoTemperatures['Dec'],
|
||||
window_height=1.6,
|
||||
opening_length=0.2,
|
||||
),
|
||||
|
|
@ -295,13 +320,13 @@ def waiting_room_mc():
|
|||
@pytest.mark.parametrize(
|
||||
"mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER",
|
||||
[
|
||||
["shared_office_mc", 6.03, 0.18, 27.22, 809],
|
||||
["classroom_mc", 9.5, 1.85, 79.98, 5624],
|
||||
["ski_cabin_mc", 16.0, 0.5, 40.25, 7966],
|
||||
["skagit_chorale_mc",65.7, 40.0, 241.28, 190422],
|
||||
["bus_ride_mc", 12.0, 8.0, 63.79, 5419],
|
||||
["gym_mc", 0.45, 0.13, 0.4852, 1145],
|
||||
["waiting_room_mc", 1.59, 0.22, 7.23, 737],
|
||||
["shared_office_mc", 6.03, 0.18, 3.198, 809],
|
||||
["classroom_mc", 9.5, 1.85, 9.478, 5624],
|
||||
["ski_cabin_mc", 16.0, 0.5, 17.315, 7966],
|
||||
["skagit_chorale_mc",65.7, 40.0, 102.213, 190422],
|
||||
["bus_ride_mc", 12.0, 8.0, 7.65, 5419],
|
||||
["gym_mc", 0.45, 0.13, 0.208, 1145],
|
||||
["waiting_room_mc", 1.59, 0.22, 0.821, 737],
|
||||
]
|
||||
)
|
||||
def test_report_models(mc_model, expected_pi, expected_new_cases,
|
||||
|
|
@ -312,7 +337,7 @@ def test_report_models(mc_model, expected_pi, expected_new_cases,
|
|||
expected_pi, rtol=TOLERANCE)
|
||||
npt.assert_allclose(exposure_model.expected_new_cases().mean(),
|
||||
expected_new_cases, rtol=TOLERANCE)
|
||||
npt.assert_allclose(exposure_model.exposure().mean(),
|
||||
npt.assert_allclose(exposure_model.deposited_exposure().mean(),
|
||||
expected_dose, rtol=TOLERANCE)
|
||||
npt.assert_allclose(
|
||||
exposure_model.concentration_model.infected.emission_rate_when_present().mean(),
|
||||
|
|
@ -322,10 +347,10 @@ def test_report_models(mc_model, expected_pi, expected_new_cases,
|
|||
@pytest.mark.parametrize(
|
||||
"mask_type, month, expected_pi, expected_dose, expected_ER",
|
||||
[
|
||||
["No mask", "Jul", 9.52, 84.54, 809],
|
||||
["Type I", "Jul", 1.7, 15.64, 149],
|
||||
["FFP2", "Jul", 0.51, 15.64, 149],
|
||||
["Type I", "Feb", 0.57, 4.59, 162],
|
||||
["No mask", "Jul", 9.52, 9.920, 809],
|
||||
["Type I", "Jul", 1.7, 0.913, 149],
|
||||
["FFP2", "Jul", 0.51, 0.239, 149],
|
||||
["Type I", "Feb", 0.57, 0.272, 162],
|
||||
],
|
||||
)
|
||||
def test_small_shared_office_Geneva(mask_type, month, expected_pi,
|
||||
|
|
@ -372,7 +397,7 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi,
|
|||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
npt.assert_allclose(exposure_model.infection_probability().mean(),
|
||||
expected_pi, rtol=TOLERANCE)
|
||||
npt.assert_allclose(exposure_model.exposure().mean(),
|
||||
npt.assert_allclose(exposure_model.deposited_exposure().mean(),
|
||||
expected_dose, rtol=TOLERANCE)
|
||||
npt.assert_allclose(
|
||||
exposure_model.concentration_model.infected.emission_rate_when_present().mean(),
|
||||
|
|
|
|||
Loading…
Reference in a new issue