Merge branch 'feature/MonteCarlo' into 'master'
Merge of MonteCarlo branch into master See merge request cara/cara!187
This commit is contained in:
Philip James Elson
commit
a00cfb1075
26 changed files with 910 additions and 231 deletions
|
|
@ -8,7 +8,9 @@ import numpy as np
|
|||
|
||||
from cara import models
|
||||
from cara import data
|
||||
import cara.monte_carlo as mc
|
||||
from .. import calculator
|
||||
from cara.monte_carlo.data import activity_distributions, virus_distributions
|
||||
|
||||
|
||||
LOG = logging.getLogger(__name__)
|
||||
|
|
@ -20,7 +22,7 @@ minutes_since_midnight = typing.NewType('minutes_since_midnight', int)
|
|||
# Used to declare when an attribute of a class must have a value provided, and
|
||||
# there should be no default value used.
|
||||
_NO_DEFAULT = object()
|
||||
|
||||
_SAMPLE_SIZE = 50000
|
||||
|
||||
@dataclass
|
||||
class FormData:
|
||||
|
|
@ -263,11 +265,15 @@ class FormData:
|
|||
ventilation = models.HVACMechanical(
|
||||
active=always_on, q_air_mech=self.air_supply)
|
||||
|
||||
# this is a minimal, always present source of ventilation, due
|
||||
# to the air infiltration from the outside.
|
||||
# See CERN-OPEN-2021-004, p. 12.
|
||||
infiltration_ventilation = models.AirChange(active=always_on, air_exch=0.25)
|
||||
if self.hepa_option:
|
||||
hepa = models.HEPAFilter(active=always_on, q_air_mech=self.hepa_amount)
|
||||
return models.MultipleVentilation((ventilation, hepa))
|
||||
return models.MultipleVentilation((ventilation, hepa, infiltration_ventilation))
|
||||
else:
|
||||
return ventilation
|
||||
return models.MultipleVentilation((ventilation, infiltration_ventilation))
|
||||
|
||||
def mask(self) -> models.Mask:
|
||||
# Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as
|
||||
|
|
@ -275,9 +281,9 @@ class FormData:
|
|||
mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask']
|
||||
return mask
|
||||
|
||||
def infected_population(self) -> models.InfectedPopulation:
|
||||
def infected_population(self) -> mc.InfectedPopulation:
|
||||
# Initializes the virus
|
||||
virus = models.Virus.types[self.virus_type]
|
||||
virus = virus_distributions[self.virus_type]
|
||||
|
||||
scenario_activity_and_expiration = {
|
||||
'office': (
|
||||
|
|
@ -315,12 +321,12 @@ class FormData:
|
|||
}
|
||||
|
||||
[activity_defn, expiration_defn] = scenario_activity_and_expiration[self.activity_type]
|
||||
activity = models.Activity.types[activity_defn]
|
||||
activity = activity_distributions[activity_defn]
|
||||
expiration = build_expiration(expiration_defn)
|
||||
|
||||
infected_occupants = self.infected_people
|
||||
|
||||
infected = models.InfectedPopulation(
|
||||
infected = mc.InfectedPopulation(
|
||||
number=infected_occupants,
|
||||
virus=virus,
|
||||
presence=self.infected_present_interval(),
|
||||
|
|
@ -330,7 +336,7 @@ class FormData:
|
|||
)
|
||||
return infected
|
||||
|
||||
def exposed_population(self) -> models.Population:
|
||||
def exposed_population(self) -> mc.Population:
|
||||
scenario_activity = {
|
||||
'office': 'Seated',
|
||||
'controlroom-day': 'Seated',
|
||||
|
|
@ -345,14 +351,14 @@ class FormData:
|
|||
}
|
||||
|
||||
activity_defn = scenario_activity[self.activity_type]
|
||||
activity = models.Activity.types[activity_defn]
|
||||
activity = activity_distributions[activity_defn]
|
||||
|
||||
infected_occupants = self.infected_people
|
||||
# The number of exposed occupants is the total number of occupants
|
||||
# minus the number of infected occupants.
|
||||
exposed_occupants = self.total_people - infected_occupants
|
||||
|
||||
exposed = models.Population(
|
||||
exposed = mc.Population(
|
||||
number=exposed_occupants,
|
||||
presence=self.exposed_present_interval(),
|
||||
activity=activity,
|
||||
|
|
@ -560,14 +566,14 @@ def model_from_form(form: FormData) -> models.ExposureModel:
|
|||
room = models.Room(volume=volume, humidity=humidity)
|
||||
|
||||
# Initializes and returns a model with the attributes defined above
|
||||
return models.ExposureModel(
|
||||
concentration_model=models.ConcentrationModel(
|
||||
return mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=room,
|
||||
ventilation=form.ventilation(),
|
||||
infected=form.infected_population(),
|
||||
),
|
||||
exposed=form.exposed_population()
|
||||
)
|
||||
).build_model(size=_SAMPLE_SIZE)
|
||||
|
||||
|
||||
def baseline_raw_form_data():
|
||||
|
|
|
|||
|
|
@ -38,12 +38,13 @@ def calculate_report_data(model: models.ExposureModel):
|
|||
|
||||
t_start, t_end = model_start_end(model)
|
||||
times = list(np.linspace(t_start, t_end, resolution))
|
||||
concentrations = [model.concentration_model.concentration(time) for time in times]
|
||||
concentrations = [np.mean(model.concentration_model.concentration(time))
|
||||
for time in times]
|
||||
highest_const = max(concentrations)
|
||||
prob = model.infection_probability()
|
||||
er = model.concentration_model.infected.emission_rate_when_present()
|
||||
prob = np.mean(model.infection_probability())
|
||||
er = np.mean(model.concentration_model.infected.emission_rate_when_present())
|
||||
exposed_occupants = model.exposed.number
|
||||
expected_new_cases = model.expected_new_cases()
|
||||
expected_new_cases = np.mean(model.expected_new_cases())
|
||||
|
||||
repeated_events = []
|
||||
for n in [1, 2, 3, 4, 5]:
|
||||
|
|
@ -52,8 +53,8 @@ def calculate_report_data(model: models.ExposureModel):
|
|||
repeated_events.append(
|
||||
RepeatEvents(
|
||||
repeats=n,
|
||||
probability_of_infection=repeat_model.infection_probability(),
|
||||
expected_new_cases=repeat_model.expected_new_cases(),
|
||||
probability_of_infection=np.mean(repeat_model.infection_probability()),
|
||||
expected_new_cases=np.mean(repeat_model.expected_new_cases()),
|
||||
)
|
||||
)
|
||||
|
||||
|
|
@ -127,13 +128,13 @@ def plot(times, concentrations, model: models.ExposureModel):
|
|||
fig = plt.figure()
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times]
|
||||
ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Concentration')
|
||||
ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean concentration')
|
||||
ax.spines['right'].set_visible(False)
|
||||
ax.spines['top'].set_visible(False)
|
||||
|
||||
ax.set_xlabel('Time of day')
|
||||
ax.set_ylabel('Concentration ($q/m^3$)')
|
||||
ax.set_title('Concentration of infectious quanta')
|
||||
ax.set_ylabel('Mean concentration ($q/m^3$)')
|
||||
ax.set_title('Mean concentration of infectious quanta')
|
||||
ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M"))
|
||||
|
||||
# Plot presence of exposed person
|
||||
|
|
@ -243,7 +244,8 @@ def comparison_plot(scenarios: typing.Dict[str, models.ExposureModel]):
|
|||
t_start, t_end = model_start_end(model)
|
||||
times = np.linspace(t_start, t_end, resolution)
|
||||
datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times]
|
||||
concentrations = [model.concentration_model.concentration(time) for time in times]
|
||||
concentrations = [np.mean(model.concentration_model.concentration(time))
|
||||
for time in times]
|
||||
|
||||
if name in dash_styled_scenarios:
|
||||
ax.plot(datetimes, concentrations, label=name, linestyle='--')
|
||||
|
|
@ -257,8 +259,8 @@ def comparison_plot(scenarios: typing.Dict[str, models.ExposureModel]):
|
|||
ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M"))
|
||||
|
||||
ax.set_xlabel('Time of day')
|
||||
ax.set_ylabel('Concentration ($q/m^3$)')
|
||||
ax.set_title('Concentration of infectious quanta')
|
||||
ax.set_ylabel('Mean concentration ($q/m^3$)')
|
||||
ax.set_title('Mean concentration of infectious quanta')
|
||||
|
||||
return fig
|
||||
|
||||
|
|
@ -267,8 +269,8 @@ def comparison_report(scenarios: typing.Dict[str, models.ExposureModel]):
|
|||
statistics = {}
|
||||
for name, model in scenarios.items():
|
||||
statistics[name] = {
|
||||
'probability_of_infection': model.infection_probability(),
|
||||
'expected_new_cases': model.expected_new_cases(),
|
||||
'probability_of_infection': np.mean(model.infection_probability()),
|
||||
'expected_new_cases': np.mean(model.expected_new_cases()),
|
||||
}
|
||||
return {
|
||||
'plot': img2base64(_figure2bytes(comparison_plot(scenarios))),
|
||||
|
|
|
|||
|
|
@ -209,35 +209,12 @@
|
|||
<p class="result_title">Results:</p>
|
||||
<p class="data_text">
|
||||
{% block report_summary %}
|
||||
In this scenario, the estimated probability of one exposed occupant getting infected P(i) is {{ prob_inf | non_zero_percentage }} and the expected number of new cases is {{ expected_new_cases | float_format }}.
|
||||
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><a href="#section1">[*]</a> and the <b>expected number of new cases is {{ expected_new_cases | float_format }}</b>.
|
||||
<p id="section1">[*] The results are based on the parameters and assumptions published in the CERN Open Report <a href="https://cds.cern.ch/record/2756083"> CERN-OPEN-2021-004</a></p>
|
||||
{% endblock report_summary %}
|
||||
</p>
|
||||
<p class="data_title">Exposure graph:</p>
|
||||
<img id="scenario_concentration_plot" src="{{ scenario_plot_src }}">
|
||||
|
||||
<p class="data_title">Repeated events:</p>
|
||||
<p class="data_text">
|
||||
The P(i) and expected number of new cases if repeating this scenario event - provided the infected person emits the same amount of viruses each day and the exposed person is subject to the same daily exposure time:
|
||||
|
||||
<table class="table table-striped w-auto">
|
||||
<thead class="thead-light">
|
||||
<tr>
|
||||
<th># of repeated events</th>
|
||||
<th>P(i)</th>
|
||||
<th>Expected new cases</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
{% for repeat_event in repeated_events %}
|
||||
<tr>
|
||||
<td>{{ repeat_event.repeats }}</td>
|
||||
<td>{{ repeat_event.probability_of_infection | non_zero_percentage }}</td>
|
||||
<td style="text-align:right">{{ repeat_event.expected_new_cases | float_format }}</td>
|
||||
</tr>
|
||||
{% endfor %}
|
||||
</tbody>
|
||||
</table>
|
||||
</p>
|
||||
<img id="scenario_concentration_plot" src="{{ scenario_plot_src }}">
|
||||
|
||||
<p class="data_title">Alternative scenarios:</p>
|
||||
<p class="data_text">
|
||||
|
|
@ -248,7 +225,7 @@
|
|||
<thead class="thead-light">
|
||||
<tr>
|
||||
<th>Scenario</th>
|
||||
<th>P(i)</th>
|
||||
<th>P(I)</th>
|
||||
<th>Expected new cases</th>
|
||||
</tr>
|
||||
</thead>
|
||||
|
|
@ -269,7 +246,7 @@
|
|||
<p class="data_text"> <strong> Notes for alternative scenarios: </strong><br>
|
||||
<ol>
|
||||
<li>This graph shows the concentration of infectious quanta in the air. The filtration of Type I and FFP2 masks, if worn, applies not only to the emission rate but also to the individual exposure (i.e. inhalation).
|
||||
For this reason, scenarios with different types of mask will show the same concentration on the graph but have different Pi values.</li>
|
||||
For this reason, scenarios with different types of mask will show the same concentration on the graph but have different absorbed doses and infection probabilities.</li>
|
||||
<li>If you have selected more sophisticated options, such as HEPA filtration or FFP2 masks, this will be indicated in the plot as the "base scenario", representing the inputs inserted in the form.<br>
|
||||
The other alternative scenarios shown for comparison will not include either HEPA filtration or FFP2 masks.</li>
|
||||
</ol>
|
||||
|
|
|
|||
|
|
@ -8,7 +8,7 @@
|
|||
<p>This is a guide to help you use the calculator app.
|
||||
If you are using the expert version of the tool, you should look at the expert notes.</p>
|
||||
<p>For more information on the Airborne Transmission of SARS-CoV-2, feel free to check out the HSE Seminar: <a href="https://cds.cern.ch/record/2743403">https://cds.cern.ch/record/2743403</a></p>
|
||||
<p>The methodology, mathematical equations and parameters of the model are described here: <a href="https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf">https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf</a></p>
|
||||
<p>The methodology, mathematical equations and parameters of the model are described here in the CERN Report: <a href="https://cds.cern.ch/record/2756083"> CERN-OPEN-2021-004</a></p>
|
||||
<h2>Disclaimer</h2>
|
||||
|
||||
<p>
|
||||
|
|
@ -172,26 +172,25 @@ Please check what are the applicable rules, before deciding which assumptions ar
|
|||
Please confirm what are the applicable rules, before deciding which assumptions are used for the simulation</p>
|
||||
<p>For the time being only the Type 1 surgical and FFP2 masks can be selected.</p>
|
||||
<h2>Generate Report</h2>
|
||||
<p>When you have entered all the necessary information, please click on the Generate Report button to execute the model.</p>
|
||||
<p>When you have entered all the necessary information, please click on the Generate Report button to execute the model. With the implementation of Monte Carlo simulations, the browser might take a few secounds to react.</p>
|
||||
<h1>Report</h1>
|
||||
<p>The report will open in your web browser.
|
||||
It contains a summary of all the input data, which will allow the simulation to be repeated if required in the future as we improve the model.</p>
|
||||
<h2>Results</h2>
|
||||
<p>This part of the report shows the <code>P(i)</code> or probability of one exposed person getting infected.
|
||||
<p>This part of the report shows the <code>P(I)</code> or probability of one exposed person getting infected.
|
||||
It is estimated based on the emission rate of virus into the simulated volume, and the amount which is inhaled by exposed individuals.
|
||||
This probability is valid for the simulation duration - i.e. if you have simulated one day and plan to work 5 days in these conditions and the infected person emits the same amount of virus each day, the cumulative probability of infection is <code>(1-(1-P(i))^5)</code>.
|
||||
This probability is valid for the simulation duration - i.e. the start and end time.
|
||||
If you are using the natural ventilation option, the simulation is only valid for the selected month, because the following or preceding month will have a different average temperature profile.
|
||||
The <code>expected number of new cases</code> for the simulation is calculated based on the probability of infection, multiplied by the number of exposed occupants.</p>
|
||||
<h3>Exposure graph</h3>
|
||||
<p>The graph shows the variation in the concentration of infectious quanta (one quanta is the amount of inhaled virus that can cause infection in 63% of the exposed occupants) within the simulated volume.
|
||||
<p>The graph shows the variation in the concentration of infectious viruses within the simulated volume.
|
||||
It is determined by:</p>
|
||||
<ul>
|
||||
<li>The presence of the infected person, who emits airborne viruses in the volume.</li>
|
||||
<li>The emission rate is related to the type of activity of the infected person (sitting, light exercise), their level of vocalisation (breathing, whispering or talking).</li>
|
||||
<li>The emission rate is related to the type of activity of the infected person (sitting, light exercise), their level of vocalisation (breathing, talking or shouting).</li>
|
||||
<li>The accumulation of infectious quanta in the volume, which is driven, among other factors, by ventilation (if applicable).<ul>
|
||||
<li>In a mechanical ventilation scenario, the removal rate is constant, based on fresh airflow supply in and out of the simulated space.</li>
|
||||
<li>Under natural ventilation conditions, the effectiveness of ventilation relies upon the hourly temperature difference between the inside and outside air temperature.</li>
|
||||
<li>A HEPA filter removes infectious quanta from the air at a constant rate and is modelled in the same way as mechanical ventilation, however air passed through a HEPA filter is recycled (i.e. it is not fresh air).</li>
|
||||
<li>A HEPA filter removes infectious virus from the air at a constant rate and is modelled in the same way as mechanical ventilation, however air passed through a HEPA filter is recycled (i.e. it is not fresh air).</li>
|
||||
</ul>
|
||||
</li>
|
||||
</ul>
|
||||
|
|
@ -204,7 +203,7 @@ This allows for:</p>
|
|||
</ul>
|
||||
<h1>Conclusion</h1>
|
||||
<p>This tool provides informative comparisons for COVID-19 (long-range) airborne risk only - see Disclaimer
|
||||
If you have any comments on your experience with the app, or feedback for potential improvements, please share them with the development team at cara-dev@cern.ch.</p>
|
||||
If you have any comments on your experience with the app, or feedback for potential improvements, please share them with the development team <a href="mailto:cara-dev@cern.ch">Send email</a>.</p>
|
||||
|
||||
|
||||
{% endblock contents %}
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@
|
|||
|
||||
{% block report_preamble %}
|
||||
<p><strong>Applicable rules: <br>
|
||||
Please ensure that this scenario conforms to current <a href="https://hse.cern/covid-19-information"> CERN HSE rules</a> (minimum ventilation requirements, mask wearing and the maximum number of people permitted in a space).</strong> <br>
|
||||
Please ensure that this scenario conforms to current COVID-related <a href="https://hse.cern/covid-19-information"> Health & Safety requirements</a>, under the applicable COVID Scale and measures in force at the time of the CARA assessment.</strong> <br>
|
||||
The results of this simulation are colour coded according to the risk values authorized at CERN (approved in December 2020):
|
||||
<ul><li>Events with a <span class="green_bkg">P(i) less than 5%</span> may go ahead without further mitigation measures.</li>
|
||||
<li>Events with a <span class="yellow_bkg">P(i) between 5% and 15%</span> shall be subject to ALARA principles (see footnote) to minimise the risk before proceeding.</li>
|
||||
|
|
|
|||
|
|
@ -17,7 +17,7 @@ CARA stands for COVID Airborne Risk Assessment and was developed in the spring o
|
|||
</ul>
|
||||
|
||||
The mathematical and physical model simulate the long-range airborne spread of SARS-CoV-2 virus in a finite volume, assuming a homogenous mixture, and estimates the risk of COVID-19 infection therein. The results DO NOT include short-range airborne exposure (where the physical distance plays a factor) nor the other known modes of SARS-CoV-2 transmission. Hence, the output from this model is only valid when the other recommended public health & safety instructions are observed, such as adequate physical distancing, good hand hygiene and other barrier measures.<br>
|
||||
<p>The methodology, mathematical equations and parameters of the model are described here: <a href="https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf">https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf</a>.</p>
|
||||
<p>The methodology, mathematical equations and parameters of the model are described here in the CERN Report: <a href="https://cds.cern.ch/record/2756083"> CERN-OPEN-2021-004</a>.</p>
|
||||
|
||||
The model used is based on scientific publications relating to airborne transmission of infectious diseases, virology, epidemiology and aerosol science. It can be used to compare the effectiveness of different airborne-related risk mitigation measures.
|
||||
|
||||
|
|
|
|||
|
|
@ -38,7 +38,7 @@ GenevaTemperatures_hourly = {
|
|||
}
|
||||
# same temperatures on a finer temperature mesh
|
||||
GenevaTemperatures = {
|
||||
month: GenevaTemperatures_hourly[month].refine(refine_factor=10)
|
||||
month: GenevaTemperatures_hourly[month].refine(refine_factor=4)
|
||||
for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items()
|
||||
}
|
||||
|
||||
|
|
|
|||
136
cara/models.py
136
cara/models.py
|
|
@ -37,6 +37,7 @@ import typing
|
|||
|
||||
import numpy as np
|
||||
from scipy.interpolate import interp1d
|
||||
import scipy.integrate
|
||||
|
||||
if not typing.TYPE_CHECKING:
|
||||
from memoization import cached
|
||||
|
|
@ -473,38 +474,47 @@ Virus.types = {
|
|||
|
||||
@dataclass(frozen=True)
|
||||
class Mask:
|
||||
#: Filtration efficiency. (In %/100)
|
||||
η_exhale: _VectorisedFloat
|
||||
|
||||
#: Leakage through side of masks.
|
||||
η_leaks: _VectorisedFloat
|
||||
|
||||
#: Filtration efficiency of masks when inhaling.
|
||||
η_inhale: _VectorisedFloat
|
||||
|
||||
#: Global factor applied to filtration efficiency of masks when exhaling.
|
||||
factor_exhale: float = 1.
|
||||
|
||||
#: Pre-populated examples of Masks.
|
||||
types: typing.ClassVar[typing.Dict[str, "Mask"]]
|
||||
|
||||
def exhale_efficiency(self, diameter: float) -> _VectorisedFloat:
|
||||
# Overall efficiency with the effect of the leaks for aerosol emission
|
||||
# Gammaitoni et al (1997). Diameter is in cm.
|
||||
if diameter < 3e-4:
|
||||
"""
|
||||
Overall exhale efficiency, including the effect of the leaks.
|
||||
See CERN-OPEN-2021-004 (doi: 10.17181/CERN.1GDQ.5Y75), and Ref.
|
||||
therein (Asadi 2020).
|
||||
Obtained from measurements of filtration efficiency and of
|
||||
the leakage through the sides.
|
||||
Diameter is in microns.
|
||||
"""
|
||||
if diameter < 0.5:
|
||||
eta_out = 0.
|
||||
elif diameter < 0.94614:
|
||||
eta_out = 0.5893 * diameter + 0.1546
|
||||
elif diameter < 3.:
|
||||
eta_out = 0.0509 * diameter + 0.664
|
||||
else:
|
||||
eta_out = self.η_exhale * (1 - self.η_leaks)
|
||||
return eta_out
|
||||
eta_out = 0.8167
|
||||
return eta_out*self.factor_exhale
|
||||
|
||||
def inhale_efficiency(self) -> _VectorisedFloat:
|
||||
"""
|
||||
Overall inhale efficiency, including the effect of the leaks.
|
||||
"""
|
||||
return self.η_inhale
|
||||
|
||||
|
||||
Mask.types = {
|
||||
'No mask': Mask(0, 0, 0),
|
||||
'No mask': Mask(0, 0),
|
||||
'Type I': Mask(
|
||||
η_exhale=0.95,
|
||||
η_leaks=0.15, # (Huang 2007)
|
||||
η_inhale=0.3, # (Browen 2010)
|
||||
η_inhale=0.5, # (CERN-OPEN-2021-004)
|
||||
),
|
||||
'FFP2': Mask(
|
||||
η_exhale=0.95, # (same outward effect as type 1 - Asadi 2020)
|
||||
η_leaks=0.15, # (same outward effect as type 1 - Asadi 2020)
|
||||
η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149)
|
||||
),
|
||||
}
|
||||
|
|
@ -516,35 +526,59 @@ class _ExpirationBase:
|
|||
Represents the expiration of aerosols by a person.
|
||||
Subclasses of _ExpirationBase represent different models.
|
||||
"""
|
||||
#: Pre-populated examples of Masks.
|
||||
#: Pre-populated examples of Expirations.
|
||||
types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
|
||||
|
||||
def aerosols(self, mask: Mask):
|
||||
# total volume of aerosols expired per volume of air (mL/cm^3).
|
||||
"""
|
||||
total volume of aerosols expired per volume of air (mL/cm^3).
|
||||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Expiration(_ExpirationBase):
|
||||
"""
|
||||
Simple model based on four different sizes of particles emitted,
|
||||
with different ejection factors. See Fig. 4 in L. Morawska et al,
|
||||
Size distribution and sites of origin of droplets expelled from the
|
||||
human respiratory tract during expiratory activities,
|
||||
Aerosol Science 40 (2009) pp. 256 - 269.
|
||||
BLO model for the expiration (G. Johnson et al., Modality of human
|
||||
expired aerosol size distributions, Journal of Aerosol Science,
|
||||
vol. 42, no. 12, pp. 839 – 851, 2011,
|
||||
https://doi.org/10.1016/j.jaerosci.2011.07.009).
|
||||
Here all diameters (d) are in microns.
|
||||
"""
|
||||
ejection_factor: typing.Tuple[float, ...]
|
||||
particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4) # In cm.
|
||||
#: factors assigned to resp. the B, L and O modes. They are
|
||||
# charateristics of the kind of expiratory activity (e.g. breathing,
|
||||
# speaking, singing, or shouting).
|
||||
BLO_factors: typing.Tuple[float, float, float]
|
||||
|
||||
@cached()
|
||||
def aerosols(self, mask: Mask):
|
||||
def volume(diameter):
|
||||
return (4 * np.pi * (diameter/2)**3) / 3
|
||||
total = 0
|
||||
for diameter, factor in zip(self.particle_sizes, self.ejection_factor):
|
||||
contribution = (volume(diameter) * factor *
|
||||
(1 - mask.exhale_efficiency(diameter)))
|
||||
total += contribution
|
||||
return total
|
||||
""" Result is in mL.cm^-3 """
|
||||
def volume(d):
|
||||
return (np.pi * d**3) / 6.
|
||||
|
||||
def _Bmode(d: float) -> float:
|
||||
# B-mode (see ref. above).
|
||||
return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) *
|
||||
np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2)))
|
||||
|
||||
def _Lmode(d: float) -> float:
|
||||
# L-mode (see ref. above).
|
||||
return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.506818)) *
|
||||
np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2)))
|
||||
|
||||
def _Omode(d: float) -> float:
|
||||
# O-mode (see ref. above).
|
||||
return ( (1 / d) * (0.0010008 / (np.sqrt(2 * np.pi) * 0.585005)) *
|
||||
np.exp(-1 * (np.log(d) - 4.97673) ** 2 / (2 * 0.585005 ** 2)))
|
||||
|
||||
def integrand(d: float) -> float:
|
||||
return (self.BLO_factors[0] * _Bmode(d) +
|
||||
self.BLO_factors[1] * _Lmode(d) +
|
||||
self.BLO_factors[2] * _Omode(d)
|
||||
) * volume(d) * (1 - mask.exhale_efficiency(d))
|
||||
|
||||
# final result converted from microns^3/cm3 to mL/cm^3
|
||||
return scipy.integrate.quad(integrand, 0.1, 30.)[0]*1e-12
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
|
|
@ -572,17 +606,20 @@ class MultipleExpiration(_ExpirationBase):
|
|||
|
||||
|
||||
_ExpirationBase.types = {
|
||||
'Breathing': Expiration((0.084, 0.009, 0.003, 0.002)),
|
||||
'Whispering': Expiration((0.11, 0.014, 0.004, 0.002)),
|
||||
'Talking': Expiration((0.236, 0.068, 0.007, 0.011)),
|
||||
'Unmodulated Vocalization': Expiration((0.751, 0.139, 0.0139, 0.059)),
|
||||
'Superspreading event': Expiration((np.inf, np.inf, np.inf, np.inf)),
|
||||
'Breathing': Expiration((1., 0., 0.)),
|
||||
'Talking': Expiration((1., 1., 1.)),
|
||||
'Shouting': Expiration((1., 5., 5.)),
|
||||
'Singing': Expiration((1., 5., 5.)),
|
||||
'Superspreading event': Expiration((np.inf, 0., 0.)),
|
||||
}
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Activity:
|
||||
#: Inhalation rate in m^3/h
|
||||
inhalation_rate: _VectorisedFloat
|
||||
|
||||
#: Exhalation rate in m^3/h
|
||||
exhalation_rate: _VectorisedFloat
|
||||
|
||||
#: Pre-populated examples of activities.
|
||||
|
|
@ -637,6 +674,8 @@ class InfectedPopulation(Population):
|
|||
|
||||
"""
|
||||
# Emission Rate (infectious quantum / 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
|
||||
aerosols = self.expiration.aerosols(self.mask)
|
||||
|
||||
ER = (self.virus.viral_load_in_sputum *
|
||||
|
|
@ -672,7 +711,6 @@ class InfectedPopulation(Population):
|
|||
|
||||
return self.emission_rate_when_present()
|
||||
|
||||
@cached()
|
||||
def emission_rate(self, time) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate of the entire population.
|
||||
|
|
@ -692,8 +730,8 @@ class ConcentrationModel:
|
|||
return self.infected.virus
|
||||
|
||||
def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat:
|
||||
# Particle deposition on the floor
|
||||
vg = 1 * 10 ** -4
|
||||
# Particle deposition on the floor (value from CERN-OPEN-2021-04)
|
||||
vg = 1.88e-4
|
||||
# Height of the emission source to the floor - i.e. mouth/nose (m)
|
||||
h = 1.5
|
||||
# Deposition rate (h^-1)
|
||||
|
|
@ -702,7 +740,6 @@ class ConcentrationModel:
|
|||
return k + self.virus.decay_constant(self.room.humidity
|
||||
) + self.ventilation.air_exchange(self.room, time)
|
||||
|
||||
@cached()
|
||||
def _concentration_limit(self, time: float) -> _VectorisedFloat:
|
||||
"""
|
||||
Provides a constant that represents the theoretical asymptotic
|
||||
|
|
@ -714,7 +751,6 @@ class ConcentrationModel:
|
|||
|
||||
return (self.infected.emission_rate(time)) / (IVRR * V)
|
||||
|
||||
@cached()
|
||||
def state_change_times(self):
|
||||
"""
|
||||
All time dependent entities on this model must provide information about
|
||||
|
|
@ -759,6 +795,9 @@ class ConcentrationModel:
|
|||
return (self.last_state_change(stop) <= start)
|
||||
|
||||
@cached()
|
||||
def _concentration_at_state_change(self, time: float) -> _VectorisedFloat:
|
||||
return self.concentration(time)
|
||||
|
||||
def concentration(self, time: float) -> _VectorisedFloat:
|
||||
"""
|
||||
Virus quanta concentration, as a function of time.
|
||||
|
|
@ -777,7 +816,7 @@ class ConcentrationModel:
|
|||
concentration_limit = self._concentration_limit(next_state_change_time)
|
||||
|
||||
t_last_state_change = self.last_state_change(time)
|
||||
concentration_at_last_state_change = self.concentration(t_last_state_change)
|
||||
concentration_at_last_state_change = self._concentration_at_state_change(t_last_state_change)
|
||||
|
||||
delta_time = time - t_last_state_change
|
||||
fac = np.exp(-IVRR * delta_time)
|
||||
|
|
@ -821,6 +860,9 @@ class ExposureModel:
|
|||
#: The number of times the exposure event is repeated (default 1).
|
||||
repeats: int = 1
|
||||
|
||||
#: The fraction of viruses actually deposited in the respiratory tract
|
||||
fraction_deposited: _VectorisedFloat = 0.6
|
||||
|
||||
def quanta_exposure(self) -> _VectorisedFloat:
|
||||
"""The number of virus quanta per meter^3."""
|
||||
exposure = 0.0
|
||||
|
|
@ -835,8 +877,8 @@ class ExposureModel:
|
|||
|
||||
inf_aero = (
|
||||
self.exposed.activity.inhalation_rate *
|
||||
(1 - self.exposed.mask.η_inhale) *
|
||||
exposure
|
||||
(1 - self.exposed.mask.inhale_efficiency()) *
|
||||
exposure * self.fraction_deposited
|
||||
)
|
||||
|
||||
# Probability of infection.
|
||||
|
|
|
|||
56
cara/monte_carlo/data.py
Normal file
56
cara/monte_carlo/data.py
Normal file
|
|
@ -0,0 +1,56 @@
|
|||
import numpy as np
|
||||
|
||||
import cara.monte_carlo as mc
|
||||
from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel
|
||||
|
||||
|
||||
# From CERN-OPEN-2021-04 and refererences therein
|
||||
activity_distributions = {
|
||||
'Seated': mc.Activity(LogNormal(-0.6872121723362303, 0.10498338229297108),
|
||||
LogNormal(-0.6872121723362303, 0.10498338229297108)),
|
||||
|
||||
'Standing': mc.Activity(LogNormal(-0.5742377578494785, 0.09373162411398223),
|
||||
LogNormal(-0.5742377578494785, 0.09373162411398223)),
|
||||
|
||||
'Light activity': mc.Activity(LogNormal(0.21380242785625422,0.09435378091059601),
|
||||
LogNormal(0.21380242785625422,0.09435378091059601)),
|
||||
|
||||
'Moderate activity': mc.Activity(LogNormal(0.551771330362601, 0.1894616357138137),
|
||||
LogNormal(0.551771330362601, 0.1894616357138137)),
|
||||
|
||||
'Heavy exercise': mc.Activity(LogNormal(1.1644665696723049, 0.21744554768657565),
|
||||
LogNormal(1.1644665696723049, 0.21744554768657565)),
|
||||
}
|
||||
|
||||
|
||||
# From CERN-OPEN-2021-04 and refererences therein
|
||||
symptomatic_vl_frequencies = LogCustomKernel(
|
||||
np.array((2.46032, 2.67431, 2.85434, 3.06155, 3.25856, 3.47256, 3.66957, 3.85979, 4.09927, 4.27081,
|
||||
4.47631, 4.66653, 4.87204, 5.10302, 5.27456, 5.46478, 5.6533, 5.88428, 6.07281, 6.30549,
|
||||
6.48552, 6.64856, 6.85407, 7.10373, 7.30075, 7.47229, 7.66081, 7.85782, 8.05653, 8.27053,
|
||||
8.48453, 8.65607, 8.90573, 9.06878, 9.27429, 9.473, 9.66152, 9.87552)),
|
||||
np.array((0.001206885, 0.007851618, 0.008078144, 0.01502491, 0.013258014, 0.018528495, 0.020053765,
|
||||
0.021896167, 0.022047184, 0.018604005, 0.01547796, 0.018075445, 0.021503523, 0.022349217,
|
||||
0.025097721, 0.032875078, 0.030594727, 0.032573045, 0.034717482, 0.034792991,
|
||||
0.033267721, 0.042887485, 0.036846816, 0.03876473, 0.045016819, 0.040063473, 0.04883754,
|
||||
0.043944602, 0.048142864, 0.041588741, 0.048762031, 0.027921732, 0.033871788,
|
||||
0.022122693, 0.016927718, 0.008833228, 0.00478598, 0.002807662)),
|
||||
kernel_bandwidth=0.1
|
||||
)
|
||||
|
||||
|
||||
# From CERN-OPEN-2021-04 and refererences therein
|
||||
virus_distributions = {
|
||||
'SARS_CoV_2': mc.SARSCoV2(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
quantum_infectious_dose=100,
|
||||
),
|
||||
'SARS_CoV_2_B117': mc.SARSCoV2(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
quantum_infectious_dose=60,
|
||||
),
|
||||
'SARS_CoV_2_P1': mc.SARSCoV2(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
quantum_infectious_dose=100/2.25,
|
||||
),
|
||||
}
|
||||
|
|
@ -1,6 +1,7 @@
|
|||
import typing
|
||||
|
||||
import numpy as np
|
||||
from sklearn.neighbors import KernelDensity # type: ignore
|
||||
|
||||
import cara.models
|
||||
|
||||
|
|
@ -16,6 +17,9 @@ class SampleableDistribution:
|
|||
|
||||
|
||||
class Normal(SampleableDistribution):
|
||||
"""
|
||||
Defines a normal (i.e. Gaussian) distribution
|
||||
"""
|
||||
def __init__(self, mean: float, standard_deviation: float):
|
||||
self.mean = mean
|
||||
self.standard_deviation = standard_deviation
|
||||
|
|
@ -24,6 +28,96 @@ class Normal(SampleableDistribution):
|
|||
return np.random.normal(self.mean, self.standard_deviation, size=size)
|
||||
|
||||
|
||||
class LogNormal(SampleableDistribution):
|
||||
"""
|
||||
Defines a lognormal distribution (i.e. Gaussian distribution vs. the
|
||||
natural logarithm of the random variable)
|
||||
"""
|
||||
|
||||
def __init__(self, mean_gaussian: float, standard_deviation_gaussian: float):
|
||||
# these are resp. the mean and std. deviation of the underlying
|
||||
# Gaussian distribution
|
||||
self.mean_gaussian = mean_gaussian
|
||||
self.standard_deviation_gaussian = standard_deviation_gaussian
|
||||
|
||||
def generate_samples(self, size: int) -> float_array_size_n:
|
||||
return np.random.lognormal(self.mean_gaussian,
|
||||
self.standard_deviation_gaussian,
|
||||
size=size)
|
||||
|
||||
|
||||
class Custom(SampleableDistribution):
|
||||
"""
|
||||
Defines a distribution which follows a custom curve vs. the random
|
||||
variable. Uses a simple algorithm. This is appropriate for a smooth
|
||||
distribution function (one should know its maximum).
|
||||
"""
|
||||
def __init__(self, bounds: typing.Tuple[float, float],
|
||||
function: typing.Callable, max_function: float):
|
||||
self.bounds = bounds
|
||||
self.function = function
|
||||
self.max_function = max_function
|
||||
|
||||
def generate_samples(self, size: int) -> float_array_size_n:
|
||||
fvalue = np.random.uniform(0,self.max_function,size)
|
||||
x = np.random.uniform(*self.bounds,size)
|
||||
invalid = np.where(fvalue>self.function(x))[0]
|
||||
while len(invalid)>0:
|
||||
fvalue[invalid] = np.random.uniform(0,self.max_function,len(invalid))
|
||||
x[invalid] = np.random.uniform(*self.bounds,len(invalid))
|
||||
invalid = np.where(fvalue>self.function(x))[0]
|
||||
|
||||
return x
|
||||
|
||||
|
||||
class CustomKernel(SampleableDistribution):
|
||||
"""
|
||||
Defines a distribution which follows a custom curve vs. the
|
||||
random variable. Uses a Gaussian kernel density fit. This is more
|
||||
appropriate for a noisy distribution function.
|
||||
"""
|
||||
def __init__(self, variable: float_array_size_n,
|
||||
frequencies: float_array_size_n,
|
||||
kernel_bandwidth: float):
|
||||
# these are resp. the random variable, the distribution
|
||||
# frequencies at these values, and the bandwidth of the Gaussian
|
||||
# kernel
|
||||
self.variable = variable
|
||||
self.frequencies = frequencies
|
||||
self.kernel_bandwidth = kernel_bandwidth
|
||||
|
||||
def generate_samples(self, size: int) -> float_array_size_n:
|
||||
kde_model = KernelDensity(kernel='gaussian',
|
||||
bandwidth=self.kernel_bandwidth)
|
||||
kde_model.fit(self.variable.reshape(-1, 1),
|
||||
sample_weight=self.frequencies)
|
||||
return kde_model.sample(n_samples=size)[:, 0]
|
||||
|
||||
|
||||
class LogCustomKernel(SampleableDistribution):
|
||||
"""
|
||||
Defines a distribution which follows a custom curve vs. the log
|
||||
(in base 10) of the random variable. Uses a Gaussian kernel density
|
||||
fit. This is more appropriate for a noisy distribution function.
|
||||
"""
|
||||
def __init__(self, log_variable: float_array_size_n,
|
||||
frequencies: float_array_size_n,
|
||||
kernel_bandwidth: float):
|
||||
# these are resp. the log of the random variable, the distribution
|
||||
# frequencies at these values, and the bandwidth of the Gaussian
|
||||
# kernel
|
||||
self.log_variable = log_variable
|
||||
self.frequencies = frequencies
|
||||
self.kernel_bandwidth = kernel_bandwidth
|
||||
|
||||
def generate_samples(self, size: int) -> float_array_size_n:
|
||||
kde_model = KernelDensity(kernel='gaussian',
|
||||
bandwidth=self.kernel_bandwidth)
|
||||
kde_model.fit(self.log_variable.reshape(-1, 1),
|
||||
sample_weight=self.frequencies)
|
||||
return 10 ** kde_model.sample(n_samples=size)[:, 0]
|
||||
|
||||
|
||||
_VectorisedFloatOrSampleable = typing.Union[
|
||||
SampleableDistribution, cara.models._VectorisedFloat,
|
||||
]
|
||||
|
|
|
|||
|
|
@ -50,7 +50,7 @@ def test_ventilation_slidingwindow(baseline_form: model_generator.FormData):
|
|||
baseline_form.opening_distance = 0.6
|
||||
|
||||
ts = np.linspace(8, 16, 100)
|
||||
np.testing.assert_allclose([window.air_exchange(room, t) for t in ts],
|
||||
np.testing.assert_allclose([window.air_exchange(room, t)+0.25 for t in ts],
|
||||
[baseline_form.ventilation().air_exchange(room, t) for t in ts])
|
||||
|
||||
|
||||
|
|
@ -73,7 +73,7 @@ def test_ventilation_hingedwindow(baseline_form: model_generator.FormData):
|
|||
baseline_form.opening_distance = 0.6
|
||||
|
||||
ts = np.linspace(8, 16, 100)
|
||||
np.testing.assert_allclose([window.air_exchange(room, t) for t in ts],
|
||||
np.testing.assert_allclose([window.air_exchange(room, t)+0.25 for t in ts],
|
||||
[baseline_form.ventilation().air_exchange(room, t) for t in ts])
|
||||
|
||||
|
||||
|
|
@ -88,7 +88,7 @@ def test_ventilation_mechanical(baseline_form: model_generator.FormData):
|
|||
baseline_form.air_supply = 500.
|
||||
|
||||
ts = np.linspace(8, 16, 100)
|
||||
np.testing.assert_allclose([mech.air_exchange(room, t) for t in ts],
|
||||
np.testing.assert_allclose([mech.air_exchange(room, t)+0.25 for t in ts],
|
||||
[baseline_form.ventilation().air_exchange(room, t) for t in ts])
|
||||
|
||||
|
||||
|
|
@ -103,7 +103,7 @@ def test_ventilation_airchanges(baseline_form: model_generator.FormData):
|
|||
baseline_form.air_changes = 3.
|
||||
|
||||
ts = np.linspace(8, 16, 100)
|
||||
np.testing.assert_allclose([airchange.air_exchange(room, t) for t in ts],
|
||||
np.testing.assert_allclose([airchange.air_exchange(room, t)+0.25 for t in ts],
|
||||
[baseline_form.ventilation().air_exchange(room, t) for t in ts])
|
||||
|
||||
|
||||
|
|
@ -131,7 +131,7 @@ def test_ventilation_window_hepa(baseline_form: model_generator.FormData):
|
|||
baseline_form.hepa_option = True
|
||||
|
||||
ts = np.linspace(9, 17, 100)
|
||||
np.testing.assert_allclose([ventilation.air_exchange(room, t) for t in ts],
|
||||
np.testing.assert_allclose([ventilation.air_exchange(room, t)+0.25 for t in ts],
|
||||
[baseline_form.ventilation().air_exchange(room, t) for t in ts])
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -11,7 +11,7 @@ def test_generate_report(baseline_form):
|
|||
# generate a report for it. Because this is what happens in the cara
|
||||
# calculator, we confirm that the generation happens within a reasonable
|
||||
# time threshold.
|
||||
time_limit: float = 5.0 # seconds
|
||||
time_limit: float = 20.0 # seconds
|
||||
|
||||
start = time.perf_counter()
|
||||
|
||||
|
|
|
|||
|
|
@ -6,6 +6,7 @@ import tornado.testing
|
|||
import cara.apps.calculator
|
||||
from cara.apps.calculator.report_generator import generate_qr_code
|
||||
|
||||
_TIMEOUT = 20.
|
||||
|
||||
@pytest.fixture
|
||||
def app():
|
||||
|
|
@ -45,11 +46,12 @@ class TestBasicApp(tornado.testing.AsyncHTTPTestCase):
|
|||
def get_app(self):
|
||||
return cara.apps.calculator.make_app()
|
||||
|
||||
@tornado.testing.gen_test(timeout=_TIMEOUT)
|
||||
def test_report(self):
|
||||
response = self.fetch('/calculator/baseline-model/result')
|
||||
response = yield self.http_client.fetch(self.get_url('/calculator/baseline-model/result'))
|
||||
self.assertEqual(response.code, 200)
|
||||
assert 'CERN HSE rules' not in response.body.decode()
|
||||
assert 'the expected number of new cases is' in response.body.decode()
|
||||
assert 'CERN HSE' not in response.body.decode()
|
||||
assert 'expected number of new cases is' in response.body.decode()
|
||||
|
||||
|
||||
class TestCernApp(tornado.testing.AsyncHTTPTestCase):
|
||||
|
|
@ -57,11 +59,12 @@ class TestCernApp(tornado.testing.AsyncHTTPTestCase):
|
|||
cern_theme = Path(cara.apps.calculator.__file__).parent / 'themes' / 'cern'
|
||||
return cara.apps.calculator.make_app(theme_dir=cern_theme)
|
||||
|
||||
@tornado.testing.gen_test(timeout=_TIMEOUT)
|
||||
def test_report(self):
|
||||
response = self.fetch('/calculator/baseline-model/result')
|
||||
response = yield self.http_client.fetch(self.get_url('/calculator/baseline-model/result'))
|
||||
self.assertEqual(response.code, 200)
|
||||
assert 'CERN HSE rules' in response.body.decode()
|
||||
assert 'the expected number of new cases is' in response.body.decode()
|
||||
assert 'CERN HSE' in response.body.decode()
|
||||
assert 'expected number of new cases is' in response.body.decode()
|
||||
|
||||
|
||||
async def test_qrcode_urls(http_server_client, baseline_form):
|
||||
|
|
|
|||
|
|
@ -7,11 +7,9 @@ import pytest
|
|||
def baseline_model():
|
||||
model = models.ConcentrationModel(
|
||||
room=models.Room(volume=75),
|
||||
ventilation=models.SlidingWindow(
|
||||
active=models.PeriodicInterval(period=120, duration=120),
|
||||
inside_temp=models.PiecewiseConstant((0,24),(293,)),
|
||||
outside_temp=models.PiecewiseConstant((0,24),(283,)),
|
||||
window_height=1.6, opening_length=0.6,
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=30.,
|
||||
),
|
||||
infected=models.InfectedPopulation(
|
||||
number=1,
|
||||
|
|
@ -19,7 +17,7 @@ def baseline_model():
|
|||
presence=models.SpecificInterval(((0, 4), (5, 8))),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Light activity'],
|
||||
expiration=models.Expiration.types['Unmodulated Vocalization'],
|
||||
expiration=models.Expiration.types['Superspreading event'],
|
||||
),
|
||||
)
|
||||
return model
|
||||
|
|
@ -30,9 +28,10 @@ def baseline_exposure_model(baseline_model):
|
|||
return models.ExposureModel(
|
||||
baseline_model,
|
||||
exposed=models.Population(
|
||||
number=10,
|
||||
number=1000,
|
||||
presence=baseline_model.infected.presence,
|
||||
activity=baseline_model.infected.activity,
|
||||
mask=baseline_model.infected.mask,
|
||||
)
|
||||
),
|
||||
fraction_deposited = 1.,
|
||||
)
|
||||
|
|
|
|||
|
|
@ -1,6 +1,7 @@
|
|||
import re
|
||||
|
||||
import numpy as np
|
||||
import numpy.testing as npt
|
||||
import pytest
|
||||
|
||||
from cara import models
|
||||
|
|
@ -13,8 +14,6 @@ from cara import models
|
|||
{'air_change': np.array([100, 120])},
|
||||
{'viral_load_in_sputum': np.array([5e8, 1e9])},
|
||||
{'quantum_infectious_dose': np.array([50, 20])},
|
||||
{'η_exhale': np.array([0.92, 0.95])},
|
||||
{'η_leaks': np.array([0.15, 0.20])},
|
||||
]
|
||||
)
|
||||
def test_concentration_model_vectorisation(override_params):
|
||||
|
|
@ -24,8 +23,6 @@ def test_concentration_model_vectorisation(override_params):
|
|||
'air_change': 100,
|
||||
'viral_load_in_sputum': 1e9,
|
||||
'quantum_infectious_dose': 50,
|
||||
'η_exhale': 0.95,
|
||||
'η_leaks': 0.15,
|
||||
}
|
||||
defaults.update(override_params)
|
||||
|
||||
|
|
@ -37,8 +34,7 @@ def test_concentration_model_vectorisation(override_params):
|
|||
number=1,
|
||||
presence=always,
|
||||
mask=models.Mask(
|
||||
η_exhale=defaults['η_exhale'],
|
||||
η_leaks=defaults['η_leaks'],
|
||||
factor_exhale=0.95,
|
||||
η_inhale=0.3,
|
||||
),
|
||||
activity=models.Activity(
|
||||
|
|
@ -49,9 +45,7 @@ def test_concentration_model_vectorisation(override_params):
|
|||
viral_load_in_sputum=defaults['viral_load_in_sputum'],
|
||||
quantum_infectious_dose=defaults['quantum_infectious_dose'],
|
||||
),
|
||||
expiration=models.Expiration(
|
||||
ejection_factor=(0.084, 0.009, 0.003, 0.002),
|
||||
),
|
||||
expiration=models.Expiration((1., 0., 0.)),
|
||||
)
|
||||
)
|
||||
concentrations = c_model.concentration(10)
|
||||
|
|
@ -128,4 +122,4 @@ def test_integrated_concentration(simple_conc_model):
|
|||
c2 = simple_conc_model.integrated_concentration(0, 1)
|
||||
c3 = simple_conc_model.integrated_concentration(1, 2)
|
||||
assert c1 != 0
|
||||
assert c1 == c2 + c3
|
||||
npt.assert_almost_equal(c1, c2 + c3, decimal=15)
|
||||
|
|
|
|||
|
|
@ -43,7 +43,7 @@ populations = [
|
|||
),
|
||||
# A population with some array component for η_inhale.
|
||||
models.Population(
|
||||
10, halftime, models.Mask(0.95, 0.15, np.array([0.3, 0.35])),
|
||||
10, halftime, models.Mask(np.array([0.3, 0.35])),
|
||||
models.Activity.types['Standing'],
|
||||
),
|
||||
# A population with some array component for inhalation_rate.
|
||||
|
|
@ -55,22 +55,25 @@ populations = [
|
|||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"population, cm, expected_exposure, expected_probability",[
|
||||
[populations[1], KnownConcentrations(lambda t: 1.2),
|
||||
"population, cm, f_dep, expected_exposure, expected_probability",[
|
||||
[populations[1], KnownConcentrations(lambda t: 1.2), 1.,
|
||||
np.array([14.4, 14.4]), np.array([99.6803184113, 99.5181053773])],
|
||||
|
||||
[populations[2], KnownConcentrations(lambda t: 1.2),
|
||||
np.array([14.4, 14.4]), np.array([99.4146994564, 99.6803184113])],
|
||||
[populations[2], KnownConcentrations(lambda t: 1.2), 1.,
|
||||
np.array([14.4, 14.4]), np.array([97.4574432074, 98.3493482895])],
|
||||
|
||||
[populations[0], KnownConcentrations(lambda t: np.array([1.2, 2.4])),
|
||||
np.array([14.4, 28.8]), np.array([99.6803184113, 99.9989780368])],
|
||||
[populations[0], KnownConcentrations(lambda t: np.array([1.2, 2.4])), 1.,
|
||||
np.array([14.4, 28.8]), np.array([98.3493482895, 99.9727534893])],
|
||||
|
||||
[populations[1], KnownConcentrations(lambda t: np.array([1.2, 2.4])),
|
||||
[populations[1], KnownConcentrations(lambda t: np.array([1.2, 2.4])), 1.,
|
||||
np.array([14.4, 28.8]), np.array([99.6803184113, 99.9976777757])],
|
||||
|
||||
[populations[0], KnownConcentrations(lambda t: 2.4), np.array([0.5, 1.]),
|
||||
28.8, np.array([98.3493482895, 99.9727534893])],
|
||||
])
|
||||
def test_exposure_model_ndarray(population, cm,
|
||||
def test_exposure_model_ndarray(population, cm, f_dep,
|
||||
expected_exposure, expected_probability):
|
||||
model = ExposureModel(cm, population)
|
||||
model = ExposureModel(cm, population, fraction_deposited = f_dep)
|
||||
np.testing.assert_almost_equal(
|
||||
model.quanta_exposure(), expected_exposure
|
||||
)
|
||||
|
|
@ -123,22 +126,22 @@ def conc_model():
|
|||
models.InfectedPopulation(
|
||||
number=1,
|
||||
presence=interesting_times,
|
||||
mask=models.Mask.types['Type I'],
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Seated'],
|
||||
virus=models.Virus.types['SARS_CoV_2'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
expiration=models.Expiration.types['Superspreading event'],
|
||||
)
|
||||
)
|
||||
|
||||
# expected quanta were computed with a trapezoidal integration, using
|
||||
# a mesh of 10'000 pts per exposed presence interval.
|
||||
@pytest.mark.parametrize("exposed_time_interval, expected_quanta", [
|
||||
[(0, 1), 0.0055680845],
|
||||
[(1, 1.01), 6.4960491e-05],
|
||||
[(1.01, 1.02), 6.3187723e-05],
|
||||
[(12, 12.01), 1.9307359e-06],
|
||||
[(12, 24), 0.079347465],
|
||||
[(0, 24), 0.086122050],
|
||||
[(0, 1), 5.3334352],
|
||||
[(1, 1.01), 0.061759078],
|
||||
[(1.01, 1.02), 0.060016487],
|
||||
[(12, 12.01), 0.0019012647],
|
||||
[(12, 24), 75.513005],
|
||||
[(0, 24), 81.956988],
|
||||
]
|
||||
)
|
||||
def test_exposure_model_integral_accuracy(exposed_time_interval,
|
||||
|
|
@ -148,5 +151,5 @@ def test_exposure_model_integral_accuracy(exposed_time_interval,
|
|||
10, presence_interval, models.Mask.types['Type I'],
|
||||
models.Activity.types['Standing'],
|
||||
)
|
||||
model = ExposureModel(conc_model, population)
|
||||
model = ExposureModel(conc_model, population, fraction_deposited=1.)
|
||||
np.testing.assert_allclose(model.quanta_exposure(), expected_quanta)
|
||||
|
|
|
|||
33
cara/tests/models/test_mask.py
Normal file
33
cara/tests/models/test_mask.py
Normal file
|
|
@ -0,0 +1,33 @@
|
|||
import numpy as np
|
||||
import numpy.testing as npt
|
||||
import pytest
|
||||
|
||||
from cara import models
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"η_inhale, expected_inhale_efficiency",
|
||||
[
|
||||
[0.5, 0.5],
|
||||
[np.array([0.3, 0.5]), np.array([0.3, 0.5])],
|
||||
],
|
||||
)
|
||||
def test_mask_inhale(η_inhale, expected_inhale_efficiency):
|
||||
mask = models.Mask(η_inhale=η_inhale)
|
||||
npt.assert_equal(mask.inhale_efficiency(),
|
||||
expected_inhale_efficiency)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"diameter, factor_exhale, expected_exhale_efficiency",
|
||||
[
|
||||
[0.3, 1., 0.],
|
||||
[0.7, 0.3, 0.56711*0.3],
|
||||
[1., 1., 0.7149],
|
||||
[4., 0.5, 0.8167*0.5],
|
||||
[5., 0., 0.],
|
||||
],
|
||||
)
|
||||
def test_mask_exhale(diameter, factor_exhale, expected_exhale_efficiency):
|
||||
mask = models.Mask(η_inhale=0.3, factor_exhale=factor_exhale)
|
||||
npt.assert_almost_equal(mask.exhale_efficiency(diameter),
|
||||
expected_exhale_efficiency)
|
||||
|
|
@ -9,7 +9,7 @@ from cara import models
|
|||
|
||||
def test_multiple_wrong_weight_size():
|
||||
weights = (1., 2., 3.)
|
||||
e_base = models.Expiration((0.084, 0.009, 0.003, 0.002))
|
||||
e_base = models.Expiration((1., 0., 0.))
|
||||
with pytest.raises(
|
||||
ValueError,
|
||||
match=re.escape("expirations and weigths should contain the"
|
||||
|
|
@ -19,11 +19,25 @@ def test_multiple_wrong_weight_size():
|
|||
|
||||
|
||||
def test_multiple():
|
||||
weights = (1., 2.)
|
||||
e1 = models.Expiration((0.03, 0.02, 0.01, 0.005))
|
||||
e2 = models.Expiration((0.05, 0.04, 0.03, 0.01))
|
||||
weights = (1., 1.)
|
||||
e1 = models.Expiration((1., 0., 0.))
|
||||
e2 = models.Expiration((4., 5., 5.))
|
||||
e_expected = models.Expiration((2.5, 2.5, 2.5))
|
||||
e = models.MultipleExpiration([e1, e2], weights)
|
||||
assert e.aerosols(models.Mask.types['No mask']) == (
|
||||
e1.aerosols(models.Mask.types['No mask'])/3. +
|
||||
2*e2.aerosols(models.Mask.types['No mask'])/3.
|
||||
)
|
||||
mask = models.Mask.types['Type I']
|
||||
npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask))
|
||||
|
||||
|
||||
# expected values obtained from another code
|
||||
@pytest.mark.parametrize(
|
||||
"BLO_weights, expected_aerosols",
|
||||
[
|
||||
[(1.,0.,0.), 1.38924e-12],
|
||||
[(1.,1.,1.), 1.07129e-10],
|
||||
[(1.,5.,5.), 5.30088e-10],
|
||||
],
|
||||
)
|
||||
def test_expiration_aerosols(BLO_weights, expected_aerosols):
|
||||
mask = models.Mask.types['No mask']
|
||||
e = models.Expiration(BLO_weights)
|
||||
npt.assert_allclose(e.aerosols(mask), expected_aerosols, rtol=1e-4)
|
||||
|
|
|
|||
|
|
@ -8,8 +8,6 @@ import cara.models
|
|||
"override_params", [
|
||||
{'viral_load_in_sputum': np.array([5e8, 1e9])},
|
||||
{'quantum_infectious_dose': np.array([50, 20])},
|
||||
{'η_exhale': np.array([0.92, 0.95])},
|
||||
{'η_leaks': np.array([0.15, 0.20])},
|
||||
{'exhalation_rate': np.array([0.75, 0.81])},
|
||||
]
|
||||
)
|
||||
|
|
@ -17,8 +15,6 @@ def test_infected_population_vectorisation(override_params):
|
|||
defaults = {
|
||||
'viral_load_in_sputum': 1e9,
|
||||
'quantum_infectious_dose': 50,
|
||||
'η_exhale': 0.95,
|
||||
'η_leaks': 0.15,
|
||||
'exhalation_rate': 0.75,
|
||||
}
|
||||
defaults.update(override_params)
|
||||
|
|
@ -28,8 +24,7 @@ def test_infected_population_vectorisation(override_params):
|
|||
number=1,
|
||||
presence=office_hours,
|
||||
mask=cara.models.Mask(
|
||||
η_exhale=defaults['η_exhale'],
|
||||
η_leaks=defaults['η_leaks'],
|
||||
factor_exhale=0.95,
|
||||
η_inhale=0.3,
|
||||
),
|
||||
activity=cara.models.Activity(
|
||||
|
|
@ -40,9 +35,7 @@ def test_infected_population_vectorisation(override_params):
|
|||
viral_load_in_sputum=defaults['viral_load_in_sputum'],
|
||||
quantum_infectious_dose=defaults['quantum_infectious_dose'],
|
||||
),
|
||||
expiration=cara.models.Expiration(
|
||||
ejection_factor=(0.084, 0.009, 0.003, 0.002),
|
||||
),
|
||||
expiration=cara.models.Expiration((1., 0., 0.)),
|
||||
)
|
||||
emission_rate = infected.emission_rate(10)
|
||||
assert isinstance(emission_rate, np.ndarray)
|
||||
|
|
|
|||
|
|
@ -6,21 +6,12 @@ import cara.models as models
|
|||
import cara.data as data
|
||||
|
||||
|
||||
def test_no_mask_aerosols(baseline_model):
|
||||
exp = models.Expiration.types['Unmodulated Vocalization']
|
||||
npt.assert_allclose(
|
||||
exp.aerosols(models.Mask.types['No mask']),
|
||||
6.077541e-12,
|
||||
rtol=1e-5,
|
||||
)
|
||||
|
||||
|
||||
def test_no_mask_emission_rate(baseline_model):
|
||||
rate = 151.938514
|
||||
def test_no_mask_superspeading_emission_rate(baseline_model):
|
||||
expected_rate = 970.
|
||||
npt.assert_allclose(
|
||||
[baseline_model.infected.emission_rate(t) for t in [0, 1, 4, 4.5, 5, 8, 9]],
|
||||
[0, rate, rate, 0, 0, rate, 0],
|
||||
rtol=1e-5
|
||||
[0, expected_rate, expected_rate, 0, 0, expected_rate, 0],
|
||||
rtol=1e-12
|
||||
)
|
||||
|
||||
|
||||
|
|
@ -48,23 +39,24 @@ def baseline_periodic_hepa():
|
|||
|
||||
|
||||
def test_concentrations(baseline_model):
|
||||
# expected concentrations were computed analytically
|
||||
ts = [0, 4, 5, 7, 10]
|
||||
concentrations = [baseline_model.concentration(t) for t in ts]
|
||||
npt.assert_allclose(
|
||||
concentrations,
|
||||
[0.000000e+00, 2.619538e-01, 1.146999e-04, 2.619537e-01, 5.022289e-08],
|
||||
rtol=1e-5
|
||||
[0.000000e+00, 0.41611256, 1.3205628e-14, 0.41611256, 4.1909001e-28],
|
||||
rtol=1e-6
|
||||
)
|
||||
|
||||
|
||||
def test_smooth_concentrations(baseline_model):
|
||||
# We don't care about the actual concentrations in this test, but rather
|
||||
# that the curve itself is smooth.
|
||||
dx = 0.1
|
||||
dy_limit = dx * 2 # Anything more than this is a bit steep.
|
||||
dx = 0.002
|
||||
dy_limit = 0.2 # Anything more than this (in relative) is a bit steep.
|
||||
ts = np.arange(0, 10, dx)
|
||||
concentrations = [baseline_model.concentration(t) for t in ts]
|
||||
assert np.abs(np.diff(concentrations)).max() < dy_limit
|
||||
assert np.abs(np.diff(concentrations)).max()/np.mean(concentrations) < dy_limit
|
||||
|
||||
|
||||
def build_model(interval_duration):
|
||||
|
|
@ -80,7 +72,7 @@ def build_model(interval_duration):
|
|||
presence=models.SpecificInterval(((0, 4), (5, 8))),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Light activity'],
|
||||
expiration=models.Expiration.types['Unmodulated Vocalization'],
|
||||
expiration=models.Expiration.types['Superspreading event'],
|
||||
),
|
||||
)
|
||||
return model
|
||||
|
|
@ -98,8 +90,8 @@ def test_concentrations_startup(baseline_model):
|
|||
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.
|
||||
p = baseline_exposure_model.infection_probability()
|
||||
npt.assert_allclose(p, 89.001748)
|
||||
r0 = baseline_exposure_model.reproduction_number()
|
||||
npt.assert_allclose(r0, 972.880852)
|
||||
|
||||
|
||||
def test_periodic_window(baseline_periodic_window, baseline_room):
|
||||
|
|
@ -183,15 +175,6 @@ def test_multiple_ventilation_HEPA_HVAC_AirChange(volume, expected_value):
|
|||
expected_value,rtol=1e-5)
|
||||
|
||||
|
||||
def test_expiration_aerosols():
|
||||
mask = models.Mask.types['Type I']
|
||||
exp1 = models.Expiration((0.751, 0.139, 0.0139, 0.059),
|
||||
particle_sizes = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4))
|
||||
exp2 = models.Expiration((0.059, 0.0139, 0.751, 0.139),
|
||||
particle_sizes = (5.5e-4, 3.5e-4, 0.8e-4, 1.8e-4))
|
||||
npt.assert_allclose(exp1.aerosols(mask), exp2.aerosols(mask), rtol=1e-5)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"time, expected_value",
|
||||
[
|
||||
|
|
@ -239,7 +222,7 @@ def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),),
|
|||
presence=models.SpecificInterval(intervals_presence_infected),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Light activity'],
|
||||
expiration=models.Expiration.types['Unmodulated Vocalization'],
|
||||
expiration=models.Expiration.types['Superspreading event'],
|
||||
),
|
||||
)
|
||||
return model
|
||||
|
|
@ -260,7 +243,7 @@ def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)):
|
|||
presence=models.SpecificInterval(((0, 4), (5, 7.5))),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Light activity'],
|
||||
expiration=models.Expiration.types['Unmodulated Vocalization'],
|
||||
expiration=models.Expiration.types['Superspreading event'],
|
||||
),
|
||||
)
|
||||
return model
|
||||
|
|
@ -287,7 +270,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5
|
|||
presence=models.SpecificInterval(((0, 4), (5, 7.5))),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Light activity'],
|
||||
expiration=models.Expiration.types['Unmodulated Vocalization'],
|
||||
expiration=models.Expiration.types['Superspreading event'],
|
||||
),
|
||||
)
|
||||
return model
|
||||
|
|
@ -366,20 +349,21 @@ def build_exposure_model(concentration_model):
|
|||
presence=infected.presence,
|
||||
activity=infected.activity,
|
||||
mask=infected.mask,
|
||||
)
|
||||
),
|
||||
fraction_deposited = 1.,
|
||||
)
|
||||
|
||||
|
||||
# expected r0 were computed with a trapezoidal integration, using
|
||||
# expected quanta were computed with a trapezoidal integration, using
|
||||
# a mesh of 100'000 pts per exposed presence interval.
|
||||
@pytest.mark.parametrize(
|
||||
"month, expected_r0",
|
||||
"month, expected_quanta",
|
||||
[
|
||||
['Jan', 86.258533],
|
||||
['Jun', 99.972103],
|
||||
['Jan', 9.930854],
|
||||
['Jun', 37.962708],
|
||||
],
|
||||
)
|
||||
def test_r0_hourly_dep(month,expected_r0):
|
||||
def test_quanta_hourly_dep(month,expected_quanta):
|
||||
m = build_exposure_model(
|
||||
build_hourly_dependent_model(
|
||||
month,
|
||||
|
|
@ -387,19 +371,20 @@ def test_r0_hourly_dep(month,expected_r0):
|
|||
intervals_presence_infected=((8, 12), (13, 17))
|
||||
)
|
||||
)
|
||||
p = m.infection_probability()
|
||||
npt.assert_allclose(p, expected_r0)
|
||||
quanta = m.quanta_exposure()
|
||||
npt.assert_allclose(quanta, expected_quanta)
|
||||
|
||||
# expected r0 were computed with a trapezoidal integration, using
|
||||
# a mesh of 100'000 pts per exposed presence interval.
|
||||
# expected quanta 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_r0",
|
||||
"month, expected_quanta",
|
||||
[
|
||||
['Jan', 86.422736],
|
||||
['Jun', 99.982323],
|
||||
['Jan', 9.993842],
|
||||
['Jun', 40.151985],
|
||||
],
|
||||
)
|
||||
def test_r0_hourly_dep_refined(month,expected_r0):
|
||||
def test_quanta_hourly_dep_refined(month,expected_quanta):
|
||||
m = build_exposure_model(
|
||||
build_hourly_dependent_model(
|
||||
month,
|
||||
|
|
@ -408,5 +393,5 @@ def test_r0_hourly_dep_refined(month,expected_r0):
|
|||
temperatures=data.GenevaTemperatures,
|
||||
)
|
||||
)
|
||||
p = m.infection_probability()
|
||||
npt.assert_allclose(p, expected_r0)
|
||||
quanta = m.quanta_exposure()
|
||||
npt.assert_allclose(quanta, expected_quanta, rtol=0.02)
|
||||
|
|
|
|||
|
|
@ -53,7 +53,7 @@ def baseline_mc_model() -> cara.monte_carlo.ConcentrationModel:
|
|||
presence=cara.models.SpecificInterval(((0, 4), (5, 8))),
|
||||
mask=cara.models.Mask.types['No mask'],
|
||||
activity=cara.models.Activity.types['Light activity'],
|
||||
expiration=cara.models.Expiration.types['Unmodulated Vocalization'],
|
||||
expiration=cara.models.Expiration.types['Breathing'],
|
||||
),
|
||||
)
|
||||
return mc_model
|
||||
|
|
@ -76,7 +76,9 @@ def test_build_concentration_model(baseline_mc_model: cara.monte_carlo.Concentra
|
|||
model = baseline_mc_model.build_model(7)
|
||||
assert isinstance(model, cara.models.ConcentrationModel)
|
||||
assert isinstance(model.concentration(time=0), float)
|
||||
assert model.concentration(time=1).shape == (7, )
|
||||
conc = model.concentration(time=1)
|
||||
assert isinstance(conc, np.ndarray)
|
||||
assert conc.shape == (7, )
|
||||
|
||||
|
||||
def test_build_exposure_model(baseline_mc_exposure_model: cara.monte_carlo.ExposureModel):
|
||||
|
|
|
|||
316
cara/tests/test_monte_carlo_full_models.py
Normal file
316
cara/tests/test_monte_carlo_full_models.py
Normal file
|
|
@ -0,0 +1,316 @@
|
|||
import numpy as np
|
||||
import numpy.testing as npt
|
||||
import pytest
|
||||
|
||||
import cara.monte_carlo as mc
|
||||
from cara import models,data
|
||||
from cara.monte_carlo.data import activity_distributions, virus_distributions
|
||||
|
||||
# TODO: seed better the random number generators
|
||||
np.random.seed(2000)
|
||||
SAMPLE_SIZE = 50000
|
||||
TOLERANCE = 0.05
|
||||
|
||||
# references values for infection_probability and expected new cases
|
||||
# in the following tests, were obtained from the feature/mc branch
|
||||
|
||||
@pytest.fixture
|
||||
def shared_office_mc():
|
||||
"""
|
||||
Corresponds to the 1st line of Table 5 in CERN-OPEN-2021-04, but
|
||||
speaking 30% of the time (instead of 1/3)
|
||||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=50, humidity=0.3),
|
||||
ventilation=models.MultipleVentilation(
|
||||
(
|
||||
models.SlidingWindow(
|
||||
active=models.PeriodicInterval(period=120, duration=10),
|
||||
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
|
||||
outside_temp=models.PiecewiseConstant((0, 24), (283,)),
|
||||
window_height=1.6, opening_length=0.6,
|
||||
),
|
||||
models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
),
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=virus_distributions['SARS_CoV_2_B117'],
|
||||
presence=mc.SpecificInterval(((0, 2), (2.1, 4), (5, 7), (7.1, 9))),
|
||||
mask=models.Mask(η_inhale=0.3),
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.MultipleExpiration(
|
||||
expirations=(models.Expiration.types['Talking'],
|
||||
models.Expiration.types['Breathing']),
|
||||
weights=(0.3, 0.7)),
|
||||
),
|
||||
)
|
||||
return mc.ExposureModel(
|
||||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=3,
|
||||
presence=concentration_mc.infected.presence,
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=concentration_mc.infected.mask,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def classroom_mc():
|
||||
"""
|
||||
Corresponds to the 2nd line of Table 5 in CERN-OPEN-2021-04
|
||||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=160, humidity=0.3),
|
||||
ventilation=models.MultipleVentilation(
|
||||
(
|
||||
models.SlidingWindow(
|
||||
active=models.PeriodicInterval(period=120, duration=10),
|
||||
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
|
||||
outside_temp=models.PiecewiseConstant((0, 24), (283,)),
|
||||
window_height=1.6, opening_length=0.6,
|
||||
),
|
||||
models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
),
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=virus_distributions['SARS_CoV_2_B117'],
|
||||
presence=mc.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=activity_distributions['Light activity'],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
),
|
||||
)
|
||||
return mc.ExposureModel(
|
||||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=19,
|
||||
presence=concentration_mc.infected.presence,
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=concentration_mc.infected.mask,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def ski_cabin_mc():
|
||||
"""
|
||||
Corresponds to the 3rd line of Table 5 in CERN-OPEN-2021-04
|
||||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=10, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=0,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=virus_distributions['SARS_CoV_2_B117'],
|
||||
presence=mc.SpecificInterval(((0, 1/3),)),
|
||||
mask=models.Mask(η_inhale=0.3),
|
||||
activity=activity_distributions['Moderate activity'],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
),
|
||||
)
|
||||
return mc.ExposureModel(
|
||||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=3,
|
||||
presence=concentration_mc.infected.presence,
|
||||
activity=models.Activity.types['Moderate activity'],
|
||||
mask=concentration_mc.infected.mask,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def gym_mc():
|
||||
"""
|
||||
Corresponds to the 4th line of Table 5 in CERN-OPEN-2021-04,
|
||||
but there the expected number of cases is wrongly reported as 0.56
|
||||
while it should be 0.63.
|
||||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=300, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=6,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=2,
|
||||
virus=virus_distributions['SARS_CoV_2_B117'],
|
||||
presence=mc.SpecificInterval(((0, 1),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Heavy exercise'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
),
|
||||
)
|
||||
return mc.ExposureModel(
|
||||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=28,
|
||||
presence=concentration_mc.infected.presence,
|
||||
activity=models.Activity.types['Heavy exercise'],
|
||||
mask=concentration_mc.infected.mask,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def waiting_room_mc():
|
||||
"""
|
||||
Corresponds to the 5th line of Table 5 in CERN-OPEN-2021-04, but
|
||||
speaking 30% of the time (instead of 20%)
|
||||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=virus_distributions['SARS_CoV_2_B117'],
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.MultipleExpiration(
|
||||
expirations=(models.Expiration.types['Talking'],
|
||||
models.Expiration.types['Breathing']),
|
||||
weights=(0.3, 0.7)),
|
||||
),
|
||||
)
|
||||
return mc.ExposureModel(
|
||||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=concentration_mc.infected.presence,
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=concentration_mc.infected.mask,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def skagit_chorale_mc():
|
||||
"""
|
||||
Corresponds to the 6th line of Table 5 in CERN-OPEN-2021-04, but
|
||||
infection probability should be 29.8% instead of 32%, and number of
|
||||
new cases 17.9 instead of 21.
|
||||
"""
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=810, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=0.7,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=virus_distributions['SARS_CoV_2'],
|
||||
presence=mc.SpecificInterval(((0, 2.5),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Light activity'],
|
||||
expiration=models.Expiration((5., 5., 5.)),
|
||||
),
|
||||
)
|
||||
return mc.ExposureModel(
|
||||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=60,
|
||||
presence=concentration_mc.infected.presence,
|
||||
activity=models.Activity.types['Moderate activity'],
|
||||
mask=concentration_mc.infected.mask,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"mc_model, expected_pi, expected_new_cases, expected_dose, expected_qR",
|
||||
[
|
||||
["shared_office_mc", 10.7, 0.32, 0.954, 10.9],
|
||||
["classroom_mc", 36.1, 6.85, 13.0, 474.4],
|
||||
["ski_cabin_mc", 16.3, 0.49, 0.599, 123.4],
|
||||
["gym_mc", 2.25, 0.63, 0.01307, 16.4],
|
||||
["waiting_room_mc", 9.72, 1.36, 0.571, 58.9],
|
||||
["skagit_chorale_mc",29.9, 17.9, 1.90, 1414],
|
||||
]
|
||||
)
|
||||
def test_report_models(mc_model, expected_pi, expected_new_cases,
|
||||
expected_dose, expected_qR, request):
|
||||
mc_model = request.getfixturevalue(mc_model)
|
||||
exposure_model = mc_model.build_model(size=SAMPLE_SIZE)
|
||||
npt.assert_allclose(exposure_model.infection_probability().mean(),
|
||||
expected_pi, rtol=TOLERANCE)
|
||||
npt.assert_allclose(exposure_model.expected_new_cases().mean(),
|
||||
expected_new_cases, rtol=TOLERANCE)
|
||||
npt.assert_allclose(exposure_model.quanta_exposure().mean(),
|
||||
expected_dose, rtol=TOLERANCE)
|
||||
npt.assert_allclose(
|
||||
exposure_model.concentration_model.infected.emission_rate_when_present().mean(),
|
||||
expected_qR, rtol=TOLERANCE)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"mask_type, month, expected_pi, expected_dose, expected_qR",
|
||||
[
|
||||
["No mask", "Jul", 30.0, 6.764, 64.9],
|
||||
["Type I", "Jul", 10.2, 1.223, 11.7],
|
||||
["FFP2", "Jul", 4.0, 1.223, 11.7],
|
||||
["Type I", "Feb", 4.25, 0.357, 11.7],
|
||||
],
|
||||
)
|
||||
def test_small_shared_office_Geneva(mask_type, month, expected_pi,
|
||||
expected_dose, expected_qR):
|
||||
concentration_mc = mc.ConcentrationModel(
|
||||
room=models.Room(volume=33, humidity=0.5),
|
||||
ventilation=models.MultipleVentilation(
|
||||
(
|
||||
models.SlidingWindow(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
|
||||
outside_temp=data.GenevaTemperatures[month],
|
||||
window_height=1.5, opening_length=0.2,
|
||||
),
|
||||
models.AirChange(
|
||||
active=models.SpecificInterval(((0,24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
),
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=virus_distributions['SARS_CoV_2_B117'],
|
||||
presence=mc.SpecificInterval(((9, 10+2/3), (10+5/6, 12.5), (13.5, 15+2/3), (15+5/6, 18))),
|
||||
mask=models.Mask.types[mask_type],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.MultipleExpiration(
|
||||
expirations=(models.Expiration.types['Talking'],
|
||||
models.Expiration.types['Breathing']),
|
||||
weights=(0.33, 0.67)),
|
||||
),
|
||||
)
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=concentration_mc,
|
||||
exposed=mc.Population(
|
||||
number=1,
|
||||
presence=concentration_mc.infected.presence,
|
||||
activity=activity_distributions['Seated'],
|
||||
mask=concentration_mc.infected.mask,
|
||||
),
|
||||
)
|
||||
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.quanta_exposure().mean(),
|
||||
expected_dose, rtol=TOLERANCE)
|
||||
npt.assert_allclose(
|
||||
exposure_model.concentration_model.infected.emission_rate_when_present().mean(),
|
||||
expected_qR, rtol=TOLERANCE)
|
||||
47
cara/tests/test_predefined_distributions.py
Normal file
47
cara/tests/test_predefined_distributions.py
Normal file
|
|
@ -0,0 +1,47 @@
|
|||
import numpy as np
|
||||
import numpy.testing as npt
|
||||
import pytest
|
||||
|
||||
from cara.monte_carlo.data import activity_distributions, virus_distributions
|
||||
|
||||
# TODO: seed better the random number generators
|
||||
np.random.seed(2000)
|
||||
|
||||
|
||||
# mean & std deviations from CERN-OPEN-2021-04 (Table 4)
|
||||
# NOTE: a mistake was corrected for the std deviation of the
|
||||
# "Moderate exercise" case (0.37 in the report, but should be 0.34)
|
||||
@pytest.mark.parametrize(
|
||||
"distribution, mean, std",[
|
||||
['Seated', 0.51, 0.053],
|
||||
|
||||
['Standing', 0.57, 0.053],
|
||||
|
||||
['Light activity', 1.24, 0.12],
|
||||
|
||||
['Moderate activity', 1.77, 0.34],
|
||||
|
||||
['Heavy exercise', 3.28, 0.72,],
|
||||
]
|
||||
)
|
||||
def test_activity_distributions(distribution, mean, std):
|
||||
activity = activity_distributions[distribution].build_model(size=1000000)
|
||||
npt.assert_allclose(activity.inhalation_rate.mean(), mean, atol=0.01)
|
||||
npt.assert_allclose(activity.inhalation_rate.std(), std, atol=0.01)
|
||||
|
||||
|
||||
# mean & std deviations from CERN-OPEN-2021-04 (Table 4) - with a refined
|
||||
# precision on the values
|
||||
@pytest.mark.parametrize(
|
||||
"distribution, mean, std",[
|
||||
['SARS_CoV_2', 6.59, 1.74],
|
||||
|
||||
['SARS_CoV_2_B117', 6.59, 1.74],
|
||||
|
||||
['SARS_CoV_2_P1', 6.59, 1.74],
|
||||
]
|
||||
)
|
||||
def test_viral_load_logdistribution(distribution, mean, std):
|
||||
virus = virus_distributions[distribution].build_model(size=1000000)
|
||||
npt.assert_allclose(np.log10(virus.viral_load_in_sputum).mean(), mean, atol=0.01)
|
||||
npt.assert_allclose(np.log10(virus.viral_load_in_sputum).std(), std, atol=0.01)
|
||||
112
cara/tests/test_sampleable_distribution.py
Normal file
112
cara/tests/test_sampleable_distribution.py
Normal file
|
|
@ -0,0 +1,112 @@
|
|||
import numpy as np
|
||||
import numpy.testing as npt
|
||||
import pytest
|
||||
|
||||
from cara.monte_carlo import sampleable
|
||||
|
||||
# TODO: seed better the random number generators
|
||||
np.random.seed(2000)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"mean, std",[
|
||||
[1., 0.5],
|
||||
]
|
||||
)
|
||||
def test_normal(mean, std):
|
||||
# test that the sample has approximately the right mean,
|
||||
# std deviation and distribution function.
|
||||
sample_size = 2000000
|
||||
samples = sampleable.Normal(mean, std).generate_samples(sample_size)
|
||||
histogram, bins = np.histogram(samples,bins=100, density=True)
|
||||
selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip(
|
||||
(bins[1:]+bins[:-1])/2,histogram) if b>=0.25 and b<=1.75])
|
||||
exact_dist = 1/(np.sqrt(2*np.pi)*std) * np.exp(
|
||||
-((np.array(selected_bins)-mean)/std)**2/2)
|
||||
|
||||
assert len(samples) == sample_size
|
||||
npt.assert_allclose([samples.mean(), samples.std()], [mean, std], rtol=0.01)
|
||||
npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"mean_gaussian, std_gaussian",[
|
||||
[-0.6872121723362303, 0.10498338229297108],
|
||||
]
|
||||
)
|
||||
def test_lognormal(mean_gaussian, std_gaussian):
|
||||
# test that the sample has approximately the right mean,
|
||||
# std deviation and distribution function.
|
||||
sample_size = 2000000
|
||||
samples = sampleable.LogNormal(mean_gaussian, std_gaussian
|
||||
).generate_samples(sample_size)
|
||||
histogram, bins = np.histogram(samples,bins=50, density=True)
|
||||
selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip(
|
||||
(bins[1:]+bins[:-1])/2,histogram) if b>=0.4 and b<=0.6])
|
||||
selected_bins = np.array(selected_bins)
|
||||
exact_dist = ( 1/(selected_bins*np.sqrt(2*np.pi)*std_gaussian) *
|
||||
np.exp(-((np.log(selected_bins)-mean_gaussian)/std_gaussian)**2/2) )
|
||||
exact_mean = np.exp(mean_gaussian + std_gaussian**2/2)
|
||||
exact_std = np.sqrt( (np.exp(std_gaussian**2)-1) *
|
||||
np.exp(2*mean_gaussian + std_gaussian**2) )
|
||||
|
||||
assert len(samples) == sample_size
|
||||
npt.assert_allclose([samples.mean(), samples.std()],
|
||||
[exact_mean, exact_std], rtol=0.01)
|
||||
npt.assert_allclose(selected_histogram, exact_dist, rtol=0.03)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"use_kernel",
|
||||
[False, True],
|
||||
)
|
||||
def test_custom(use_kernel):
|
||||
# test that the sample has approximately the right distribution
|
||||
# function, with both Custom and CustomKernel method. The latter
|
||||
# is less accurate for smooth functions.
|
||||
# the distribution function is an inverted parabola, with maximum 0.15,
|
||||
# which is 0 at the bounds (0,10) (normalized)
|
||||
norm = 500/3.
|
||||
function = lambda x: (-(5 - x)**2 + 25)/norm
|
||||
max_function = 0.15
|
||||
sample_size = 2000000
|
||||
|
||||
if use_kernel:
|
||||
variable = np.linspace(0.1,9.9,100)
|
||||
frequencies = function(variable)
|
||||
samples = sampleable.CustomKernel(variable, frequencies,
|
||||
kernel_bandwidth=0.1
|
||||
).generate_samples(sample_size)
|
||||
else:
|
||||
samples = sampleable.Custom((0, 10), function, max_function
|
||||
).generate_samples(sample_size)
|
||||
|
||||
histogram, bins = np.histogram(samples, bins=100, density=True)
|
||||
selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip(
|
||||
(bins[1:]+bins[:-1])/2,histogram) if b>=1 and b<=9])
|
||||
correct_dist = function(np.array(selected_bins))
|
||||
assert len(samples) == sample_size
|
||||
npt.assert_allclose(selected_histogram, correct_dist, rtol=0.05)
|
||||
|
||||
|
||||
def test_logcustomkernel():
|
||||
# test that the sample has approximately the right distribution
|
||||
# function, for the LogCustomKernel.
|
||||
# the distribution function is an inverted parabola vs. the log of
|
||||
# the variable (normalized)
|
||||
norm = 500/3.
|
||||
function = lambda x: (-(5 - x)**2 + 25)/norm
|
||||
sample_size = 2000000
|
||||
|
||||
log_variable = np.linspace(0.1,9.9,100)
|
||||
frequencies = function(log_variable)
|
||||
samples = sampleable.LogCustomKernel(log_variable, frequencies,
|
||||
kernel_bandwidth=0.1
|
||||
).generate_samples(sample_size)
|
||||
|
||||
histogram, bins = np.histogram(np.log10(samples), bins=100, density=True)
|
||||
selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip(
|
||||
(bins[1:]+bins[:-1])/2,histogram) if b>=1 and b<=9])
|
||||
correct_dist = function(np.array(selected_bins))
|
||||
assert len(samples) == sample_size
|
||||
npt.assert_allclose(selected_histogram, correct_dist, rtol=0.05)
|
||||
|
|
@ -62,6 +62,7 @@ python-dateutil==2.8.1
|
|||
pyzmq==22.0.3
|
||||
qrcode==6.1
|
||||
scipy==1.5.4
|
||||
scikit_learn==0.23.1
|
||||
Send2Trash==1.5.0
|
||||
six==1.15.0
|
||||
sniffio==1.2.0
|
||||
|
|
|
|||
1
setup.py
1
setup.py
|
|
@ -29,6 +29,7 @@ REQUIREMENTS: dict = {
|
|||
'numpy',
|
||||
'qrcode[pil]',
|
||||
'scipy',
|
||||
'sklearn',
|
||||
'tornado',
|
||||
'voila >=0.2.4',
|
||||
],
|
||||
|
|
|
|||
Loading…
Reference in a new issue