diff --git a/caimira/calculator/models/models.py b/caimira/calculator/models/models.py index b0707a4a..b49f6692 100644 --- a/caimira/calculator/models/models.py +++ b/caimira/calculator/models/models.py @@ -687,13 +687,8 @@ def particle(self) -> Particle: def aerosols(self, mask: Mask): """ - Total volume of aerosols expired per volume of exhaled air (mL/cm^3). - """ - raise NotImplementedError("Subclass must implement") - - def jet_origin_concentration(self): - """ - Concentration of viruses at the jet origin (mL/m3). + Total volume of aerosols expired per volume of exhaled air + considering the outward mask efficiency (mL/cm^3). """ raise NotImplementedError("Subclass must implement") @@ -723,8 +718,9 @@ def particle(self) -> Particle: @cached() def aerosols(self, mask: Mask): """ - Total volume of aerosols expired per volume of exhaled air. - Result is in mL.cm^-3 + Total volume of aerosols expired per volume + of exhaled air considering the outward mask + efficiency. Result is in mL.cm^-3. """ def volume(d): return (np.pi * d**3) / 6. @@ -733,14 +729,6 @@ def volume(d): return self.cn * (volume(self.diameter) * (1 - mask.exhale_efficiency(self.diameter))) * 1e-12 - @cached() - def jet_origin_concentration(self): - def volume(d): - return (np.pi * d**3) / 6. - - # Final result converted from microns^3/cm3 to mL/m3 - return self.cn * volume(self.diameter) * 1e-6 - @dataclass(frozen=True) class MultipleExpiration(_ExpirationBase): @@ -889,9 +877,9 @@ def aerosols(self): 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, - 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. + The emission rate of infectious respiratory particles (IRP) in the expired air per + mL of respiratory fluid, if the infected population is present, in (virions.cm^3)/(mL.h). + This method returns only the diameter-independent variables within the emission rate. It should not be a function of time. """ raise NotImplementedError("Subclass must implement") @@ -900,12 +888,12 @@ def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat: def emission_rate_per_person_when_present(self) -> _VectorisedFloat: """ The emission rate if the infected population is present, per person - (in virions / h). + (in virions/h). """ return (self.emission_rate_per_aerosol_per_person_when_present() * self.aerosols()) - def emission_rate(self, time) -> _VectorisedFloat: + def emission_rate(self, time: float) -> _VectorisedFloat: """ The emission rate of the population vs time. """ @@ -945,13 +933,13 @@ def aerosols(self): @method_cache 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, - 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. + The emission rate of infectious respiratory particles (IRP) in the expired air per + mL of respiratory fluid, if the infected population is present, in (virions.cm^3)/(mL.h). + This method returns only the diameter-independent variables within the emission rate. It should not be a function of time. """ return self.known_individual_emission_rate - + @dataclass(frozen=True) class InfectedPopulation(_PopulationWithVirus): @@ -974,17 +962,19 @@ def aerosols(self): @method_cache 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). - This method includes only the diameter-independent variables within the emission rate. + The emission rate of infectious respiratory particles (IRP) in the expired air per + mL of respiratory fluid, if the infected population is present, in (virions.cm^3)/(mL.h). + This method returns only the diameter-independent variables within the emission rate. It should not be a function of time. """ - # Note on units: exhalation rate is in m^3/h -> 1e6 conversion factor - # Returns the emission rate times the number of infected hosts in the room + # Conversion factor explanation: + # The exhalation rate is in m^3/h, therefore the 1e6 conversion factor + # is to convert m^3/h into cm^3/h to return (virions.cm^3)/(mL.h), + # so that we can then multiply by aerosols (mL/cm^3). ER = (self.virus.viral_load_in_sputum * self.activity.exhalation_rate * self.fraction_of_infectious_virus() * - 10 ** 6) + 10**6) return ER @property @@ -1335,7 +1325,7 @@ def dilution_factor(self) -> _VectorisedFloat: ''' The dilution factor for the respective expiratory activity type. ''' - _dilution_factor = self.data_registry.short_range_model['dilution_factor'] + _dilution_factor = self.data_registry.short_range_model['dilution_factor'] # Average mouth opening diameter (m) mouth_diameter: float = _dilution_factor['mouth_diameter'] # type: ignore @@ -1386,29 +1376,38 @@ def dilution_factor(self) -> _VectorisedFloat: xstar[distances >= xstar])/𝛽r1/(xstar[distances >= xstar] + x0))**3 return factors + + def _normed_jet_origin_concentration(self) -> _VectorisedFloat: + """ + The initial jet concentration at the source origin (mouth/nose), normalized by + normalization_factor in the ShortRange class (corresponding to the diameter-independent + variables). Results in mL.cm^-3. + """ + # The short range origin concentration does not consider the mask contribution. + return self.expiration.aerosols(mask=Mask.types['No mask']) def _long_range_normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: """ - Virus long-range exposure concentration normalized by the - virus viral load, as function of time. + Virus long-range exposure concentration normalized by normalization_factor in the + ShortRange class, as function of time. Results in mL.cm^-3. """ - return (concentration_model.concentration(time) / - concentration_model.virus.viral_load_in_sputum) + return (concentration_model.concentration(time) / self.normalization_factor(concentration_model.infected)) def _normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: """ Virus short-range exposure concentration, as a function of time. If the given time falls within a short-range interval it returns the - short-range concentration normalized by the virus viral load. Otherwise - it returns 0. + short-range concentration normalized by normalization_factor in the + Short-range class. Otherwise it returns 0. Results in mL.cm^-3. """ start, stop = self.presence.boundaries()[0] # Verifies if the given time falls within a short-range interaction if start <= time <= stop: dilution = self.dilution_factor() - jet_origin_concentration = self.expiration.jet_origin_concentration() - # Long-range concentration normalized by the virus viral load + # Jet origin concentration normalized by the emission rate (except the BR) + normed_jet_origin_concentration = self._normed_jet_origin_concentration() + # Long-range concentration normalized by the emission rate (except the BR) long_range_normed_concentration = self._long_range_normed_concentration(concentration_model, time) # The long-range concentration values are then approximated using interpolation: @@ -1420,15 +1419,32 @@ def _normed_concentration(self, concentration_model: ConcentrationModel, time: f # Short-range concentration formula. The long-range concentration is added in the concentration method (ExposureModel). # based on continuum model proposed by Jia et al (2022) - https://doi.org/10.1016/j.buildenv.2022.109166 - return ((1/dilution)*(jet_origin_concentration - long_range_normed_concentration_interpolated)) + return ((1/dilution)*(normed_jet_origin_concentration - long_range_normed_concentration_interpolated)) return 0. - + + def normalization_factor(self, infected: InfectedPopulation) -> _VectorisedFloat: + """ + The normalization factor applied to the short-range results. It refers to the emission + rate per aerosol without accounting for the exhalation rate (viral load and f_inf). + Result in (virions.cm^3)/(mL.m^3). + """ + # Re-use the emission rate method divided by the BR contribution. + return infected.emission_rate_per_aerosol_per_person_when_present() / infected.activity.exhalation_rate + + def jet_origin_concentration(self, infected: InfectedPopulation) -> _VectorisedFloat: + """ + The initial jet concentration at the source origin (mouth/nose). + Returns the full result with the diameter dependent and independent variables, in virions/m^3. + """ + return self._normed_jet_origin_concentration() * self.normalization_factor(infected) + def short_range_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: """ Virus short-range exposure concentration, as a function of time. + Factor of normalization applied back here. Results in virions/m^3. """ - return (self._normed_concentration(concentration_model, time) * - concentration_model.virus.viral_load_in_sputum) + return (self._normed_concentration(concentration_model, time) * + self.normalization_factor(concentration_model.infected)) @method_cache def _normed_short_range_concentration_cached(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: @@ -1462,16 +1478,16 @@ def extract_between_bounds(self, time1: float, time2: float) -> typing.Union[Non return start, stop def _normed_jet_exposure_between_bounds(self, - concentration_model: ConcentrationModel, time1: float, time2: float): """ Get the part of the integrated short-range concentration of viruses in the air, between the times start and stop, coming - from the jet concentration, normalized by the viral load, and - without dilution. + from the jet concentration, normalized by normalization_factor, + and without dilution. """ start, stop = self.extract_between_bounds(time1, time2) - jet_origin = self.expiration.jet_origin_concentration() + # Note the conversion factor mL.cm^-3 -> mL.m^-3 + jet_origin = self._normed_jet_origin_concentration() * 10**6 return jet_origin * (stop - start) def _normed_interpolated_longrange_exposure_between_bounds( @@ -1479,20 +1495,24 @@ def _normed_interpolated_longrange_exposure_between_bounds( time1: float, time2: float): """ Get the part of the integrated short-range concentration due - to the background concentration, normalized by the viral load - and the breathing rate, and without dilution. + to the background concentration, normalized by normalization_factor + together with breathing rate, and without dilution. One needs to interpolate the integrated long-range concentration for the particle diameters defined here. - TODO: make sure any potential extrapolation has a - negligible effect. """ start, stop = self.extract_between_bounds(time1, time2) if stop<=start: return 0. + # Note that for the correct integration one needs to isolate those parameters + # that are diameter-dependent from those that are diameter independent. + # Therefore, the diameter-independent parameters (viral load, f_ind and BR) + # are removed for the interpolation, and added back once the integration over + # the new aerosol diameters (done with the mean) is completed. normed_int_concentration = ( concentration_model.integrated_concentration(start, stop) /concentration_model.virus.viral_load_in_sputum + /concentration_model.infected.fraction_of_infectious_virus() /concentration_model.infected.activity.exhalation_rate ) normed_int_concentration_interpolated = np.interp( @@ -1566,7 +1586,8 @@ def fun(x): @dataclass(frozen=True) class ExposureModel: """ - Represents the exposure to a concentration of virions in the air. + Represents the exposure to a concentration of + infectious respiratory particles (IRP) in the air. """ data_registry: DataRegistry @@ -1693,7 +1714,6 @@ def long_range_deposited_exposure_between_bounds(self, time1: float, time2: floa 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. return deposited_exposure def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: @@ -1709,8 +1729,7 @@ def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _Vect deposited_exposure: _VectorisedFloat = 0. for interaction in self.short_range: start, stop = interaction.extract_between_bounds(time1, time2) - short_range_jet_exposure = interaction._normed_jet_exposure_between_bounds( - self.concentration_model, start, stop) + short_range_jet_exposure = interaction._normed_jet_exposure_between_bounds(start, stop) short_range_lr_exposure = interaction._normed_interpolated_longrange_exposure_between_bounds( self.concentration_model, start, stop) dilution = interaction.dilution_factor() @@ -1741,11 +1760,12 @@ def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _Vect interaction.activity.inhalation_rate /dilution) - # Then we multiply by diameter-independent quantities: viral load - # and fraction of infected virions + # Then we multiply by the emission rate without the BR contribution (and conversion factor), + # and parameters of the vD equation (i.e. n_in). deposited_exposure *= ( - self.concentration_model.virus.viral_load_in_sputum - * (1 - self.exposed.mask.inhale_efficiency())) + (self.concentration_model.infected.emission_rate_per_aerosol_per_person_when_present() / ( + self.concentration_model.infected.activity.exhalation_rate * 10**6)) * + (1 - self.exposed.mask.inhale_efficiency())) # Long-range concentration deposited_exposure += self.long_range_deposited_exposure_between_bounds(time1, time2) @@ -1823,7 +1843,9 @@ def expected_new_cases(self) -> _VectorisedFloat: The expected_new_cases may provide one or two different outputs: 1) Long-range exposure: take the infection_probability and multiply by the occupants exposed to long-range. 2) Short- and long-range exposure: take the infection_probability of long-range multiplied by the occupants exposed to long-range only, - plus the infection_probability of short- and long-range multiplied by the occupants exposed to short-range only. """ + plus the infection_probability of short- and long-range multiplied by the occupants exposed to short-range only. + """ + if self.short_range != (): new_cases_long_range = nested_replace(self, {'short_range': (),}).infection_probability() * (self.exposed.number - self.exposed_to_short_range) return (new_cases_long_range + (self.infection_probability() * self.exposed_to_short_range)) / 100 diff --git a/caimira/calculator/report/report_generator.py b/caimira/calculator/report/report_generator.py index 9ed8316a..7c6f5cea 100644 --- a/caimira/calculator/report/report_generator.py +++ b/caimira/calculator/report/report_generator.py @@ -256,11 +256,14 @@ def uncertainties_plot(infection_probability: models._VectorisedFloat, lower_percentiles: models._VectorisedFloat, upper_percentiles: models._VectorisedFloat): - fig, axs = plt.subplots(2, 3, + fig, axes = plt.subplots(2, 3, gridspec_kw={'width_ratios': [5, 0.5] + [1], 'height_ratios': [3, 1], 'wspace': 0}, sharey='row', sharex='col') + + # Type hint for axs + axs: np.ndarray = np.array(axes) for y, x in [(0, 1)] + [(1, i + 1) for i in range(2)]: axs[y, x].axis('off') diff --git a/caimira/calculator/tests/models/test_short_range_model.py b/caimira/calculator/tests/models/test_short_range_model.py index f762a777..0382369a 100644 --- a/caimira/calculator/tests/models/test_short_range_model.py +++ b/caimira/calculator/tests/models/test_short_range_model.py @@ -49,7 +49,7 @@ def test_short_range_model_ndarray(concentration_model, short_range_model): model = short_range_model.build_model(SAMPLE_SIZE) assert isinstance(model._normed_concentration(concentration_model, 10.75), np.ndarray) assert isinstance(model.short_range_concentration(concentration_model, 10.75), np.ndarray) - assert isinstance(model._normed_jet_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) + assert isinstance(model._normed_jet_exposure_between_bounds(10.75, 10.85), np.ndarray) assert isinstance(model._normed_interpolated_longrange_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) assert isinstance(model.short_range_concentration(concentration_model, 14.0), float) @@ -106,9 +106,9 @@ def test_extract_between_bounds(short_range_model, time1, time2, @pytest.mark.parametrize( "time, expected_short_range_concentration", [ [8.5, 0.], - [10.5, 11.266605], - [10.6, 11.266605], - [11.0, 11.266605], + [10.5, 5.6333025], + [10.6, 5.6333025], + [11.0, 5.6333025], [12.0, 0.], ] ) diff --git a/caimira/calculator/tests/test_full_algorithm.py b/caimira/calculator/tests/test_full_algorithm.py index fb5f1c78..adb93df0 100644 --- a/caimira/calculator/tests/test_full_algorithm.py +++ b/caimira/calculator/tests/test_full_algorithm.py @@ -16,6 +16,7 @@ expiration_BLO_factors,short_range_expiration_distributions, short_range_distances,virus_distributions,activity_distributions) + SAMPLE_SIZE = 1_000_000 TOLERANCE = 0.04 @@ -55,7 +56,7 @@ class SimpleConcentrationModel: #: Number of infected people num_infected: int = 1 - #: Fraction of infected viruses (viable to RNA ratio) + #: Viable to RNA ratio viable_to_RNA: _VectorisedFloat = 0.5 #: Host immunity factor (0. for not immune) @@ -97,6 +98,12 @@ def removal_rate(self) -> _VectorisedFloat: return (self.lambda_ventilation + ln2/(np.where(hl_calc <= 0, 6.43, np.minimum(6.43, hl_calc)))) + def fraction_of_infectious_virus(self) -> _VectorisedFloat: + """ + The fraction of infectious virus. + """ + return self.viable_to_RNA * (1 - self.HI) + @method_cache def deposition_removal_coefficient(self) -> float: """ @@ -181,8 +188,8 @@ def concentration(self,t: float) -> _VectorisedFloat: return ( ( (0 if not self.infected_presence.triggered(t) else self.f(lambda_rate,0)) + result * np.exp(-lambda_rate*(t-ti)) ) - * self.num_infected * self.viable_to_RNA - * (1. - self.HI) / self.room_volume) + * self.num_infected * self.fraction_of_infectious_virus() + / self.room_volume) @dataclass(frozen=True) @@ -263,6 +270,7 @@ def jet_concentration(self,conc_model: SimpleConcentrationModel) -> _VectorisedF we perform the integral of Np(d)*V(d) over diameter analytically """ vl = conc_model.viral_load + viable_to_RNA = conc_model.viable_to_RNA dmin = self.diameter_min dmax = self.diameter_max result = 0. @@ -273,7 +281,7 @@ def jet_concentration(self,conc_model: SimpleConcentrationModel) -> _VectorisedF ymax = (np.log(dmax)-mu)/(sqrt2*sigma)-3.*sigma/sqrt2 result += ( (cn * famp * d0**3)/2. * np.exp(9*sigma**2/2.) * (erf(ymax) - erf(ymin)) ) - return vl * 1e-6 * result * np.pi/6. + return vl * viable_to_RNA * 1e-6 * result * np.pi/6. def concentration(self, conc_model: SimpleConcentrationModel, time: float) -> _VectorisedFloat: """ @@ -410,8 +418,8 @@ def primitive(time): else self.f_with_fdep(lambda_rate,0,evaporation)*(t2-t1)) + (primitive(t2) * np.exp(-lambda_rate*(t2-ti)) - primitive(t1) * np.exp(-lambda_rate*(t1-ti)) ) ) - * self.num_infected * self.viable_to_RNA - * (1. - self.HI) / self.room_volume) + * self.num_infected * self.fraction_of_infectious_virus() + / self.room_volume) @method_cache def integrated_shortrange_concentration(self) -> _VectorisedFloat: @@ -430,7 +438,7 @@ def integrated_shortrange_concentration(self) -> _VectorisedFloat: res = (quad(integrand, sr_model.diameter_min,sr_model.diameter_max, epsabs=0.,limit=500)[0] - * self.viral_load * 1e-6 * (t2-t1) ) + * self.viral_load * self.fraction_of_infectious_virus() * 1e-6 * (t2-t1) ) result += sr_model.breathing_rate * ( res-self.integrated_longrange_concentration(t1,t2,evaporation) )/sr_model.dilution_factor() @@ -842,3 +850,119 @@ def test_exposure_with_shortrange_and_distributions(expo_sr_model_distr, rtol=0.03 ) + +def exposure_model_from_parameter(data_registry, short_range_models, f_inf=0.5, viral_load=1e9, BR=1.25): + virus: models.SARSCoV2 = models.SARSCoV2( + viral_load_in_sputum=viral_load, + infectious_dose=50, + viable_to_RNA_ratio=f_inf, + transmissibility_factor=0.51, + ) + c_model = mc.ConcentrationModel( + data_registry=data_registry, + room=models.Room(volume=50, humidity=0.3), + ventilation=models.AirChange(active=models.PeriodicInterval(period=120, duration=120), + air_exch=10_000_000), + infected=mc.InfectedPopulation( + data_registry=data_registry, + number=1, + presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + virus=virus, + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + expiration=expiration_distributions(data_registry)['Breathing'], + host_immunity=0., + ), + evaporation_factor=0.3, + ) + return mc.ExposureModel( + data_registry=data_registry, + concentration_model=c_model, + short_range=short_range_models, + exposed=mc.Population( + number=1, + presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + mask=models.Mask.types['No mask'], + activity=models.Activity(inhalation_rate=BR, exhalation_rate=1.25), + host_immunity=0., + ), + geographical_data=models.Cases(), + ).build_model(SAMPLE_SIZE) + + +@retry(tries=10) +def test_exposure_scale_with_f_inf(data_registry, sr_models): + """ + Exposure scaling test for the fraction of infectious virus. + """ + e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, f_inf=0.5) + e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, f_inf=1) + np.testing.assert_allclose( + 2*e_model_1.deposited_exposure().mean(), + e_model_2.deposited_exposure().mean(), rtol=0.02 + ) + + +@retry(tries=10) +def test_exposure_scale_with_viral_load(data_registry, sr_models): + """ + Exposure scaling test for the viral load. + """ + e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, viral_load=1e9) + e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, viral_load=2e9) + np.testing.assert_allclose( + 2*e_model_1.deposited_exposure().mean(), + e_model_2.deposited_exposure().mean(), rtol=0.02 + ) + + +@retry(tries=10) +def test_lr_exposure_scale_with_breathing_rate(data_registry): + """ + Exposure scaling test for the breathing rate when there are only long-range + interactions defined. Only the inhalation rate of the infected takes place + at the deposited exposure level. + """ + e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=(), BR=1.25) + e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=(), BR=2.5) + np.testing.assert_allclose( + 2*e_model_1.deposited_exposure().mean(), + e_model_2.deposited_exposure().mean(), rtol=0.02 + ) + + +@retry(tries=10) +def test_exposure_scale_with_breathing_rate(data_registry, sr_models): + """ + Exposure scaling test for the breathing rate when long- and short-range + interactions are defined. We need to apply the multiplication factor + to the inhalation rate of the infected (long-range), but also for + each short-range interaction. + """ + e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, BR=1.25) + + seated_act = models.Activity.types['Seated'] + heavy_exercise_act = models.Activity.types['Heavy exercise'] + sr_models_activity = ( + mc.ShortRangeModel( + data_registry = data_registry, + expiration = short_range_expiration_distributions(data_registry)['Speaking'], + activity = models.Activity(inhalation_rate=seated_act.inhalation_rate * 2, + exhalation_rate=seated_act.exhalation_rate), + presence = interaction_intervals[0], + distance = 0.854, + ), + mc.ShortRangeModel( + data_registry = data_registry, + expiration = short_range_expiration_distributions(data_registry)['Breathing'], + activity = models.Activity(inhalation_rate=heavy_exercise_act.inhalation_rate * 2, + exhalation_rate=heavy_exercise_act.inhalation_rate), + presence = interaction_intervals[1], + distance = 0.854, + ), + ) + e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models_activity, BR=2.5) + np.testing.assert_allclose( + 2*e_model_1.deposited_exposure().mean(), + e_model_2.deposited_exposure().mean(), rtol=0.02 + )