From 739864cc28407cab84b1a5abd10a7064887004fb Mon Sep 17 00:00:00 2001
From: Phil Elson <philip.elson@cern.ch>
Date: Mon, 26 Oct 2020 19:10:53 +0100
Subject: [PATCH] Add infection probability calculation to the model, and
 introduce it into the model output in a rudimentary form.

---
 cara/apps.py                        | 13 ++++++++++++-
 cara/models.py                      | 26 +++++++++++++++++++++++++-
 cara/tests/test_known_quantities.py |  8 ++++++--
 3 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/cara/apps.py b/cara/apps.py
index 96edf037..09fd55f4 100644
--- a/cara/apps.py
+++ b/cara/apps.py
@@ -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),)
diff --git a/cara/models.py b/cara/models.py
index 69f642a6..998d025c 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -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 = {
@@ -294,3 +294,27 @@ class Model:
             end_concentration = self.concentration(t1)
             fac = np.exp(IVRR * (t1 - t))
             return end_concentration * fac
+
+    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
+
diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py
index b87a318b..9dc36f38 100644
--- a/cara/tests/test_known_quantities.py
+++ b/cara/tests/test_known_quantities.py
@@ -1,4 +1,3 @@
-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.196908)
+
+
 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]