add MCExposureModel class

cara/montecarlo.py

@@ -175,6 +175,63 @@ class MCConcentrationModel:
         return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac
 
 
+@dataclass(frozen=True)
+class MCExposureModel:
+    #: The virus concentration model which this exposure model should consider.
+    concentration_model: MCConcentrationModel
+
+    #: The population of non-infected people to be used in the model.
+    exposed: models.Population
+
+    #: The number of times the exposure event is repeated (default 1).
+    repeats: int = 1
+
+    def quanta_exposure(self) -> np.ndarray:
+        """The number of virus quanta per meter^3."""
+        exposure = np.zeros(self.concentration_model.infected.samples)
+
+        for start, stop in self.exposed.presence.boundaries():
+            concentrations = np.asarray([self.concentration_model.concentration(t) for t in np.linspace(start, stop)])
+            integrals = np.trapz(concentrations, axis=0)
+            exposure += integrals
+
+        return exposure * self.repeats
+
+    def infection_probability(self):
+        exposure = self.quanta_exposure()
+
+        inf_aero = (
+            self.exposed.activity.inhalation_rate *
+            (1 - self.exposed.mask.η_inhale) *
+            exposure
+        )
+
+        # Probability of infection.
+        return (1 - np.exp(-inf_aero)) * 100
+
+    def expected_new_cases(self):
+        prob = self.infection_probability()
+        exposed_occupants = self.exposed.number
+        return prob * exposed_occupants / 100
+
+    # def reproduction_number(self):
+    #     """
+    #     The reproduction number can be thought of as the expected number of
+    #     cases directly generated by one infected case in a population.
+    #
+    #     """
+    #     if self.concentration_model.infected.number == 1:
+    #         return self.expected_new_cases()
+    #
+    #     # Create an equivalent exposure model but with precisely
+    #     # one infected case.
+    #     single_exposure_model = nested_replace(
+    #         self, {'concentration_model.infected.number': 1}
+    #     )
+    #
+    #     return single_exposure_model.expected_new_cases()
+
+
 def logscale_hist(x: typing.Iterable, bins: int) -> None:
     """
     Plots the data of x as a log-scale histogram