Commit 31d2e220 authored by Philip Elson's avatar Philip Elson 🐍
Browse files

Add infection probability calculation to the model, and introduce it into the...

Add infection probability calculation to the model, and introduce it into the model output in a rudimentary form.
parent 568f7ca5
......@@ -77,7 +77,6 @@ class WidgetView:
self.widget = widgets.VBox([])
self.widgets = {}
self.out = widgets.Output()
self.widget.children += (self.out, )
self.plots = []
self.construct_widgets()
# Trigger the first result.
......@@ -95,6 +94,7 @@ class WidgetView:
self.widgets['results'] = collapsible([
widgets.HBox([
concentration.figure.canvas,
self.out,
])
], 'Results', start_collapsed=False)
......@@ -109,6 +109,17 @@ class WidgetView:
for plot in self.plots:
plot.update(model)
self.out.clear_output()
with self.out:
P = model.infection_probability()
print(f'Emission rate (quanta/hr): {model.infected.emission_rate(0)}')
print(f'Probability of infection: {np.round(P, 0)}%')
print(f'Number of exposed: {model.exposed_occupants}')
R0 = np.round(P / 100 * model.exposed_occupants, 1)
print(f'Number of expected new cases (R0): {R0}')
def _build_widget(self, node):
self.widget.children += (self._build_room(node.room),)
self.widget.children += (self._build_ventilation(node.ventilation),)
......
......@@ -132,7 +132,7 @@ class Mask:
def exhale_efficiency(self):
# Overall efficiency with the effect of the leaks for aerosol emission
# Gammaitoni et al (1997)
return self.η_exhale - (self.η_exhale * self.η_leaks)
return self.η_exhale * (1 - self.η_leaks)
Mask.types = {
......@@ -287,3 +287,26 @@ class Model:
# Concentration while infected not present.
end_concentration = self.concentration(t1)
return (end_concentration + ((np.exp(IVRR * t1) - 1) * end_concentration)) * np.exp(-IVRR * t)
def infection_probability(self):
# Infection probability
# Probability of COVID-19 Infection
exposure = 0.0 # q/m3*h
def integrate(fn, start, stop):
values = np.linspace(start, stop)
return np.trapz([fn(v) for v in values], values)
# TODO: Have this for exposed not infected.
for start, stop in self.infected.present_times:
exposure += (integrate(self.concentration, start, stop))
inf_aero = (
self.exposed_activity.inhalation_rate *
(1 - self.infected.mask.η_inhale) *
exposure
)
# Probability of infection.
return (1 - np.exp(-inf_aero)) * 100
import numpy as np
import numpy.testing as npt
import pytest
......@@ -59,7 +58,7 @@ def baseline_periodic_hepa():
return models.PeriodicHEPA(period=120, duration=15, q_air_mech=514.74)
def test_r0(baseline_model):
def test_concentrations(baseline_model):
ts = [0, 4, 5, 7, 10]
concentrations = [baseline_model.concentration(t) for t in ts]
npt.assert_allclose(
......@@ -69,6 +68,11 @@ def test_r0(baseline_model):
)
def test_r0(baseline_model):
p = baseline_model.infection_probability()
npt.assert_allclose(p, 93.305049)
def test_periodic_window(baseline_periodic_window, baseline_room):
# Interesting transition times for a period of 120 and duration of 15.
ts = [0, 14/60, 15/60, 16/60, 119/60, 120/60, 121/60, 239/60, 240/60]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment