diff --git a/cara/montecarlo.py b/cara/montecarlo.py index 1302ea06..34132c83 100644 --- a/cara/montecarlo.py +++ b/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