Merge branch 'feature/piecewise_presence' into 'master'
Dynamic Number of Occupants (Infected) Closes #297 See merge request caimira/caimira!433
This commit is contained in:
Nicolas Mounet
commit
74322c46bb
11 changed files with 370 additions and 115 deletions
|
|
@ -35,7 +35,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 CAiMIRA version (found at ``caimira.__version__``).
|
||||
__version__ = "4.8"
|
||||
__version__ = "4.9"
|
||||
|
||||
|
||||
class BaseRequestHandler(RequestHandler):
|
||||
|
|
|
|||
|
|
@ -20,10 +20,10 @@ from ... import dataclass_utils
|
|||
|
||||
|
||||
def model_start_end(model: models.ExposureModel):
|
||||
t_start = min(model.exposed.presence.boundaries()[0][0],
|
||||
model.concentration_model.infected.presence.boundaries()[0][0])
|
||||
t_end = max(model.exposed.presence.boundaries()[-1][1],
|
||||
model.concentration_model.infected.presence.boundaries()[-1][1])
|
||||
t_start = min(model.exposed.presence_interval().boundaries()[0][0],
|
||||
model.concentration_model.infected.presence_interval().boundaries()[0][0])
|
||||
t_end = max(model.exposed.presence_interval().boundaries()[-1][1],
|
||||
model.concentration_model.infected.presence_interval().boundaries()[-1][1])
|
||||
return t_start, t_end
|
||||
|
||||
|
||||
|
|
@ -137,15 +137,16 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing
|
|||
prob = np.array(model.infection_probability())
|
||||
prob_dist_count, prob_dist_bins = np.histogram(prob/100, bins=100, density=True)
|
||||
prob_probabilistic_exposure = np.array(model.total_probability_rule()).mean()
|
||||
er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean()
|
||||
er = np.array(model.concentration_model.infected.emission_rate_per_person_when_present()).mean()
|
||||
exposed_occupants = model.exposed.number
|
||||
expected_new_cases = np.array(model.expected_new_cases()).mean()
|
||||
uncertainties_plot_src = img2base64(_figure2bytes(uncertainties_plot(model))) if form.conditional_probability_plot else None
|
||||
exposed_presence_intervals = [list(interval) for interval in model.exposed.presence_interval().boundaries()]
|
||||
|
||||
return {
|
||||
"model_repr": repr(model),
|
||||
"times": list(times),
|
||||
"exposed_presence_intervals": [list(interval) for interval in model.exposed.presence.boundaries()],
|
||||
"exposed_presence_intervals": exposed_presence_intervals,
|
||||
"short_range_intervals": short_range_intervals,
|
||||
"short_range_expirations": short_range_expirations,
|
||||
"concentrations": concentrations,
|
||||
|
|
@ -154,12 +155,12 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing
|
|||
"cumulative_doses": list(cumulative_doses),
|
||||
"long_range_cumulative_doses": list(long_range_cumulative_doses),
|
||||
"prob_inf": prob.mean(),
|
||||
"prob_inf_sd": np.std(prob),
|
||||
"prob_inf_sd": prob.std(),
|
||||
"prob_dist": list(prob),
|
||||
"prob_hist_count": list(prob_dist_count),
|
||||
"prob_hist_bins": list(prob_dist_bins),
|
||||
"prob_probabilistic_exposure": prob_probabilistic_exposure,
|
||||
"emission_rate": er,
|
||||
"emission_rate_per_person": er,
|
||||
"exposed_occupants": exposed_occupants,
|
||||
"expected_new_cases": expected_new_cases,
|
||||
"uncertainties_plot_src": uncertainties_plot_src,
|
||||
|
|
|
|||
|
|
@ -148,8 +148,9 @@ class ExposureModelResult(View):
|
|||
|
||||
def update_plot(self, model: models.ExposureModel):
|
||||
resolution = 600
|
||||
ts = np.linspace(sorted(model.concentration_model.infected.presence.transition_times())[0],
|
||||
sorted(model.concentration_model.infected.presence.transition_times())[-1], resolution)
|
||||
infected_presence = model.concentration_model.infected.presence_interval()
|
||||
ts = np.linspace(sorted(infected_presence.transition_times())[0],
|
||||
sorted(infected_presence.transition_times())[-1], resolution)
|
||||
concentration = [model.concentration(t) for t in ts]
|
||||
|
||||
cumulative_doses = np.cumsum([
|
||||
|
|
@ -164,16 +165,18 @@ class ExposureModelResult(View):
|
|||
self.ax.ignore_existing_data_limits = False
|
||||
self.concentration_line.set_data(ts, concentration)
|
||||
|
||||
exposed_presence = model.exposed.presence_interval()
|
||||
|
||||
if self.concentration_area is None:
|
||||
self.concentration_area = self.ax.fill_between(x = ts, y1=0, y2=concentration, color="#96cbff",
|
||||
where = ((model.exposed.presence.boundaries()[0][0] < ts) & (ts < model.exposed.presence.boundaries()[0][1]) |
|
||||
(model.exposed.presence.boundaries()[1][0] < ts) & (ts < model.exposed.presence.boundaries()[1][1])))
|
||||
where = ((exposed_presence.boundaries()[0][0] < ts) & (ts < exposed_presence.boundaries()[0][1]) |
|
||||
(exposed_presence.boundaries()[1][0] < ts) & (ts < exposed_presence.boundaries()[1][1])))
|
||||
|
||||
else:
|
||||
self.concentration_area.remove()
|
||||
self.concentration_area = self.ax.fill_between(x = ts, y1=0, y2=concentration, color="#96cbff",
|
||||
where = ((model.exposed.presence.boundaries()[0][0] < ts) & (ts < model.exposed.presence.boundaries()[0][1]) |
|
||||
(model.exposed.presence.boundaries()[1][0] < ts) & (ts < model.exposed.presence.boundaries()[1][1])))
|
||||
where = ((exposed_presence.boundaries()[0][0] < ts) & (ts < exposed_presence.boundaries()[0][1]) |
|
||||
(exposed_presence.boundaries()[1][0] < ts) & (ts < exposed_presence.boundaries()[1][1])))
|
||||
|
||||
if self.cumulative_line is None:
|
||||
[self.cumulative_line] = self.ax2.plot(ts[:-1], cumulative_doses, color='#0000c8', linestyle='dotted')
|
||||
|
|
@ -187,8 +190,8 @@ class ExposureModelResult(View):
|
|||
cumulative_top = max(cumulative_doses)
|
||||
self.ax2.set_ylim(bottom=0., top=cumulative_top)
|
||||
|
||||
self.ax.set_xlim(left = min(min(model.concentration_model.infected.presence.boundaries()[0]), min(model.exposed.presence.boundaries()[0])),
|
||||
right = max(max(model.concentration_model.infected.presence.boundaries()[1]), max(model.exposed.presence.boundaries()[1])))
|
||||
self.ax.set_xlim(left = min(min(infected_presence.boundaries()[0]), min(exposed_presence.boundaries()[0])),
|
||||
right = max(max(infected_presence.boundaries()[1]), max(exposed_presence.boundaries()[1])))
|
||||
|
||||
figure_legends = [mlines.Line2D([], [], color='#3530fe', markersize=15, label='Mean concentration'),
|
||||
mlines.Line2D([], [], color='#0000c8', markersize=15, ls="dotted", label='Cumulative dose'),
|
||||
|
|
@ -200,7 +203,7 @@ class ExposureModelResult(View):
|
|||
def update_textual_result(self, model: models.ExposureModel):
|
||||
lines = []
|
||||
P = np.array(model.infection_probability()).mean()
|
||||
lines.append(f'Emission rate (virus/hr): {np.round(model.concentration_model.infected.emission_rate_when_present(),0)}')
|
||||
lines.append(f'Emission rate per infected person (virus/hr): {np.round(model.concentration_model.infected.emission_rate_per_person_when_present(),0)}')
|
||||
lines.append(f'Probability of infection: {np.round(P, 0)}%')
|
||||
|
||||
lines.append(f'Number of exposed: {model.exposed.number}')
|
||||
|
|
@ -649,21 +652,10 @@ class ModelWidgets(View):
|
|||
number.observe(on_exposed_number_change, names=['value'])
|
||||
|
||||
return widgets.HBox([widgets.Label('Number of exposed people in the room '), number], layout=widgets.Layout(justify_content='space-between'))
|
||||
|
||||
def generate_presence_widget(self, min, max, node):
|
||||
options = list(pd.date_range(min, max, freq="1min").strftime('%H:%M'))
|
||||
start_hour = float(node[0])
|
||||
end_hour = float(node[1])
|
||||
start_hour_datetime = datetime.time(hour = int(start_hour), minute=int(start_hour%1*60))
|
||||
end_hour_datetime = datetime.time(hour = int(end_hour), minute=int(end_hour%1*60))
|
||||
return widgets.SelectionRangeSlider(
|
||||
options=options,
|
||||
index=(options.index(str(start_hour_datetime)[:-3]), options.index(str(end_hour_datetime)[:-3])),
|
||||
)
|
||||
|
||||
def _build_exposed_presence(self, node):
|
||||
presence_start = self.generate_presence_widget(min='00:00', max='13:00', node=node.present_times[0])
|
||||
presence_finish = self.generate_presence_widget(min='13:00', max='23:59', node=node.present_times[1])
|
||||
presence_start = generate_presence_widget(min='00:00', max='13:00', node=node.present_times[0])
|
||||
presence_finish = generate_presence_widget(min='13:00', max='23:59', node=node.present_times[1])
|
||||
|
||||
def on_presence_start_change(change):
|
||||
new_value = tuple([int(time[:-3])+float(time[3:])/60 for time in change['new']])
|
||||
|
|
@ -719,8 +711,8 @@ class ModelWidgets(View):
|
|||
return widgets.HBox([widgets.Label("Viral load (copies/ml)"), viral_load_in_sputum], layout=widgets.Layout(justify_content='space-between'))
|
||||
|
||||
def _build_infected_presence(self, node, ventilation_node):
|
||||
presence_start = self.generate_presence_widget(min='00:00', max='13:00', node=node.present_times[0])
|
||||
presence_finish = self.generate_presence_widget(min='13:00', max='23:59', node=node.present_times[1])
|
||||
presence_start = generate_presence_widget(min='00:00', max='13:00', node=node.present_times[0])
|
||||
presence_finish = generate_presence_widget(min='13:00', max='23:59', node=node.present_times[1])
|
||||
|
||||
def on_presence_start_change(change):
|
||||
new_value = tuple([int(time[:-3])+float(time[3:])/60 for time in change['new']])
|
||||
|
|
@ -1119,6 +1111,18 @@ def models_start_end(models: typing.Sequence[models.ExposureModel]) -> typing.Tu
|
|||
Returns the earliest start and latest end time of a collection of ConcentrationModel objects
|
||||
|
||||
"""
|
||||
infected_start = min(model.concentration_model.infected.presence.boundaries()[0][0] for model in models)
|
||||
infected_finish = min(model.concentration_model.infected.presence.boundaries()[-1][1] for model in models)
|
||||
infected_start = min(model.concentration_model.infected.presence_interval().boundaries()[0][0] for model in models)
|
||||
infected_finish = min(model.concentration_model.infected.presence_interval().boundaries()[-1][1] for model in models)
|
||||
return infected_start, infected_finish
|
||||
|
||||
|
||||
def generate_presence_widget(min, max, node):
|
||||
options = list(pd.date_range(min, max, freq="1min").strftime('%H:%M'))
|
||||
start_hour = float(node[0])
|
||||
end_hour = float(node[1])
|
||||
start_hour_datetime = datetime.time(hour = int(start_hour), minute=int(start_hour%1*60))
|
||||
end_hour_datetime = datetime.time(hour = int(end_hour), minute=int(end_hour%1*60))
|
||||
return widgets.SelectionRangeSlider(
|
||||
options=options,
|
||||
index=(options.index(str(start_hour_datetime)[:-3]), options.index(str(end_hour_datetime)[:-3])),
|
||||
)
|
||||
|
|
|
|||
|
|
@ -8,7 +8,7 @@ import matplotlib
|
|||
import matplotlib.figure
|
||||
import matplotlib.lines as mlines
|
||||
import matplotlib.patches as patches
|
||||
from .expert import collapsible, ipympl_canvas, WidgetGroup, CAIMIRAStateBuilder
|
||||
from .expert import generate_presence_widget, collapsible, ipympl_canvas, WidgetGroup, CAIMIRAStateBuilder
|
||||
|
||||
|
||||
baseline_model = models.CO2ConcentrationModel(
|
||||
|
|
@ -86,8 +86,8 @@ class ExposureModelResult(View):
|
|||
|
||||
def update_plot(self, model: models.CO2ConcentrationModel):
|
||||
resolution = 600
|
||||
ts = np.linspace(sorted(model.CO2_emitters.presence.transition_times())[0],
|
||||
sorted(model.CO2_emitters.presence.transition_times())[-1], resolution)
|
||||
ts = np.linspace(sorted(model.CO2_emitters.presence_interval().transition_times())[0],
|
||||
sorted(model.CO2_emitters.presence_interval().transition_times())[-1], resolution)
|
||||
concentration = [model.concentration(t) for t in ts]
|
||||
|
||||
if self.concentration_line is None:
|
||||
|
|
@ -99,19 +99,19 @@ class ExposureModelResult(View):
|
|||
|
||||
if self.concentration_area is None:
|
||||
self.concentration_area = self.ax.fill_between(x = ts, y1=0, y2=concentration, color="#96cbff",
|
||||
where = ((model.CO2_emitters.presence.boundaries()[0][0] < ts) & (ts < model.CO2_emitters.presence.boundaries()[0][1]) |
|
||||
(model.CO2_emitters.presence.boundaries()[1][0] < ts) & (ts < model.CO2_emitters.presence.boundaries()[1][1])))
|
||||
where = ((model.CO2_emitters.presence_interval().boundaries()[0][0] < ts) & (ts < model.CO2_emitters.presence_interval().boundaries()[0][1]) |
|
||||
(model.CO2_emitters.presence_interval().boundaries()[1][0] < ts) & (ts < model.CO2_emitters.presence_interval().boundaries()[1][1])))
|
||||
|
||||
else:
|
||||
self.concentration_area.remove()
|
||||
self.concentration_area = self.ax.fill_between(x = ts, y1=0, y2=concentration, color="#96cbff",
|
||||
where = ((model.CO2_emitters.presence.boundaries()[0][0] < ts) & (ts < model.CO2_emitters.presence.boundaries()[0][1]) |
|
||||
(model.CO2_emitters.presence.boundaries()[1][0] < ts) & (ts < model.CO2_emitters.presence.boundaries()[1][1])))
|
||||
where = ((model.CO2_emitters.presence_interval().boundaries()[0][0] < ts) & (ts < model.CO2_emitters.presence_interval().boundaries()[0][1]) |
|
||||
(model.CO2_emitters.presence_interval().boundaries()[1][0] < ts) & (ts < model.CO2_emitters.presence_interval().boundaries()[1][1])))
|
||||
|
||||
concentration_top = max(np.array(concentration))
|
||||
self.ax.set_ylim(bottom=model.CO2_atmosphere_concentration * 0.9, top=concentration_top*1.1)
|
||||
self.ax.set_xlim(left = min(model.CO2_emitters.presence.boundaries()[0])*0.95,
|
||||
right = max(model.CO2_emitters.presence.boundaries()[1])*1.05)
|
||||
self.ax.set_xlim(left = min(model.CO2_emitters.presence_interval().boundaries()[0])*0.95,
|
||||
right = max(model.CO2_emitters.presence_interval().boundaries()[1])*1.05)
|
||||
|
||||
figure_legends = [mlines.Line2D([], [], color='#3530fe', markersize=15, label='CO₂ concentration'),
|
||||
mlines.Line2D([], [], color='salmon', markersize=15, label='Insufficient level', linestyle='--'),
|
||||
|
|
@ -122,7 +122,10 @@ class ExposureModelResult(View):
|
|||
self.ax.set_ylim(top=concentration_top*1.1)
|
||||
else:
|
||||
self.ax.set_ylim(top=1550)
|
||||
self.ax.hlines([800, 1500], xmin=min(model.CO2_emitters.presence.boundaries()[0])*0.95, xmax=max(model.CO2_emitters.presence.boundaries()[1])*1.05, colors=['limegreen', 'salmon'], linestyles='dashed')
|
||||
self.ax.hlines([800, 1500], xmin=min(model.CO2_emitters.presence_interval().boundaries()[0])*0.95,
|
||||
xmax=max(model.CO2_emitters.presence_interval().boundaries()[1])*1.05,
|
||||
colors=['limegreen', 'salmon'],
|
||||
linestyles='dashed')
|
||||
self.figure.canvas.draw()
|
||||
|
||||
|
||||
|
|
@ -368,20 +371,26 @@ class ModelWidgets(View):
|
|||
return widgets.HBox([widgets.Label('Number of people in the room '), number], layout=widgets.Layout(justify_content='space-between'))
|
||||
|
||||
def _build_population_presence(self, node, ventilation_node):
|
||||
presence_start = widgets.FloatRangeSlider(value = node.present_times[0], min = 8., max=13., step=0.1)
|
||||
presence_finish = widgets.FloatRangeSlider(value = node.present_times[1], min = 13., max=18., step=0.1)
|
||||
presence_start = generate_presence_widget(min='00:00', max='13:00', node=node.present_times[0])
|
||||
presence_finish = generate_presence_widget(min='13:00', max='23:59', node=node.present_times[1])
|
||||
|
||||
def on_presence_start_change(change):
|
||||
ventilation_node.active.start = change['new'][0] - ventilation_node.active.duration / 60
|
||||
node.present_times = (change['new'], presence_finish.value)
|
||||
new_value = tuple([int(time[:-3])+float(time[3:])/60 for time in change['new']])
|
||||
ventilation_node.active.start = new_value[0] - ventilation_node.active.duration / 60
|
||||
node.present_times = (new_value, node.present_times[1])
|
||||
|
||||
def on_presence_finish_change(change):
|
||||
node.present_times = (presence_start.value, change['new'])
|
||||
new_value = tuple([int(time[:-3])+float(time[3:])/60 for time in change['new']])
|
||||
node.present_times = (node.present_times[0], new_value)
|
||||
|
||||
presence_start.observe(on_presence_start_change, names=['value'])
|
||||
presence_finish.observe(on_presence_finish_change, names=['value'])
|
||||
|
||||
return widgets.HBox([widgets.Label('Population presence'), presence_start, presence_finish], layout = widgets.Layout(justify_content='space-between'))
|
||||
return widgets.VBox([
|
||||
widgets.Label('Exposed presence:'),
|
||||
widgets.HBox([widgets.Label('Morning:', layout=widgets.Layout(width='15%')), presence_start]),
|
||||
widgets.HBox([widgets.Label('Afternoon:', layout=widgets.Layout(width='15%')), presence_finish])
|
||||
])
|
||||
|
||||
def present(self):
|
||||
return self.widget
|
||||
|
|
@ -782,6 +791,6 @@ def models_start_end(models: typing.Sequence[models.CO2ConcentrationModel]) -> t
|
|||
Returns the earliest start and latest end time of a collection of v objects
|
||||
|
||||
"""
|
||||
emitters_start = min(model.CO2_emitters.presence.boundaries()[0][0] for model in models)
|
||||
emitters_finish = min(model.CO2_emitters.presence.boundaries()[-1][1] for model in models)
|
||||
emitters_start = min(model.CO2_emitters.presence_interval().boundaries()[0][0] for model in models)
|
||||
emitters_finish = min(model.CO2_emitters.presence_interval().boundaries()[-1][1] for model in models)
|
||||
return emitters_start, emitters_finish
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@ provided the sample size is large enough. Example of the MC integration over the
|
|||
|
||||
It is important to distinguish between 1) Monte-Carlo random variables (which are vectorised independently on its diameter-dependence) and 2) numerical Monte-Carlo integration for the diameter-dependence.
|
||||
Since the integral of the diameter-dependent variables are solved when computing the dose -- :math:`\mathrm{vD^{total}}` -- while performing some of the intermediate calculations,
|
||||
we normalize the results by *dividing* by the Monte-Carlo variables that are diameter-independent, so that they are not considered in the Monte-Carlo integration (e.g. the **viral load** parameter, or the result of the :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_when_present` method).
|
||||
we normalize the results by *dividing* by the Monte-Carlo variables that are diameter-independent, so that they are not considered in the Monte-Carlo integration (e.g. the **viral load** parameter, or the result of the :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_per_person_when_present` method).
|
||||
|
||||
Expiration
|
||||
==========
|
||||
|
|
@ -62,9 +62,9 @@ Note that :math:`D_{\mathrm{max}}` value will differ, depending on the type of e
|
|||
|
||||
In the code, for a given Expiration, we use different methods to perform the calculations *step-by-step*:
|
||||
|
||||
1. Calculate the non aerosol-dependent quantities in the emission rate, which is the multiplication of the diameter-**independent** variables: :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_when_present`. This corresponds to the :math:`\mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` part of the :math:`\mathrm{vR}(D)` equation.
|
||||
1. Calculate the non aerosol-dependent quantities in the emission rate per person infected, which is the multiplication of the diameter-**independent** variables: :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_per_person_when_present`. This corresponds to the :math:`\mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` part of the :math:`\mathrm{vR}(D)` equation.
|
||||
2. Calculate the diameter-**dependent** variable :meth:`caimira.models.InfectedPopulation.aerosols`, which is the result of :math:`E_{c,j}(D) = N_p(D) \cdot V_p(D) \cdot (1 − η_\mathrm{out}(D))` (in mL/(m\ :sup:`3` \.µm)), with :math:`N_p(D)` being the product of the BLO distribution by the scaling factor :math:`cn`. Note that this result is not integrated over the diameters at this stage, thus the units are still *'per aerosol diameter'*.
|
||||
3. Calculate the full emission rate, which is the multiplication of the two previous methods, and corresponds to :math:`\mathrm{vR(D)}`: :meth:`caimira.models._PopulationWithVirus.emission_rate_when_present`.
|
||||
3. Calculate the full emission rate (per person infected), which is the multiplication of the two previous methods, and corresponds to :math:`\mathrm{vR(D)}`: :meth:`caimira.models._PopulationWithVirus.emission_rate_per_person_when_present`.
|
||||
|
||||
Note that the diameter-dependence is kept at this stage. Since other parameters downstream in code are also diameter-dependent, the Monte-Carlo integration over the aerosol sizes is computed at the level of the dose :math:`\mathrm{vD^{total}}`.
|
||||
In case one would like to have intermediate results for emission rate, perform the Monte-Carlo integration of :math:`E_{c, j}^{\mathrm{total}}` and compute :math:`\mathrm{vR^{total}} =\mathrm{vl_{in}} \cdot E_{c, j}^{\mathrm{total}} \cdot \mathrm{BR_k}`.
|
||||
|
|
@ -97,8 +97,8 @@ the code computes first a normalized version of the concentration, i.e. divided
|
|||
|
||||
To summarize, we can split the concentration in two different formulations:
|
||||
|
||||
* Normalized concentration :meth:`caimira.models._ConcentrationModelBase._normed_concentration`: :math:`\mathrm{C_\mathrm{LR, normed}}(t, D)` that computes the concentration without including the emission rate.
|
||||
* Concentration :meth:`caimira.models._ConcentrationModelBase.concentration` : :math:`C_{\mathrm{LR}}(t, D) = \mathrm{C_\mathrm{LR, normed}}(t, D) \cdot \mathrm{vR}(D)`, where :math:`\mathrm{vR}(D)` is the result of the :meth:`caimira.models._PopulationWithVirus.emission_rate_when_present` method.
|
||||
* Normalized concentration :meth:`caimira.models._ConcentrationModelBase._normed_concentration`: :math:`\mathrm{C_\mathrm{LR, normed}}(t, D)` that computes the concentration without including the emission rate per person infected.
|
||||
* Concentration :meth:`caimira.models._ConcentrationModelBase.concentration` : :math:`C_{\mathrm{LR}}(t, D) = \mathrm{C_\mathrm{LR, normed}}(t, D) \cdot \mathrm{vR}(D)`, where :math:`\mathrm{vR}(D)` is the result of the :meth:`caimira.models._PopulationWithVirus.emission_rate_per_person_when_present` method.
|
||||
|
||||
Note that in order to get the total concentration value in this stage, the final result should be averaged over the particle diameters (i.e. Monte-Carlo integration over diameters, see above).
|
||||
For the calculator app report, the total concentration (MC integral over the diameter) is performed only when generating the plot.
|
||||
|
|
@ -106,8 +106,8 @@ Otherwise, the diameter-dependence continues until we compute the inhaled dose i
|
|||
|
||||
The following methods calculate the integrated concentration between two times. They are mostly used when calculating the **dose**:
|
||||
|
||||
* :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`, :math:`\mathrm{C_\mathrm{normed}}(D)` that returns the integrated long-range concentration of viruses in the air, between any two times, normalized by the emission rate. Note that this method performs the integral between any two times of the previously mentioned :meth:`caimira.models._ConcentrationModelBase._normed_concentration` method.
|
||||
* :meth:`caimira.models._ConcentrationModelBase.integrated_concentration`, :math:`C(D)`, that returns the same result as the previous one, but multiplied by the emission rate.
|
||||
* :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`, :math:`\mathrm{C_\mathrm{normed}}(D)` that returns the integrated long-range concentration of viruses in the air, between any two times, normalized by the emission rate per person infected. Note that this method performs the integral between any two times of the previously mentioned :meth:`caimira.models._ConcentrationModelBase._normed_concentration` method.
|
||||
* :meth:`caimira.models._ConcentrationModelBase.integrated_concentration`, :math:`C(D)`, that returns the same result as the previous one, but multiplied by the emission rate (per person infected).
|
||||
|
||||
The integral over the exposure times is calculated directly in the class (integrated methods).
|
||||
|
||||
|
|
@ -221,14 +221,14 @@ Long-range approach
|
|||
*******************
|
||||
|
||||
Regarding the concentration part of the long-range exposure (concentration integrated over time, :math:`\int_{t1}^{t2}C_{\mathrm{LR}}(t, D)\;\mathrm{d}t`), the respective method is :meth:`caimira.models.ExposureModel._long_range_normed_exposure_between_bounds`,
|
||||
which uses the long-range exposure (concentration) between two bounds (time1 and time2), normalized by the emission rate of the infected population, calculated from :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`.
|
||||
which uses the long-range exposure (concentration) between two bounds (time1 and time2), normalized by the emission rate of the infected population (per person infected), calculated from :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`.
|
||||
The former method filters out the given bounds considering the breaks through the day (i.e. the time intervals during which there is no exposition to the virus) and retrieves the integrated long-range concentration of viruses in the air between any two times.
|
||||
|
||||
After the calculations of the integrated concentration over the time, in order to calculate the final dose, we have to compute the remaining factors in the above equation.
|
||||
Note that the **Monte-Carlo integration over the diameters is performed at this stage**, where all the diameter-dependent parameters are grouped together to calculate the final average (:code:`np.mean()`).
|
||||
|
||||
Since, in the previous chapters, the quantities where normalised by the emission rate, one will need to re-incorporate it in the equations before performing the MC integrations over :math:`D`.
|
||||
For that we need to split :math:`\mathrm{vR}(D)` (:meth:`caimira.models._PopulationWithVirus.emission_rate_when_present`) in diameter-dependent and diameter-independent quantities:
|
||||
Since, in the previous chapters, the quantities where normalised by the emission rate per person infected, one will need to re-incorporate it in the equations before performing the MC integrations over :math:`D`.
|
||||
For that we need to split :math:`\mathrm{vR}(D)` (:meth:`caimira.models._PopulationWithVirus.emission_rate_per_person_when_present`) in diameter-dependent and diameter-independent quantities:
|
||||
|
||||
:math:`\mathrm{vR}(D) = \mathrm{vR}(D-\mathrm{dependent}) \times \mathrm{vR}(D-\mathrm{independent})`
|
||||
|
||||
|
|
@ -236,7 +236,7 @@ with
|
|||
|
||||
:math:`\mathrm{vR}(D-\mathrm{dependent}) = \mathrm{cn} \cdot V_p(D) \cdot (1 − \mathrm{η_{out}}(D))` - :meth:`caimira.models.InfectedPopulation.aerosols`
|
||||
|
||||
:math:`\mathrm{vR}(D-\mathrm{independent}) = \mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` - :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_when_present`
|
||||
:math:`\mathrm{vR}(D-\mathrm{independent}) = \mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` - :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_per_person_when_present`
|
||||
|
||||
|
||||
In other words, in the code the procedure is the following (all performed in :meth:`caimira.models.ExposureModel.long_range_deposited_exposure_between_bounds` method):
|
||||
|
|
@ -248,7 +248,7 @@ In other words, in the code the procedure is the following (all performed in :me
|
|||
* in order to complete the equation, multiply by the remaining diameter-independent variables in :math:`\mathrm{vD}` to obtain the total value: :math:`\mathrm{vD^{total}} = \mathrm{vD_{emission-rate}} \cdot \mathrm{BR}_{\mathrm{k}} \cdot (1-\eta_{\mathrm{in}}) \cdot f_{\mathrm{inf}}`;
|
||||
* in the end, the dose is a vectorized float used in the probability of infection formula.
|
||||
|
||||
**Note**: The aerosol volume concentration (*aerosols*) is introduced because the integrated concentration over the time was previously normalized by the emission rate.
|
||||
**Note**: The aerosol volume concentration (*aerosols*) is introduced because the integrated concentration over the time was previously normalized by the emission rate (per person).
|
||||
Here, to calculate the integral over the diameters we also need to consider the diameter-dependent variables that are on the emission rate, represented by the aerosol volume concentration which depends on the diameter and on the mask type:
|
||||
|
||||
:math:`\mathrm{aerosols} = \mathrm{cn} \cdot V_p(D) \cdot (1 − \mathrm{η_{out}}(D))` .
|
||||
|
|
@ -339,4 +339,4 @@ REFERENCES
|
|||
==========
|
||||
|
||||
.. [1] Jia, Wei, et al. "Exposure and respiratory infection risk via the short-range airborne route." Building and environment 219 (2022): 109166. `doi.org/10.1016/j.buildenv.2022.109166 <https://doi.org/10.1016/j.buildenv.2022.109166>`_
|
||||
.. [2] Henriques, Andre, et al. "Modelling airborne transmission of SARS-CoV-2 using CARA: risk assessment for enclosed spaces." Interface Focus 12.2 (2022): 20210076. `doi.org/10.1098/rsfs.2021.0076 <https://doi.org/10.1098/rsfs.2021.0076>`_
|
||||
.. [2] Henriques, Andre, et al. "Modelling airborne transmission of SARS-CoV-2 using CARA: risk assessment for enclosed spaces." Interface Focus 12.2 (2022): 20210076. `doi.org/10.1098/rsfs.2021.0076 <https://doi.org/10.1098/rsfs.2021.0076>`_
|
||||
|
|
|
|||
|
|
@ -186,6 +186,23 @@ class PiecewiseConstant:
|
|||
)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class IntPiecewiseConstant(PiecewiseConstant):
|
||||
|
||||
#: values of the function between transitions
|
||||
values: typing.Tuple[int, ...]
|
||||
|
||||
def value(self, time) -> _VectorisedFloat:
|
||||
if time <= self.transition_times[0] or time > self.transition_times[-1]:
|
||||
return 0
|
||||
|
||||
for t1, t2, value in zip(self.transition_times[:-1],
|
||||
self.transition_times[1:], self.values):
|
||||
if t1 < time <= t2:
|
||||
break
|
||||
return value
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Room:
|
||||
#: The total volume of the room
|
||||
|
|
@ -779,10 +796,10 @@ class Population:
|
|||
|
||||
"""
|
||||
#: How many in the population.
|
||||
number: int
|
||||
number: typing.Union[int, IntPiecewiseConstant]
|
||||
|
||||
#: The times in which the people are in the room.
|
||||
presence: Interval
|
||||
presence: typing.Union[None, Interval]
|
||||
|
||||
#: The kind of mask being worn by the people.
|
||||
mask: Mask
|
||||
|
|
@ -795,8 +812,33 @@ class Population:
|
|||
# which might render inactive some RNA copies (virions).
|
||||
host_immunity: float
|
||||
|
||||
def person_present(self, time):
|
||||
return self.presence.triggered(time)
|
||||
def __post_init__(self):
|
||||
if isinstance(self.number, int):
|
||||
if not isinstance(self.presence, Interval):
|
||||
raise TypeError(f'The presence argument must be an "Interval". Got {type(self.presence)}')
|
||||
else:
|
||||
if self.presence is not None:
|
||||
raise TypeError(f'The presence argument must be None for a IntPiecewiseConstant number')
|
||||
|
||||
def presence_interval(self):
|
||||
if isinstance(self.presence, Interval):
|
||||
return self.presence
|
||||
elif isinstance(self.number, IntPiecewiseConstant):
|
||||
return self.number.interval()
|
||||
|
||||
def person_present(self, time: float):
|
||||
# Allow back-compatibility
|
||||
if isinstance(self.number, int) and isinstance(self.presence, Interval):
|
||||
return self.presence.triggered(time)
|
||||
elif isinstance(self.number, IntPiecewiseConstant):
|
||||
return self.number.value(time) != 0
|
||||
|
||||
def people_present(self, time: float):
|
||||
# Allow back-compatibility
|
||||
if isinstance(self.number, int):
|
||||
return self.number * self.person_present(time)
|
||||
else:
|
||||
return int(self.number.value(time))
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
|
|
@ -818,22 +860,22 @@ class _PopulationWithVirus(Population):
|
|||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
|
||||
def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat:
|
||||
def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate of virions in the expired air per mL of respiratory fluid,
|
||||
if the infected population is present, in (virion.cm^3)/(mL.h).
|
||||
per person, if the infected population is present, in (virion.cm^3)/(mL.h).
|
||||
This method includes only the diameter-independent variables within the emission rate.
|
||||
It should not be a function of time.
|
||||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
|
||||
@method_cache
|
||||
def emission_rate_when_present(self) -> _VectorisedFloat:
|
||||
def emission_rate_per_person_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate if the infected population is present
|
||||
The emission rate if the infected population is present, per person
|
||||
(in virions / h).
|
||||
"""
|
||||
return (self.emission_rate_per_aerosol_when_present() *
|
||||
return (self.emission_rate_per_aerosol_per_person_when_present() *
|
||||
self.aerosols())
|
||||
|
||||
def emission_rate(self, time) -> _VectorisedFloat:
|
||||
|
|
@ -850,8 +892,7 @@ class _PopulationWithVirus(Population):
|
|||
# itself a function of time. Any change in rate must be accompanied
|
||||
# with a declaration of state change time, as is the case for things
|
||||
# like Ventilation.
|
||||
|
||||
return self.emission_rate_when_present()
|
||||
return self.emission_rate_per_person_when_present() * self.people_present(time)
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
|
|
@ -875,14 +916,14 @@ class EmittingPopulation(_PopulationWithVirus):
|
|||
return 1.
|
||||
|
||||
@method_cache
|
||||
def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat:
|
||||
def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate of virions in the expired air per mL of respiratory fluid,
|
||||
if the infected population is present, in (virion.cm^3)/(mL.h).
|
||||
per person, if the infected population is present, in (virion.cm^3)/(mL.h).
|
||||
This method includes only the diameter-independent variables within the emission rate.
|
||||
It should not be a function of time.
|
||||
"""
|
||||
return self.known_individual_emission_rate * self.number
|
||||
return self.known_individual_emission_rate
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
|
|
@ -904,7 +945,7 @@ class InfectedPopulation(_PopulationWithVirus):
|
|||
return self.expiration.aerosols(self.mask)
|
||||
|
||||
@method_cache
|
||||
def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat:
|
||||
def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate of virions in the expired air per mL of respiratory fluid,
|
||||
if the infected population is present, in (virion.cm^3)/(mL.h).
|
||||
|
|
@ -917,8 +958,7 @@ class InfectedPopulation(_PopulationWithVirus):
|
|||
ER = (self.virus.viral_load_in_sputum *
|
||||
self.activity.exhalation_rate *
|
||||
10 ** 6)
|
||||
|
||||
return ER * self.number
|
||||
return ER
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
|
|
@ -1005,13 +1045,11 @@ class _ConcentrationModelBase:
|
|||
can be put back in front of the concentration after the time
|
||||
dependence has been solved for.
|
||||
"""
|
||||
if not self.population.person_present(time):
|
||||
return self.min_background_concentration()/self.normalization_factor()
|
||||
V = self.room.volume
|
||||
RR = self.removal_rate(time)
|
||||
|
||||
return (1. / (RR * V) + self.min_background_concentration()/
|
||||
self.normalization_factor())
|
||||
return (self.population.people_present(time) / (RR * V) +
|
||||
self.min_background_concentration()/self.normalization_factor())
|
||||
|
||||
@method_cache
|
||||
def state_change_times(self) -> typing.List[float]:
|
||||
|
|
@ -1020,8 +1058,9 @@ class _ConcentrationModelBase:
|
|||
the times at which their state changes.
|
||||
"""
|
||||
state_change_times = {0.}
|
||||
state_change_times.update(self.population.presence.transition_times())
|
||||
state_change_times.update(self.population.presence_interval().transition_times())
|
||||
state_change_times.update(self.ventilation.transition_times(self.room))
|
||||
|
||||
return sorted(state_change_times)
|
||||
|
||||
@method_cache
|
||||
|
|
@ -1029,8 +1068,8 @@ class _ConcentrationModelBase:
|
|||
"""
|
||||
First presence time. Before that, the concentration is zero.
|
||||
"""
|
||||
return self.population.presence.boundaries()[0][0]
|
||||
|
||||
return self.population.presence_interval().boundaries()[0][0]
|
||||
|
||||
def last_state_change(self, time: float) -> float:
|
||||
"""
|
||||
Find the most recent/previous state change.
|
||||
|
|
@ -1083,7 +1122,9 @@ class _ConcentrationModelBase:
|
|||
# before the first presence as an optimisation.
|
||||
if time <= self._first_presence_time():
|
||||
return self.min_background_concentration()/self.normalization_factor()
|
||||
|
||||
next_state_change_time = self._next_state_change(time)
|
||||
|
||||
RR = self.removal_rate(next_state_change_time)
|
||||
# If RR is 0, conc_limit does not play a role but its computation
|
||||
# would raise an error -> we set it to zero.
|
||||
|
|
@ -1097,6 +1138,7 @@ class _ConcentrationModelBase:
|
|||
|
||||
delta_time = time - t_last_state_change
|
||||
fac = np.exp(-RR * delta_time)
|
||||
|
||||
return conc_limit * (1 - fac) + conc_at_last_state_change * fac
|
||||
|
||||
def concentration(self, time: float) -> _VectorisedFloat:
|
||||
|
|
@ -1172,7 +1214,7 @@ class ConcentrationModel(_ConcentrationModelBase):
|
|||
|
||||
def normalization_factor(self) -> _VectorisedFloat:
|
||||
# we normalize by the emission rate
|
||||
return self.infected.emission_rate_when_present()
|
||||
return self.infected.emission_rate_per_person_when_present()
|
||||
|
||||
def removal_rate(self, time: float) -> _VectorisedFloat:
|
||||
# Equilibrium velocity of particle motion toward the floor
|
||||
|
|
@ -1221,10 +1263,9 @@ class CO2ConcentrationModel(_ConcentrationModelBase):
|
|||
return self.CO2_atmosphere_concentration
|
||||
|
||||
def normalization_factor(self) -> _VectorisedFloat:
|
||||
# normalization by the CO2 exhaled.
|
||||
# normalization by the CO2 exhaled per person.
|
||||
# CO2 concentration given in ppm, hence the 1e6 factor.
|
||||
return (1e6*self.population.number
|
||||
*self.population.activity.exhalation_rate
|
||||
return (1e6*self.population.activity.exhalation_rate
|
||||
*self.CO2_fraction_exhaled)
|
||||
|
||||
|
||||
|
|
@ -1451,8 +1492,11 @@ class ExposureModel:
|
|||
and not (
|
||||
all(np.isscalar(c_model.virus.decay_constant(c_model.room.humidity, c_model.room.inside_temp.value(time)) +
|
||||
c_model.ventilation.air_exchange(c_model.room, time)) for time in c_model.state_change_times()))):
|
||||
raise ValueError("If the diameter is an array, none of the ventilation parameters "
|
||||
"or virus decay constant can be arrays at the same time.")
|
||||
raise ValueError("If the diameter is an array, none of the ventilation parameters "
|
||||
"or virus decay constant can be arrays at the same time.")
|
||||
if not isinstance(self.exposed.number, int):
|
||||
raise NotImplementedError("Cannot use dynamic occupancy for"
|
||||
" the exposed population")
|
||||
|
||||
def long_range_fraction_deposited(self) -> _VectorisedFloat:
|
||||
"""
|
||||
|
|
@ -1469,7 +1513,7 @@ class ExposureModel:
|
|||
by the emission rate of the infected population
|
||||
"""
|
||||
exposure = 0.
|
||||
for start, stop in self.exposed.presence.boundaries():
|
||||
for start, stop in self.exposed.presence_interval().boundaries():
|
||||
if stop < time1:
|
||||
continue
|
||||
elif start > time2:
|
||||
|
|
@ -1499,7 +1543,8 @@ class ExposureModel:
|
|||
def long_range_deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
|
||||
deposited_exposure = 0.
|
||||
|
||||
emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present()
|
||||
emission_rate_per_aerosol_per_person = \
|
||||
self.concentration_model.infected.emission_rate_per_aerosol_per_person_when_present()
|
||||
aerosols = self.concentration_model.infected.aerosols()
|
||||
f_inf = self.concentration_model.infected.fraction_of_infectious_virus()
|
||||
fdep = self.long_range_fraction_deposited()
|
||||
|
|
@ -1519,10 +1564,11 @@ class ExposureModel:
|
|||
# one should not take any mean at this stage.
|
||||
dep_exposure_integrated = self._long_range_normed_exposure_between_bounds(time1, time2)*aerosols*fdep
|
||||
|
||||
# Then we multiply by the diameter-independent quantity emission_rate_per_aerosol,
|
||||
# Then we multiply by the diameter-independent quantity emission_rate_per_aerosol_per_person,
|
||||
# and parameters of the vD equation (i.e. BR_k and n_in).
|
||||
deposited_exposure += (dep_exposure_integrated * emission_rate_per_aerosol *
|
||||
self.exposed.activity.inhalation_rate *
|
||||
deposited_exposure += (dep_exposure_integrated *
|
||||
emission_rate_per_aerosol_per_person *
|
||||
self.exposed.activity.inhalation_rate *
|
||||
(1 - self.exposed.mask.inhale_efficiency()))
|
||||
|
||||
# In the end we multiply the final results by the fraction of infectious virus of the vD equation.
|
||||
|
|
@ -1589,8 +1635,7 @@ class ExposureModel:
|
|||
The number of virus per m^3 deposited on the respiratory tract.
|
||||
"""
|
||||
deposited_exposure: _VectorisedFloat = 0.0
|
||||
|
||||
for start, stop in self.exposed.presence.boundaries():
|
||||
for start, stop in self.exposed.presence_interval().boundaries():
|
||||
deposited_exposure += self.deposited_exposure_between_bounds(start, stop)
|
||||
|
||||
return deposited_exposure * self.repeats
|
||||
|
|
@ -1607,11 +1652,17 @@ class ExposureModel:
|
|||
self.concentration_model.virus.transmissibility_factor)))) * 100
|
||||
|
||||
def total_probability_rule(self) -> _VectorisedFloat:
|
||||
if (isinstance(self.concentration_model.infected.number, IntPiecewiseConstant) or
|
||||
isinstance(self.exposed.number, IntPiecewiseConstant)):
|
||||
raise NotImplementedError("Cannot compute total probability "
|
||||
"(including incidence rate) with dynamic occupancy")
|
||||
|
||||
if (self.geographical_data.geographic_population != 0 and self.geographical_data.geographic_cases != 0):
|
||||
sum_probability = 0.0
|
||||
|
||||
# Create an equivalent exposure model but changing the number of infected cases.
|
||||
total_people = self.concentration_model.infected.number + self.exposed.number
|
||||
max_num_infected = (total_people if total_people < 10 else 10)
|
||||
max_num_infected = (total_people if total_people < 10 else 10)
|
||||
# The influence of a higher number of simultainious infected people (> 4 - 5) yields an almost negligible contirbution to the total probability.
|
||||
# To be on the safe side, a hard coded limit with a safety margin of 2x was set.
|
||||
# Therefore we decided a hard limit of 10 infected people.
|
||||
|
|
|
|||
|
|
@ -74,6 +74,14 @@ def _build_mc_model(model: dataclass_instance) -> typing.Type[MCModelBase[_Model
|
|||
elif new_field.type == typing.Tuple[caimira.models.SpecificInterval, ...]:
|
||||
SI = getattr(sys.modules[__name__], "SpecificInterval")
|
||||
field_type = typing.Tuple[typing.Union[caimira.models.SpecificInterval, SI], ...]
|
||||
|
||||
elif new_field.type == typing.Union[int, caimira.models.IntPiecewiseConstant]:
|
||||
IPC = getattr(sys.modules[__name__], "IntPiecewiseConstant")
|
||||
field_type = typing.Union[int, caimira.models.IntPiecewiseConstant, IPC]
|
||||
elif new_field.type == typing.Union[caimira.models.Interval, None]:
|
||||
I = getattr(sys.modules[__name__], "Interval")
|
||||
field_type = typing.Union[None, caimira.models.Interval, I]
|
||||
|
||||
else:
|
||||
# Check that we don't need to do anything with this type.
|
||||
for item in new_field.type.__args__:
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@ import numpy as np
|
|||
import numpy.testing as npt
|
||||
import pytest
|
||||
from dataclasses import dataclass
|
||||
import typing
|
||||
|
||||
from caimira import models
|
||||
|
||||
|
|
@ -106,7 +107,7 @@ def simple_conc_model():
|
|||
@pytest.fixture
|
||||
def dummy_population(simple_conc_model) -> models.Population:
|
||||
return models.Population(
|
||||
number=10,
|
||||
number=1,
|
||||
presence=simple_conc_model.infected.presence,
|
||||
mask=models.Mask.types['Type I'],
|
||||
activity=models.Activity.types['Seated'],
|
||||
|
|
@ -267,7 +268,7 @@ def test_zero_ventilation_rate(
|
|||
ventilation = simple_conc_model.ventilation,
|
||||
known_population = dummy_population,
|
||||
known_removal_rate = known_removal_rate,
|
||||
known_normalization_factor=10.,
|
||||
known_normalization_factor=1.,
|
||||
known_min_background_concentration = known_min_background_concentration)
|
||||
|
||||
normed_concentration = known_conc_model.concentration(1)
|
||||
|
|
|
|||
180
caimira/tests/models/test_dynamic_population.py
Normal file
180
caimira/tests/models/test_dynamic_population.py
Normal file
|
|
@ -0,0 +1,180 @@
|
|||
import re
|
||||
|
||||
import numpy as np
|
||||
import numpy.testing as npt
|
||||
import pytest
|
||||
|
||||
from caimira import models
|
||||
import caimira.dataclass_utils as dc_utils
|
||||
|
||||
@pytest.fixture
|
||||
def full_exposure_model():
|
||||
return models.ExposureModel(
|
||||
concentration_model=models.ConcentrationModel(
|
||||
room=models.Room(volume=100),
|
||||
ventilation=models.AirChange(
|
||||
active=models.PeriodicInterval(120, 120), air_exch=0.25),
|
||||
infected=models.InfectedPopulation(
|
||||
number=1,
|
||||
presence=models.SpecificInterval(((8, 12), (13, 17), )),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Seated'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
virus=models.Virus.types['SARS_CoV_2'],
|
||||
host_immunity=0.
|
||||
),
|
||||
),
|
||||
short_range=(),
|
||||
exposed=models.Population(
|
||||
number=10,
|
||||
presence=models.SpecificInterval(((8, 12), (13, 17), )),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Seated'],
|
||||
host_immunity=0.
|
||||
),
|
||||
geographical_data=(),
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def baseline_infected_population_number():
|
||||
return models.InfectedPopulation(
|
||||
number=models.IntPiecewiseConstant(
|
||||
(8, 12, 13, 17), (1, 0, 1)),
|
||||
presence=None,
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Seated'],
|
||||
virus=models.Virus.types['SARS_CoV_2'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
host_immunity=0.,
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def dynamic_single_exposure_model(full_exposure_model, baseline_infected_population_number):
|
||||
return dc_utils.nested_replace(full_exposure_model,
|
||||
{'concentration_model.infected': baseline_infected_population_number, })
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"time",
|
||||
[4., 8., 10., 12., 13., 14., 16., 20., 24.],
|
||||
)
|
||||
def test_population_number(full_exposure_model: models.ExposureModel,
|
||||
baseline_infected_population_number: models.InfectedPopulation, time: float):
|
||||
|
||||
int_population_number: models.InfectedPopulation = full_exposure_model.concentration_model.infected
|
||||
piecewise_population_number: models.InfectedPopulation = baseline_infected_population_number
|
||||
|
||||
with pytest.raises(
|
||||
TypeError,
|
||||
match=f'The presence argument must be an "Interval". Got {type(None)}'
|
||||
):
|
||||
dc_utils.nested_replace(
|
||||
int_population_number, {'presence': None}
|
||||
)
|
||||
|
||||
with pytest.raises(
|
||||
TypeError,
|
||||
match="The presence argument must be None for a IntPiecewiseConstant number"
|
||||
):
|
||||
dc_utils.nested_replace(
|
||||
piecewise_population_number, {'presence': models.SpecificInterval(((8, 12), ))}
|
||||
)
|
||||
|
||||
assert int_population_number.person_present(time) == piecewise_population_number.person_present(time)
|
||||
assert int_population_number.people_present(time) == piecewise_population_number.people_present(time)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"time",
|
||||
[4., 8., 10., 12., 13., 14., 16., 20., 24.],
|
||||
)
|
||||
def test_concentration_model_dynamic_population(full_exposure_model: models.ExposureModel,
|
||||
dynamic_single_exposure_model: models.ExposureModel,
|
||||
time: float):
|
||||
|
||||
assert full_exposure_model.concentration(time) == dynamic_single_exposure_model.concentration(time)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("number_of_infected",[1, 2, 3, 4, 5])
|
||||
@pytest.mark.parametrize("time",[9., 12.5, 16.])
|
||||
def test_linearity_with_number_of_infected(full_exposure_model: models.ExposureModel,
|
||||
dynamic_single_exposure_model: models.ExposureModel,
|
||||
time: float,
|
||||
number_of_infected: int):
|
||||
|
||||
|
||||
static_multiple_exposure_model: models.ExposureModel = dc_utils.nested_replace(
|
||||
full_exposure_model,
|
||||
{
|
||||
'concentration_model.infected.number': number_of_infected,
|
||||
}
|
||||
)
|
||||
|
||||
npt.assert_almost_equal(static_multiple_exposure_model.concentration(time), dynamic_single_exposure_model.concentration(time) * number_of_infected)
|
||||
npt.assert_almost_equal(static_multiple_exposure_model.deposited_exposure(), dynamic_single_exposure_model.deposited_exposure() * number_of_infected)
|
||||
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"time", (8., 9., 10., 11., 12., 13., 14.),
|
||||
)
|
||||
def test_dynamic_dose(full_exposure_model, time):
|
||||
|
||||
dynamic_infected: models.ExposureModel = dc_utils.nested_replace(
|
||||
full_exposure_model,
|
||||
{
|
||||
'concentration_model.infected': models.InfectedPopulation(
|
||||
number=models.IntPiecewiseConstant(
|
||||
(8, 10, 12, 13, 17), (1, 2, 0, 3)),
|
||||
presence=None,
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=models.Activity.types['Seated'],
|
||||
virus=models.Virus.types['SARS_CoV_2'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
}
|
||||
)
|
||||
|
||||
single_infected: models.ExposureModel = dc_utils.nested_replace(
|
||||
full_exposure_model,
|
||||
{
|
||||
'concentration_model.infected.number': 1,
|
||||
'concentration_model.infected.presence': models.SpecificInterval(((8, 10), )),
|
||||
}
|
||||
)
|
||||
|
||||
two_infected: models.ExposureModel = dc_utils.nested_replace(
|
||||
full_exposure_model,
|
||||
{
|
||||
'concentration_model.infected.number': 2,
|
||||
'concentration_model.infected.presence': models.SpecificInterval(((10, 12), )),
|
||||
}
|
||||
)
|
||||
|
||||
three_infected: models.ExposureModel = dc_utils.nested_replace(
|
||||
full_exposure_model,
|
||||
{
|
||||
'concentration_model.infected.number': 3,
|
||||
'concentration_model.infected.presence': models.SpecificInterval(((13, 17), )),
|
||||
}
|
||||
)
|
||||
|
||||
dynamic_concentration = dynamic_infected.concentration(time)
|
||||
dynamic_exposure = dynamic_infected.deposited_exposure()
|
||||
|
||||
static_concentration, static_exposure = zip(*[(model.concentration(time), model.deposited_exposure())
|
||||
for model in (single_infected, two_infected, three_infected)])
|
||||
|
||||
npt.assert_almost_equal(dynamic_concentration, np.sum(static_concentration))
|
||||
npt.assert_almost_equal(dynamic_exposure, np.sum(static_exposure))
|
||||
|
||||
|
||||
def test_dynamic_total_probability_rule(dynamic_single_exposure_model: models.ExposureModel):
|
||||
with pytest.raises(
|
||||
NotImplementedError,
|
||||
match=re.escape("Cannot compute total probability "
|
||||
"(including incidence rate) with dynamic occupancy")
|
||||
):
|
||||
dynamic_single_exposure_model.total_probability_rule()
|
||||
|
|
@ -24,7 +24,7 @@ class KnownNormedconcentration(models.ConcentrationModel):
|
|||
return 1.e50
|
||||
|
||||
def _normed_concentration_limit(self, time: float) -> models._VectorisedFloat:
|
||||
return self.normed_concentration_function(time)
|
||||
return self.normed_concentration_function(time) * self.infected.number
|
||||
|
||||
def state_change_times(self):
|
||||
return [0., 24.]
|
||||
|
|
@ -33,7 +33,7 @@ class KnownNormedconcentration(models.ConcentrationModel):
|
|||
return 24.
|
||||
|
||||
def _normed_concentration(self, time: float) -> models._VectorisedFloat: # noqa
|
||||
return self.normed_concentration_function(time)
|
||||
return self.normed_concentration_function(time) * self.infected.number
|
||||
|
||||
|
||||
halftime = models.PeriodicInterval(120, 60)
|
||||
|
|
@ -67,7 +67,8 @@ def known_concentrations(func):
|
|||
expiration=models.Expiration.types['Speaking'],
|
||||
host_immunity=0.,
|
||||
)
|
||||
normed_func = lambda x: func(x) / dummy_infected_population.emission_rate_when_present()
|
||||
normed_func = lambda x: (func(x) /
|
||||
dummy_infected_population.emission_rate_per_person_when_present())
|
||||
return KnownNormedconcentration(dummy_room, dummy_ventilation,
|
||||
dummy_infected_population, 0.3, normed_func)
|
||||
|
||||
|
|
|
|||
|
|
@ -313,19 +313,19 @@ def waiting_room_mc():
|
|||
|
||||
@retry(tries=10)
|
||||
@pytest.mark.parametrize(
|
||||
"mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER",
|
||||
"mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER_per_person",
|
||||
[
|
||||
["shared_office_mc", 5.38, 0.16, 3.350, 1056],
|
||||
["classroom_mc", 8.21, 1.56, 11.356, 7416],
|
||||
["ski_cabin_mc", 12.92, 0.39, 21.796, 10231],
|
||||
["skagit_chorale_mc",61.01, 36.53, 84.730, 190422],
|
||||
["bus_ride_mc", 10.59, 7.06, 6.650, 5419],
|
||||
["gym_mc", 0.52, 0.14, 0.249, 1450],
|
||||
["gym_mc", 0.52, 0.14, 0.249, 1450/2.], # there are two infected in this case
|
||||
["waiting_room_mc", 1.53, 0.21, 0.844, 929],
|
||||
]
|
||||
)
|
||||
def test_report_models(mc_model, expected_pi, expected_new_cases,
|
||||
expected_dose, expected_ER, request):
|
||||
expected_dose, expected_ER_per_person, request):
|
||||
mc_model = request.getfixturevalue(mc_model)
|
||||
exposure_model = mc_model.build_model(size=SAMPLE_SIZE)
|
||||
npt.assert_allclose(exposure_model.infection_probability().mean(),
|
||||
|
|
@ -335,8 +335,8 @@ def test_report_models(mc_model, expected_pi, expected_new_cases,
|
|||
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(),
|
||||
expected_ER, rtol=TOLERANCE)
|
||||
exposure_model.concentration_model.infected.emission_rate_per_person_when_present().mean(),
|
||||
expected_ER_per_person, rtol=TOLERANCE)
|
||||
|
||||
|
||||
@retry(tries=10)
|
||||
|
|
@ -397,5 +397,5 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi,
|
|||
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(),
|
||||
exposure_model.concentration_model.infected.emission_rate_per_person_when_present().mean(),
|
||||
expected_ER, rtol=TOLERANCE)
|
||||
|
|
|
|||
Loading…
Reference in a new issue