Scenarios reformulation and correct Cn values
This commit is contained in:
Luis Aleixo
parent
f4587988db
commit
ca1db4ef59
3 changed files with 299 additions and 846 deletions
|
|
@ -1,7 +1,11 @@
|
|||
from cara import models
|
||||
from cara.monte_carlo.data import activity_distributions, symptomatic_vl_frequencies, viable_to_RNA_ratio_distribution, infectious_dose_distribution
|
||||
from cara.monte_carlo.data import activity_distributions, symptomatic_vl_frequencies, viable_to_RNA_ratio_distribution, infectious_dose_distribution, expiration_distributions
|
||||
import cara.monte_carlo as mc
|
||||
import numpy as np
|
||||
from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel,CustomKernel,Uniform
|
||||
from cara.monte_carlo.data import BLOmodel
|
||||
import typing
|
||||
from cara.apps.calculator.model_generator import build_expiration
|
||||
|
||||
######### Scatter points (data taken: copies per hour) #########
|
||||
|
||||
|
|
@ -26,10 +30,30 @@ milton_er = [22, 220, 1120]
|
|||
yann_vl = [np.log10(7.86E+07), np.log10(2.23E+09), np.log10(1.51E+10)]
|
||||
yann_er = [8396.78166, 45324.55964, 400054.0827]
|
||||
|
||||
def cn_expiration_distribution(BLO_factors, cn_values):
|
||||
"""
|
||||
Returns an Expiration with an aerosol diameter distribution, defined
|
||||
by the BLO factors (a length-3 tuple).
|
||||
The total concentration of aerosols is computed by integrating
|
||||
the distribution between 0.1 and 30 microns - these boundaries are
|
||||
an historical choice based on previous implementations of the model
|
||||
(it limits the influence of the O-mode).
|
||||
"""
|
||||
dscan = np.linspace(0.1, 30. ,3000)
|
||||
return mc.Expiration(CustomKernel(dscan,
|
||||
BLOmodel(BLO_factors, cn_values).distribution(dscan),kernel_bandwidth=0.1),
|
||||
BLOmodel(BLO_factors, cn_values).integrate(0.1, 30.))
|
||||
|
||||
expiration_BLO_factors = {
|
||||
'Breathing': (1., 0., 0.),
|
||||
'Talking': (1., 1., 1.),
|
||||
'Singing': (1., 5., 5.),
|
||||
'Shouting': (1., 5., 5.),
|
||||
}
|
||||
|
||||
######### Standard exposure models ###########
|
||||
|
||||
######### Breathing model ###########
|
||||
def breathing_exposure(activity: str, mask: str):
|
||||
def exposure_module(activity: str, expiration: str, mask: str):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
|
|
@ -48,90 +72,22 @@ def breathing_exposure(activity: str, mask: str):
|
|||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
expiration=expiration_distributions[expiration],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types[activity],
|
||||
activity=activity_distributions[activity],
|
||||
mask=models.Mask.types[mask],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Speaking model ###########
|
||||
def speaking_exposure(activity: str, mask: str):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types[activity],
|
||||
mask=models.Mask.types[mask],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Shouting model ###########
|
||||
def shouting_exposure(activity: str, mask: str):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=models.Expiration.types['Shouting'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types[activity],
|
||||
mask=models.Mask.types[mask],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Breathing model for specific viral load ###########
|
||||
def breathing_exposure_vl(vl):
|
||||
######### Exposure model for specific viral load ###########
|
||||
def exposure_vl(activity: str, expiration: str, mask: str, vl: float):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
|
|
@ -148,24 +104,24 @@ def breathing_exposure_vl(vl):
|
|||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=expiration_distributions[expiration],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions[activity],
|
||||
mask=models.Mask.types[mask],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Talking model for specific viral load ###########
|
||||
def talking_exposure_vl(vl):
|
||||
######### Exposure model for specific viral load ###########
|
||||
def exposure_vl_cn(activity: str, expiration: str, mask: str, vl: float, cn: typing.Tuple[float, float, float]):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
|
|
@ -182,358 +138,23 @@ def talking_exposure_vl(vl):
|
|||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=cn_expiration_distribution(expiration_BLO_factors[expiration], cn),
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions[activity],
|
||||
mask=models.Mask.types[mask],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Shouting model for specific viral load ###########
|
||||
def shouting_exposure_vl(vl):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=10**vl,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types['No mask'],
|
||||
activity=activity_distributions['Light activity'],
|
||||
expiration=models.Expiration.types['Shouting'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Light activity'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Used for CDF Models ###########
|
||||
######### Breathing Models #########
|
||||
def breathing_seated_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
def breathing_light_activity_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Light activity'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Light activity'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
def breathing_heavy_exercise_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Heavy exercise'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Heavy exercise'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Speaking Models #########
|
||||
def speaking_seated_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
def speaking_light_activity_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Light activity'],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Light activity'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
def speaking_heavy_exercise_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Heavy exercise'],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Heavy exercise'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Shouting Models #########
|
||||
def shouting_seated_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.Expiration.types['Shouting'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
def shouting_light_activity_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Light activity'],
|
||||
expiration=models.Expiration.types['Shouting'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Light activity'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
def shouting_heavy_exercise_exposure():
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Heavy exercise'],
|
||||
expiration=models.Expiration.types['Shouting'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Heavy exercise'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
########## Concentration curves ###########
|
||||
########## Concentration curves for specific scenarios ###########
|
||||
def office_model_no_mask_windows_closed():
|
||||
office_model_no_vent = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
|
|
@ -552,18 +173,14 @@ def office_model_no_mask_windows_closed():
|
|||
),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.MultipleExpiration(
|
||||
expirations = (models.Expiration.types['Talking'],
|
||||
models.Expiration.types['Breathing']),
|
||||
weights=(1, 2)
|
||||
),
|
||||
expiration=build_expiration({'Talking': 0.33, 'Breathing': 0.67}),
|
||||
host_immunity=0.,
|
||||
)
|
||||
),
|
||||
exposed=models.Population(
|
||||
number=18,
|
||||
presence=models.SpecificInterval(present_times = ((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
|
||||
activity=models.Activity.types['Seated'],
|
||||
activity=activity_distributions['Seated'],
|
||||
mask=models.Mask.types['No mask'],
|
||||
host_immunity=0.,
|
||||
)
|
||||
|
|
@ -597,18 +214,14 @@ def office_model_no_mask_windows_open_breaks():
|
|||
),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.MultipleExpiration(
|
||||
expirations = (models.Expiration.types['Talking'],
|
||||
models.Expiration.types['Breathing']),
|
||||
weights=(1, 2)
|
||||
),
|
||||
expiration=build_expiration({'Talking': 0.33, 'Breathing': 0.67}),
|
||||
host_immunity=0.,
|
||||
)
|
||||
),
|
||||
exposed=models.Population(
|
||||
number=18,
|
||||
presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
|
||||
activity=models.Activity.types['Seated'],
|
||||
activity=activity_distributions['Seated'],
|
||||
mask=models.Mask.types['No mask'],
|
||||
host_immunity=0.,
|
||||
)
|
||||
|
|
@ -641,108 +254,16 @@ def office_model_no_mask_windows_open_alltimes():
|
|||
),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.MultipleExpiration(
|
||||
expirations = (models.Expiration.types['Talking'],
|
||||
models.Expiration.types['Breathing']),
|
||||
weights=(1, 2)
|
||||
),
|
||||
expiration=build_expiration({'Talking': 0.33, 'Breathing': 0.67}),
|
||||
host_immunity=0.,
|
||||
)
|
||||
),
|
||||
exposed=models.Population(
|
||||
number=18,
|
||||
presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
|
||||
activity=models.Activity.types['Seated'],
|
||||
activity=activity_distributions['Seated'],
|
||||
mask=models.Mask.types['No mask'],
|
||||
host_immunity=0.,
|
||||
)
|
||||
)
|
||||
return office_model_no_vent
|
||||
|
||||
|
||||
######### Standard exposure models ###########
|
||||
|
||||
######### Breathing model ###########
|
||||
def breathing_exposure(activity: str, mask: str):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types[activity],
|
||||
mask=models.Mask.types[mask],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Speaking model ###########
|
||||
def speaking_exposure(activity: str, mask: str):
|
||||
exposure_mc = mc.ExposureModel(
|
||||
concentration_model=mc.ConcentrationModel(
|
||||
room=models.Room(volume=100, humidity=0.5),
|
||||
ventilation=models.AirChange(
|
||||
active=models.SpecificInterval(((0, 24),)),
|
||||
air_exch=0.25,
|
||||
),
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=models.Expiration.types['Talking'],
|
||||
host_immunity=0.,
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types[activity],
|
||||
mask=models.Mask.types[mask],
|
||||
host_immunity=0.,
|
||||
),
|
||||
)
|
||||
return exposure_mc
|
||||
|
||||
######### Infected Population model ###########
|
||||
def infected_model(mask: str, activity: str, expiratory_activity: str):
|
||||
infected=mc.InfectedPopulation(
|
||||
number=1,
|
||||
virus=mc.Virus(
|
||||
viral_load_in_sputum=symptomatic_vl_frequencies,
|
||||
infectious_dose=infectious_dose_distribution,
|
||||
viable_to_RNA_ratio=viable_to_RNA_ratio_distribution,
|
||||
transmissibility_factor=1.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types[mask],
|
||||
activity=activity_distributions[activity],
|
||||
expiration=models.Expiration.types[expiratory_activity])
|
||||
|
||||
return infected
|
||||
return office_model_no_vent
|
||||
|
|
@ -13,27 +13,28 @@ from dataclasses import dataclass
|
|||
|
||||
# Exhaled virions while talking, seated #
|
||||
#print('\n<<<<<<<<<<< Vlout for Talking, seated >>>>>>>>>>>')
|
||||
#exposure_model_from_vl_talking()
|
||||
#exposure_model_from_vl(activity='Seated', expiration='Talking', mask='No mask')
|
||||
|
||||
# Exhaled virions while breathing, seated #
|
||||
#print('\n<<<<<<<<<<< Vlout for Breathing, seated >>>>>>>>>>>')
|
||||
#exposure_model_from_vl_breathing()
|
||||
#exposure_model_from_vl(activity='Seated', expiration='Breathing', mask='No mask')
|
||||
|
||||
# Exhaled virions while breathing, light activity #
|
||||
#print('\n<<<<<<<<<<< Vlout for Shouting, light activity >>>>>>>>>>>')
|
||||
#exposure_model_from_vl_shouting()
|
||||
#exposure_model_from_vl(activity='Light activity', expiration='Shouting', mask='No mask')
|
||||
|
||||
# Exhaled virions while talking according to BLO model, seated #
|
||||
#print('\n<<<<<<<<<<< Vlout for Talking, seated with chosen Cn,L >>>>>>>>>>>')
|
||||
#exposure_model_from_vl_talking_cn()
|
||||
#exposure_model_from_vl_cn(activity='Seated', expiration='Talking', mask='No mask')
|
||||
#print('\n')
|
||||
|
||||
# Exhaled virions while breathing according to BLO model, seated #
|
||||
#print('\n<<<<<<<<<<< Vlout for Breathing, seated with chosen Cn,B >>>>>>>>>>>')
|
||||
#exposure_model_from_vl_breathing_cn()
|
||||
#exposure_model_from_vl_cn(activity='Seated', expiration='Breathing', mask='No mask')
|
||||
#print('\n')
|
||||
|
||||
############ Plots with viral loads and emission rates + statistical data ############
|
||||
present_vl_er_histograms(activity='Seated', mask='No mask')
|
||||
#present_vl_er_histograms(activity='Seated', mask='No mask')
|
||||
#present_vl_er_histograms(activity='Light activity', mask='No mask')
|
||||
#present_vl_er_histograms(activity='Heavy exercise', mask='No mask')
|
||||
|
||||
|
|
@ -48,7 +49,7 @@ present_vl_er_histograms(activity='Seated', mask='No mask')
|
|||
#compare_concentration_curves()
|
||||
|
||||
############ Emission Rate Violin plot ############
|
||||
#compare_viruses_vr()
|
||||
compare_viruses_vr()
|
||||
|
||||
############ Used for testing ############
|
||||
#exposure_model_from_vl_talking_new_points()
|
||||
|
|
|
|||
|
|
@ -1,6 +1,3 @@
|
|||
from numpy.core.function_base import linspace
|
||||
from cara import models
|
||||
from cara.monte_carlo.data import activity_distributions
|
||||
from tqdm import tqdm
|
||||
from matplotlib.patches import Rectangle, Patch
|
||||
from scipy.spatial import ConvexHull
|
||||
|
|
@ -12,22 +9,50 @@ import pandas as pd
|
|||
import matplotlib.lines as mlines
|
||||
import matplotlib.patches as patches
|
||||
import matplotlib as mpl
|
||||
from random import randrange
|
||||
|
||||
|
||||
######### Plot material #########
|
||||
np.random.seed(2000)
|
||||
SAMPLE_SIZE = 250000
|
||||
viral_loads = np.linspace(2, 12, 600)
|
||||
|
||||
############# Markers (for legend) #############
|
||||
markers = [5, 'd', 4]
|
||||
|
||||
""" Exhaled virions from exposure models """
|
||||
|
||||
######### Breathing #########
|
||||
def emission_rate_when_present(exposure_model: mc.ExposureModel):
|
||||
aerosols = exposure_model.concentration_model.infected.expiration.aerosols(
|
||||
exposure_model.concentration_model.infected.mask).mean()
|
||||
exhalation_rate = exposure_model.concentration_model.infected.activity.exhalation_rate
|
||||
viral_load_in_sputum = exposure_model.concentration_model.infected.virus.viral_load_in_sputum
|
||||
return (viral_load_in_sputum * exhalation_rate * 10 ** 6 * aerosols) * exposure_model.concentration_model.infected.number
|
||||
|
||||
|
||||
def exposure_model_from_vl_breathing():
|
||||
def _normed_exposure_between_bounds(model: mc.ExposureModel, time1: float, time2: float):
|
||||
"""The number of virions per meter^3 between any two times, normalized
|
||||
by the emission rate of the infected population"""
|
||||
for start, stop in model.exposed.presence.boundaries():
|
||||
if start > time2:
|
||||
normed_exposure = 0.
|
||||
break
|
||||
elif time2 <= stop:
|
||||
normed_exposure = model.concentration_model.normed_integrated_concentration(
|
||||
time1, time2)
|
||||
break
|
||||
else:
|
||||
normed_exposure = model.concentration_model.normed_integrated_concentration(
|
||||
time1, time2)
|
||||
return normed_exposure
|
||||
|
||||
|
||||
def exposure_between_bounds(model: mc.ExposureModel, time1: float, time2: float):
|
||||
"""The number of virions per meter^3 between any two times."""
|
||||
return (_normed_exposure_between_bounds(model, time1, time2) *
|
||||
emission_rate_when_present(model))
|
||||
|
||||
######### Exhaled virions from exposure models #########
|
||||
|
||||
|
||||
def exposure_model_from_vl(activity, expiration, mask):
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
|
||||
|
|
@ -38,50 +63,80 @@ def exposure_model_from_vl_breathing():
|
|||
upper_percentiles = []
|
||||
|
||||
for vl in tqdm(viral_loads):
|
||||
exposure_mc = breathing_exposure_vl(vl)
|
||||
exposure_mc = exposure_vl(activity, expiration, mask, vl)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
# divide by 2 to have in 30min (half an hour)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present() / 2
|
||||
|
||||
if expiration == 'Breathing':
|
||||
# divide by 2 to have in 30min (half an hour)
|
||||
emission_rate = emission_rate_when_present(exposure_model) / 2
|
||||
elif expiration == 'Talking':
|
||||
# divide by 4 to have in 15min (quarter of an hour)
|
||||
emission_rate = emission_rate_when_present(exposure_model) / 4
|
||||
elif expiration == 'Shouting':
|
||||
emission_rate = emission_rate_when_present(exposure_model)
|
||||
|
||||
er_means.append(np.mean(emission_rate))
|
||||
er_medians.append(np.median(emission_rate))
|
||||
lower_percentiles.append(np.quantile(emission_rate, 0.01))
|
||||
upper_percentiles.append(np.quantile(emission_rate, 0.99))
|
||||
emission_rate_1h = exposure_model.concentration_model.infected.emission_rate_when_present()
|
||||
emission_rate_1h = emission_rate_when_present(exposure_model)
|
||||
er_means_1h.append(np.mean(emission_rate_1h))
|
||||
|
||||
# divide by 2 to have in 30min (half an hour)
|
||||
coleman_etal_er_breathing_2 = [x/2 for x in coleman_etal_er_breathing]
|
||||
milton_er_2 = [x/2 for x in milton_er]
|
||||
yann_er_2 = [x/2 for x in yann_er]
|
||||
if expiration == 'Breathing':
|
||||
# divide by 2 to have in 30min (half an hour)
|
||||
coleman_etal_er_breathing_2 = [x/2 for x in coleman_etal_er_breathing]
|
||||
milton_er_2 = [x/2 for x in milton_er]
|
||||
yann_er_2 = [x/2 for x in yann_er]
|
||||
|
||||
ratio = np.mean(10**viral_loads / er_means)
|
||||
ratio_1h = np.mean(10**viral_loads / er_means_1h)
|
||||
print('Mean swab-to-aersol vl ratio in 30min:')
|
||||
print(format(ratio, "5.1e"))
|
||||
print('Mean swab-to-aersol vl ratio emission rate per hour:')
|
||||
print(format(ratio_1h, "5.1e"))
|
||||
|
||||
############# Coleman #############
|
||||
scatter_coleman_data(coleman_etal_vl_breathing,
|
||||
coleman_etal_er_breathing_2)
|
||||
|
||||
############# Milton et al #############
|
||||
scatter_milton_data(milton_vl, milton_er_2)
|
||||
|
||||
############# Yan et al #############
|
||||
scatter_yann_data(yann_vl, yann_er_2)
|
||||
|
||||
############ Legend ############
|
||||
build_breathing_legend(fig)
|
||||
|
||||
elif expiration == 'Talking':
|
||||
# divide by 4 to have in 15min (quarter of an hour)
|
||||
coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
|
||||
|
||||
ratio = np.mean(10**viral_loads / er_means)
|
||||
ratio_1h = np.mean(10**viral_loads / er_means_1h)
|
||||
print('Mean swab-to-aersol vl ratio in 30min:')
|
||||
print(format(ratio, "5.1e"))
|
||||
print('Mean swab-to-aersol vl ratio emission rate per hour:')
|
||||
print(format(ratio_1h, "5.1e"))
|
||||
|
||||
############# Coleman #############
|
||||
scatter_coleman_data(coleman_etal_vl_talking,
|
||||
coleman_etal_er_talking_2)
|
||||
|
||||
############ Legend ############
|
||||
build_talking_legend(fig)
|
||||
|
||||
elif expiration == 'Shouting':
|
||||
ratio_1h = np.mean(10**viral_loads / er_means_1h)
|
||||
print('Mean swab-to-aersol vl ratio emission rate per hour:')
|
||||
print(format(ratio_1h, "5.1e"))
|
||||
|
||||
ax.plot(viral_loads, er_means)
|
||||
ax.fill_between(viral_loads, lower_percentiles,
|
||||
upper_percentiles, alpha=0.2)
|
||||
ax.set_yscale('log')
|
||||
|
||||
ratio = np.mean(10**viral_loads / er_means)
|
||||
ratio_1h = np.mean(10**viral_loads / er_means_1h)
|
||||
print('Mean swab-to-aersol vl ratio in 30min:')
|
||||
print(format(ratio, "5.1e"))
|
||||
print('Mean swab-to-aersol vl ratio emission rate per hour:')
|
||||
print(format(ratio_1h, "5.1e"))
|
||||
|
||||
############# Coleman #############
|
||||
scatter_coleman_data(coleman_etal_vl_breathing,
|
||||
coleman_etal_er_breathing_2)
|
||||
|
||||
############# Milton et al #############
|
||||
scatter_milton_data(milton_vl, milton_er_2)
|
||||
|
||||
############# Yan et al #############
|
||||
scatter_yann_data(yann_vl, yann_er_2)
|
||||
|
||||
############ Legend ############
|
||||
build_breathing_legend(fig)
|
||||
|
||||
############ Plot ############
|
||||
plt.title('',
|
||||
fontsize=16, fontweight="bold")
|
||||
plt.ylabel(
|
||||
'Aerosol viral load, $\mathrm{vl_{out}}$\n(RNA copies)', fontsize=14)
|
||||
plt.xticks(ticks=[i for i in range(2, 13)], labels=[
|
||||
|
|
@ -91,10 +146,9 @@ def exposure_model_from_vl_breathing():
|
|||
|
||||
|
||||
""" Variation according to the BLO model """
|
||||
############ Breathing ############
|
||||
|
||||
|
||||
def exposure_model_from_vl_breathing_cn():
|
||||
def exposure_model_from_vl_cn(activity, expiration, mask):
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
|
||||
|
|
@ -105,14 +159,13 @@ def exposure_model_from_vl_breathing_cn():
|
|||
|
||||
for cn in tqdm(cns):
|
||||
er_means = np.array([])
|
||||
for vl in viral_loads:
|
||||
exposure_mc = breathing_exposure_vl(vl)
|
||||
for vl in tqdm(viral_loads):
|
||||
exposure_mc = exposure_vl_cn(
|
||||
activity, expiration, mask, vl, (cn, 0.2, 0.0010008))
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
# divide by 2 to have in 30min (half an hour)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
|
||||
cn_B=cn, cn_L=0.2) / 2
|
||||
emission_rate = emission_rate_when_present(exposure_model) / 2
|
||||
er_means = np.append(er_means, np.mean(emission_rate))
|
||||
|
||||
# divide by 2 to have in 30min (half an hour)
|
||||
coleman_etal_er_breathing_2 = [x/2 for x in coleman_etal_er_breathing]
|
||||
milton_er_2 = [x/2 for x in milton_er]
|
||||
|
|
@ -123,19 +176,21 @@ def exposure_model_from_vl_breathing_cn():
|
|||
# The dashed line for the chosen Cn,B
|
||||
er_means = np.array([])
|
||||
for vl in viral_loads:
|
||||
exposure_mc = breathing_exposure_vl(vl)
|
||||
exposure_mc = exposure_vl_cn(
|
||||
activity, expiration, mask, vl, (0.06, 0.2, 0.0010008))
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
# divide by 2 to have in 30min (half an hour)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
|
||||
cn_B=0.06, cn_L=0.2) / 2
|
||||
emission_rate = emission_rate_when_present(exposure_model) / 2
|
||||
er_means = np.append(er_means, np.mean(emission_rate))
|
||||
ax.plot(viral_loads, er_means, color=cmap.to_rgba(
|
||||
cn, alpha=0.75), linewidth=1, ls='--')
|
||||
plt.text(viral_loads[int(len(viral_loads)*0.9)], 10**4.2,
|
||||
r"$\mathbf{c_{n,B}=0.06}$", color=cmap.to_rgba(cn), fontsize=12)
|
||||
|
||||
plt.text(viral_loads[int(len(viral_loads)*0.9)], 10**4.2, r"$\mathbf{c_{n,B}=0.06}$", color=cmap.to_rgba(cn), fontsize=12) if activity == 'Breathing' else plt.text(
|
||||
viral_loads[int(len(viral_loads)*0.93)], 10**5.5, r"$\mathbf{c_{n,L}=0.2}$", color=cmap.to_rgba(cn), fontsize=12)
|
||||
|
||||
cmap = fig.colorbar(cmap, ticks=[0.01, 0.1, 0.25, 0.5])
|
||||
cmap.set_label(label='Particle emission concentration, ${c_{n,B}}$', fontsize=12)
|
||||
cmap.set_label(
|
||||
label='Particle emission concentration, ${c_{n,B}}$', fontsize=12)
|
||||
ax.set_yscale('log')
|
||||
|
||||
############# Coleman #############
|
||||
|
|
@ -161,189 +216,29 @@ def exposure_model_from_vl_breathing_cn():
|
|||
plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
|
||||
plt.show()
|
||||
|
||||
############ Talking ############
|
||||
|
||||
|
||||
def exposure_model_from_vl_talking():
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
|
||||
er_means = []
|
||||
er_means_1h = []
|
||||
er_medians = []
|
||||
lower_percentiles = []
|
||||
upper_percentiles = []
|
||||
|
||||
for vl in tqdm(viral_loads):
|
||||
exposure_mc = talking_exposure_vl(vl)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
# divide by 4 to have in 15min (quarter of an hour)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present()/4
|
||||
er_means.append(np.mean(emission_rate))
|
||||
er_medians.append(np.median(emission_rate))
|
||||
lower_percentiles.append(np.quantile(emission_rate, 0.01))
|
||||
upper_percentiles.append(np.quantile(emission_rate, 0.99))
|
||||
emission_rate_1h = exposure_model.concentration_model.infected.emission_rate_when_present()
|
||||
er_means_1h.append(np.mean(emission_rate_1h))
|
||||
|
||||
# divide by 4 to have in 15min (quarter of an hour)
|
||||
coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
|
||||
|
||||
ax.plot(viral_loads, er_means)
|
||||
ax.fill_between(viral_loads, lower_percentiles,
|
||||
upper_percentiles, alpha=0.2)
|
||||
ax.set_yscale('log')
|
||||
|
||||
ratio = np.mean(10**viral_loads / er_means)
|
||||
ratio_1h = np.mean(10**viral_loads / er_means_1h)
|
||||
print('Mean swab-to-aersol vl ratio in 30min:')
|
||||
print(format(ratio, "5.1e"))
|
||||
print('Mean swab-to-aersol vl ratio emission rate per hour:')
|
||||
print(format(ratio_1h, "5.1e"))
|
||||
|
||||
############# Coleman #############
|
||||
scatter_coleman_data(coleman_etal_vl_talking, coleman_etal_er_talking_2)
|
||||
|
||||
############ Legend ############
|
||||
build_talking_legend(fig)
|
||||
|
||||
############ Plot ############
|
||||
plt.title('',
|
||||
fontsize=16, fontweight="bold")
|
||||
plt.ylabel(
|
||||
'Aerosol viral load, $\mathrm{vl_{out}}$\n(RNA copies)', fontsize=14)
|
||||
plt.xticks(ticks=[i for i in range(2, 13)], labels=[
|
||||
'$10^{' + str(i) + '}$' for i in range(2, 13)])
|
||||
plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
|
||||
|
||||
plt.show()
|
||||
|
||||
def exposure_model_from_vl_talking_cn():
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
|
||||
n_lines = 30
|
||||
cns = np.linspace(0.01, 2, n_lines)
|
||||
cmap = define_colormap(cns)
|
||||
|
||||
for cn in tqdm(cns):
|
||||
er_means = np.array([])
|
||||
for vl in viral_loads:
|
||||
exposure_mc = talking_exposure_vl(vl)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
# divide by 4 to have in 15min (quarter of an hour)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
|
||||
cn_B=0.1, cn_L=cn) / 4
|
||||
er_means = np.append(er_means, np.mean(emission_rate))
|
||||
|
||||
# divide by 4 to have in 15min (quarter of an hour)
|
||||
coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
|
||||
ax.plot(viral_loads, er_means, color=cmap.to_rgba(
|
||||
cn, alpha=0.75), linewidth=0.5)
|
||||
|
||||
# The dashed line for the chosen Cn,L
|
||||
er_means = np.array([])
|
||||
for vl in viral_loads:
|
||||
exposure_mc = talking_exposure_vl(vl)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
# divide by 4 to have in 15min
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
|
||||
cn_B=0.06, cn_L=0.2) / 4
|
||||
er_means = np.append(er_means, np.mean(emission_rate))
|
||||
ax.plot(viral_loads, er_means, color=cmap.to_rgba(
|
||||
cn, alpha=0.75), linewidth=1, ls='--')
|
||||
plt.text(viral_loads[int(len(viral_loads)*0.93)], 10**5.5,
|
||||
r"$\mathbf{c_{n,L}=0.2}$", color=cmap.to_rgba(cn), fontsize=12)
|
||||
|
||||
cmap = fig.colorbar(cmap, ticks=[0.01, 0.5, 1.0, 2.0])
|
||||
cmap.set_label(label='Particle emission concentration, ${c_{n,L}}$', fontsize=12)
|
||||
ax.set_yscale('log')
|
||||
|
||||
############# Coleman #############
|
||||
scatter_coleman_data(coleman_etal_vl_talking, coleman_etal_er_talking_2)
|
||||
|
||||
############ Legend ############
|
||||
build_talking_legend(fig)
|
||||
|
||||
############ Plot ############
|
||||
plt.title('',
|
||||
fontsize=16, fontweight="bold")
|
||||
plt.ylabel(
|
||||
'Aerosol viral load, $\mathrm{vl_{out}}$\n(RNA copies)', fontsize=14)
|
||||
plt.xticks(ticks=[i for i in range(2, 13)], labels=[
|
||||
'$10^{' + str(i) + '}$' for i in range(2, 13)])
|
||||
plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
|
||||
plt.show()
|
||||
|
||||
####### Shouting ########
|
||||
|
||||
|
||||
def exposure_model_from_vl_shouting():
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
|
||||
er_means_1h = []
|
||||
|
||||
for vl in tqdm(viral_loads):
|
||||
exposure_mc = shouting_exposure_vl(vl)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
emission_rate_1h = exposure_model.concentration_model.infected.emission_rate_when_present()
|
||||
er_means_1h.append(np.mean(emission_rate_1h))
|
||||
|
||||
ax.plot(viral_loads, er_means_1h)
|
||||
ax.set_yscale('log')
|
||||
|
||||
ratio_1h = np.mean(10**viral_loads / er_means_1h)
|
||||
print('Mean swab-to-aersol vl ratio emission rate per hour:')
|
||||
print(format(ratio_1h, "5.1e"))
|
||||
|
||||
############ Plot ############
|
||||
plt.title('',
|
||||
fontsize=16, fontweight="bold")
|
||||
plt.ylabel(
|
||||
'Aerosol viral load, $\mathrm{vl_{out}}$\n(RNA copies)', fontsize=14)
|
||||
plt.xticks(ticks=[i for i in range(2, 13)], labels=[
|
||||
'$10^{' + str(i) + '}$' for i in range(2, 13)])
|
||||
plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14)
|
||||
plt.show()
|
||||
|
||||
|
||||
############ Plots with viral loads and emission rates ############
|
||||
############ Statistical Data ############
|
||||
|
||||
|
||||
def get_statistical_data(activity: str, mask: str):
|
||||
############ Breathing model ############
|
||||
exposure_mc = breathing_exposure(activity, mask)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present()
|
||||
breathing_er = [np.log10(er) for er in emission_rate]
|
||||
print('\n<<<<<<<<<<< Breathing model statistics >>>>>>>>>>>')
|
||||
print_er_info(emission_rate, breathing_er)
|
||||
|
||||
############ Speaking model ############
|
||||
exposure_mc = speaking_exposure(activity, mask)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present()
|
||||
speaking_er = [np.log10(er) for er in emission_rate]
|
||||
print('\n<<<<<<<<<<< Speaking model statistics >>>>>>>>>>>')
|
||||
print_er_info(emission_rate, speaking_er)
|
||||
|
||||
############ Shouting model ############
|
||||
exposure_mc = shouting_exposure(activity, mask)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present()
|
||||
shouting_er = [np.log10(er) for er in emission_rate]
|
||||
print('\n<<<<<<<<<<< Shouting model statistics >>>>>>>>>>>')
|
||||
print_er_info(emission_rate, shouting_er)
|
||||
log10_ers = {}
|
||||
for expiration in ('Breathing', 'Talking', 'Shouting'):
|
||||
exposure_mc = exposure_module(activity, expiration, mask)
|
||||
exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
|
||||
emission_rate = emission_rate_when_present(exposure_model)
|
||||
log10_ers[expiration] = [np.log10(er) for er in emission_rate]
|
||||
print('\n<<<<<<<<<<< ' + expiration + ' model statistics >>>>>>>>>>>')
|
||||
print_er_info(emission_rate, log10_ers[expiration])
|
||||
|
||||
viral_load_in_sputum = exposure_model.concentration_model.infected.virus.viral_load_in_sputum
|
||||
|
||||
return viral_load_in_sputum, breathing_er, speaking_er, shouting_er
|
||||
return viral_load_in_sputum, log10_ers['Breathing'], log10_ers['Talking'], log10_ers['Shouting']
|
||||
|
||||
|
||||
def present_vl_er_histograms(activity: str, mask: str):
|
||||
viral_load_in_sputum, breathing_er, speaking_er, shouting_er = get_statistical_data(activity, mask)
|
||||
viral_load_in_sputum, breathing_er, speaking_er, shouting_er = get_statistical_data(
|
||||
activity, mask)
|
||||
|
||||
fig, axs = plt.subplots(1, 2, sharex=False, sharey=False)
|
||||
fig.set_figheight(5)
|
||||
|
|
@ -355,29 +250,30 @@ def present_vl_er_histograms(activity: str, mask: str):
|
|||
|
||||
viral_loads = [np.log10(vl) for vl in viral_load_in_sputum]
|
||||
|
||||
axs[0].hist(viral_loads, bins = 300, color='lightgrey')
|
||||
axs[0].hist(viral_loads, bins=300, color='lightgrey')
|
||||
axs[0].set_xlabel('vl$_{\mathrm{in}}$ (log$_{10}$ RNA copies mL$^{-1}$)')
|
||||
|
||||
mean = np.mean(viral_loads)
|
||||
axs[0].vlines(x=(mean), ymin=0, ymax=axs[0].get_ylim()[1], colors=('grey'), linestyles=('dashed'))
|
||||
axs[0].vlines(x=(mean), ymin=0, ymax=axs[0].get_ylim()[
|
||||
1], colors=('grey'), linestyles=('dashed'))
|
||||
|
||||
breathing_mean_er = np.mean(breathing_er)
|
||||
speaking_mean_er = np.mean(speaking_er)
|
||||
shouting_mean_er = np.mean(shouting_er)
|
||||
|
||||
axs[1].hist(breathing_er, bins = 300, color='lightsteelblue')
|
||||
axs[1].hist(speaking_er, bins = 300, color='wheat')
|
||||
axs[1].hist(shouting_er, bins = 300, color='darkseagreen')
|
||||
axs[1].hist(breathing_er, bins=300, color='lightsteelblue')
|
||||
axs[1].hist(speaking_er, bins=300, color='wheat')
|
||||
axs[1].hist(shouting_er, bins=300, color='darkseagreen')
|
||||
axs[1].set_xlabel('vR (log$_{10}$ virions h$^{-1}$)')
|
||||
|
||||
axs[1].vlines(x=(breathing_mean_er, speaking_mean_er, shouting_mean_er), ymin=0, ymax=axs[1].get_ylim()[1], colors=('cornflowerblue', 'goldenrod', 'olivedrab'), alpha=(0.75, 0.75, 0.75), linestyles=('dashed', 'dashed', 'dashed'))
|
||||
axs[1].vlines(x=(breathing_mean_er, speaking_mean_er, shouting_mean_er), ymin=0, ymax=axs[1].get_ylim()[1], colors=(
|
||||
'cornflowerblue', 'goldenrod', 'olivedrab'), alpha=(0.75, 0.75, 0.75), linestyles=('dashed', 'dashed', 'dashed'))
|
||||
|
||||
labels = [Patch([], [], color=color, label=label)
|
||||
for color, label in zip(['lightgrey', 'lightsteelblue', 'wheat', 'darkseagreen'],
|
||||
['Viral Load', 'Breathing', 'Speaking', 'Shouting'])]
|
||||
for color, label in zip(['lightgrey', 'lightsteelblue', 'wheat', 'darkseagreen'],
|
||||
['Viral Load', 'Breathing', 'Speaking', 'Shouting'])]
|
||||
labels.append(mlines.Line2D([], [], color='black',
|
||||
marker='', linestyle='dashed', label='Mean'))
|
||||
|
||||
marker='', linestyle='dashed', label='Mean'))
|
||||
|
||||
for x in (0, 1):
|
||||
axs[x].set_yticklabels([])
|
||||
|
|
@ -394,39 +290,44 @@ def generate_cdf_curves():
|
|||
fig, axs = plt.subplots(3, 1, sharex='all')
|
||||
|
||||
############ Breathing models ############
|
||||
br_seated = breathing_seated_exposure()
|
||||
br_seated = exposure_module('Seated', 'Breathing', 'No mask')
|
||||
br_seated_model = br_seated.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
br_light_activity = breathing_light_activity_exposure()
|
||||
br_light_activity = exposure_module(
|
||||
'Light activity', 'Breathing', 'No mask')
|
||||
br_light_activity_model = br_light_activity.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
br_heavy_exercise = breathing_heavy_exercise_exposure()
|
||||
br_heavy_exercise = exposure_module(
|
||||
'Heavy exercise', 'Breathing', 'No mask')
|
||||
br_heavy_exercise_model = br_heavy_exercise.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
############ Speaking models ############
|
||||
sp_seated = speaking_seated_exposure()
|
||||
sp_seated = exposure_module('Seated', 'Talking', 'No mask')
|
||||
sp_seated_model = sp_seated.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
sp_light_activity = speaking_light_activity_exposure()
|
||||
sp_light_activity = exposure_module('Light activity', 'Talking', 'No mask')
|
||||
sp_light_activity_model = sp_light_activity.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
sp_heavy_exercise = speaking_heavy_exercise_exposure()
|
||||
sp_heavy_exercise = exposure_module('Heavy exercise', 'Talking', 'No mask')
|
||||
sp_heavy_exercise_model = sp_heavy_exercise.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
############ Shouting models ############
|
||||
sh_seated = shouting_seated_exposure()
|
||||
sh_seated = exposure_module('Seated', 'Shouting', 'No mask')
|
||||
sh_seated_model = sh_seated.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
sh_light_activity = shouting_light_activity_exposure()
|
||||
sh_light_activity = exposure_module(
|
||||
'Light activity', 'Shouting', 'No mask')
|
||||
sh_light_activity_model = sh_light_activity.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
sh_heavy_exercise = shouting_heavy_exercise_exposure()
|
||||
sh_heavy_exercise = exposure_module(
|
||||
'Heavy exercise', 'Shouting', 'No mask')
|
||||
sh_heavy_exercise_model = sh_heavy_exercise.build_model(size=SAMPLE_SIZE)
|
||||
|
||||
er_values = [np.log10(scenario.concentration_model.infected.emission_rate_when_present()) for scenario in (br_seated_model, br_light_activity_model,
|
||||
br_heavy_exercise_model, sp_seated_model, sp_light_activity_model, sp_heavy_exercise_model, sh_seated_model, sh_light_activity_model, sh_heavy_exercise_model)]
|
||||
left = min(np.min(ers) for ers in er_values)
|
||||
right = max(np.max(ers) for ers in er_values)
|
||||
er_values = [np.log10(emission_rate_when_present(scenario)) for scenario in (br_seated_model, br_light_activity_model,
|
||||
br_heavy_exercise_model, sp_seated_model,
|
||||
sp_light_activity_model, sp_heavy_exercise_model,
|
||||
sh_seated_model, sh_light_activity_model,
|
||||
sh_heavy_exercise_model)]
|
||||
|
||||
# Colors can be changed here
|
||||
colors_breathing = ['lightsteelblue', 'cornflowerblue', 'royalblue']
|
||||
|
|
@ -444,12 +345,11 @@ def generate_cdf_curves():
|
|||
for color, label in zip(colors_shouting, breathing_rates)]
|
||||
lines = [lines_breathing, lines_speaking, lines_shouting]
|
||||
|
||||
samples: int = 250000
|
||||
for i in range(3):
|
||||
axs[i].hist(er_values[3 * i:3 * (i + 1)], bins=2000,
|
||||
histtype='step', cumulative=True, range=(-7, 6), color=colors[i])
|
||||
axs[i].set_xlim(-6, 6)
|
||||
axs[i].set_yticks([0, samples / 2, samples])
|
||||
axs[i].set_yticks([0, SAMPLE_SIZE / 2, SAMPLE_SIZE])
|
||||
axs[i].set_yticklabels(['0.0', '0.5', '1.0'])
|
||||
axs[i].yaxis.set_label_position("right")
|
||||
axs[i].set_ylabel(activities[i], fontsize=14)
|
||||
|
|
@ -467,13 +367,15 @@ def generate_cdf_curves():
|
|||
fig.set_figwidth(5)
|
||||
plt.show()
|
||||
|
||||
############ Deposition Fraction Graph #############
|
||||
############ Deposition Fraction Graph #############
|
||||
|
||||
|
||||
def calculate_deposition_factor():
|
||||
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(1, 1, 1)
|
||||
|
||||
diameters = np.linspace(0.001, 0.01, 1000) #particle diameter (μm)
|
||||
|
||||
diameters = np.linspace(0.001, 0.01, 1000) # particle diameter (μm)
|
||||
diameters = np.append(diameters, np.linspace(0.01, 0.1, 1000))
|
||||
diameters = np.append(diameters, np.linspace(0.1, 1., 1000))
|
||||
diameters = np.append(diameters, np.linspace(1., 10., 1000))
|
||||
|
|
@ -484,31 +386,38 @@ def calculate_deposition_factor():
|
|||
fractions_al = []
|
||||
fractions_df = []
|
||||
for d in diameters:
|
||||
|
||||
|
||||
IF = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8)))))
|
||||
DF_et = IF * (
|
||||
(1 / ( 1 + np.exp(6.84 + 1.183 * np.log(d)))) + (1 / (1 + np.exp(0.924 - 1.885 * np.log(d))))
|
||||
(1 / (1 + np.exp(6.84 + 1.183 * np.log(d)))) +
|
||||
(1 / (1 + np.exp(0.924 - 1.885 * np.log(d))))
|
||||
)
|
||||
fractions_et.append(DF_et)
|
||||
|
||||
DF_tb = (0.00352/d) * (np.exp(-0.234*((np.log(d) + 3.40)**2)) + (63.9 * np.exp(-0.819*((np.log(d) - 1.61)**2))))
|
||||
DF_tb = (0.00352/d) * (np.exp(-0.234*((np.log(d) + 3.40)**2)
|
||||
) + (63.9 * np.exp(-0.819*((np.log(d) - 1.61)**2))))
|
||||
fractions_tb.append(DF_tb)
|
||||
|
||||
DF_al = (0.0155/d) * (np.exp(-0.416*((np.log(d) + 2.84)**2)) + (19.11*np.exp(-0.482 * ((np.log(d) - 1.362)**2))))
|
||||
DF_al = (0.0155/d) * (np.exp(-0.416*((np.log(d) + 2.84)**2)) +
|
||||
(19.11*np.exp(-0.482 * ((np.log(d) - 1.362)**2))))
|
||||
fractions_al.append(DF_al)
|
||||
|
||||
DF = IF * (0.0587 + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d)))))
|
||||
DF = IF * (0.0587 + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))
|
||||
) + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d)))))
|
||||
fractions_df.append(DF)
|
||||
|
||||
ax.plot(diameters, fractions_df, label = 'Total Deposition', color='k')
|
||||
ax.plot(diameters, fractions_et, label = 'Extrathoracic', ls='-.', lw=0.9, color='grey')
|
||||
ax.plot(diameters, fractions_tb, label = 'Tracheobronchial', ls='--', lw=0.9, color='darkgray')
|
||||
ax.plot(diameters, fractions_al, label = 'Alveolar', ls=(0,(1,1)), lw=0.9, color='darkgray')
|
||||
|
||||
ax.plot(diameters, fractions_df, label='Total Deposition', color='k')
|
||||
ax.plot(diameters, fractions_et, label='Extrathoracic',
|
||||
ls='-.', lw=0.9, color='grey')
|
||||
ax.plot(diameters, fractions_tb, label='Tracheobronchial',
|
||||
ls='--', lw=0.9, color='darkgray')
|
||||
ax.plot(diameters, fractions_al, label='Alveolar',
|
||||
ls=(0, (1, 1)), lw=0.9, color='darkgray')
|
||||
|
||||
ax.grid(linestyle='--')
|
||||
ax.set_xscale('log')
|
||||
ax.margins(x=0, y=0)
|
||||
plt.legend(bbox_to_anchor=(1.04,1), loc="upper left")
|
||||
plt.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
|
||||
|
||||
y_ticks = [0., 0.2, 0.4, 0.6, 0.8, 1]
|
||||
x_ticks = [0.001, 0.01, 0.1, 1, 10, 100]
|
||||
|
|
@ -518,38 +427,46 @@ def calculate_deposition_factor():
|
|||
str(i) for i in x_ticks])
|
||||
plt.xlabel('Particle diameter (μm)', fontsize=14)
|
||||
plt.ylabel('Deposition fraction\nf$_{dep}$', fontsize=14)
|
||||
|
||||
|
||||
fig.set_figwidth(10)
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
############ Compare concentration curves ############
|
||||
|
||||
|
||||
def compare_concentration_curves():
|
||||
|
||||
exp_models=[office_model_no_mask_windows_closed().build_model(size=SAMPLE_SIZE),
|
||||
office_model_no_mask_windows_open_breaks().build_model(size=SAMPLE_SIZE),
|
||||
office_model_no_mask_windows_open_alltimes().build_model(size=SAMPLE_SIZE)]
|
||||
exp_models = [office_model_no_mask_windows_closed().build_model(size=SAMPLE_SIZE),
|
||||
office_model_no_mask_windows_open_breaks().build_model(size=SAMPLE_SIZE),
|
||||
office_model_no_mask_windows_open_alltimes().build_model(size=SAMPLE_SIZE)]
|
||||
|
||||
labels=['Windows closed', 'Window open during breaks', 'Window open at all times']
|
||||
colors=['tomato', 'lightskyblue', 'limegreen', '#1f77b4', 'seagreen', 'lightskyblue', 'deepskyblue']
|
||||
labels = ['Windows closed', 'Window open during breaks',
|
||||
'Window open at all times']
|
||||
colors = ['tomato', 'lightskyblue', 'limegreen',
|
||||
'#1f77b4', 'seagreen', 'lightskyblue', 'deepskyblue']
|
||||
|
||||
start=min(min(model.concentration_model.infected.presence.transition_times()) for model in exp_models)
|
||||
stop=max(max(model.concentration_model.infected.presence.transition_times()) for model in exp_models)
|
||||
start = min(min(model.concentration_model.infected.presence.transition_times())
|
||||
for model in exp_models)
|
||||
stop = max(max(model.concentration_model.infected.presence.transition_times())
|
||||
for model in exp_models)
|
||||
|
||||
TIMESTEP = 0.01
|
||||
times=np.arange(start, stop, TIMESTEP)
|
||||
|
||||
concentrations = [[np.mean(model.concentration_model.concentration(t)) for t in times] for model in exp_models]
|
||||
times = np.arange(start, stop, TIMESTEP)
|
||||
|
||||
concentrations = [[np.mean(model.concentration_model.concentration(
|
||||
t)) for t in times] for model in exp_models]
|
||||
fig, ax = plt.subplots()
|
||||
for c, label, color in zip(concentrations, labels, colors):
|
||||
ax.plot(times, c, label=label, color=color)
|
||||
|
||||
|
||||
ax.legend(loc='upper left')
|
||||
ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2)
|
||||
ax.spines["right"].set_visible(False)
|
||||
|
||||
cumulative_doses = [np.cumsum([
|
||||
np.array(model.exposure_between_bounds(float(time1), float(time2))).mean()
|
||||
np.array(exposure_between_bounds(model,
|
||||
float(time1), float(time2))).mean()
|
||||
for time1, time2 in zip(times[:-1], times[1:])
|
||||
]) for model in exp_models]
|
||||
|
||||
|
|
@ -558,85 +475,97 @@ def compare_concentration_curves():
|
|||
|
||||
ax1 = ax.twinx()
|
||||
for qd, label, color in zip(cumulative_doses, labels, colors):
|
||||
ax1.plot(times[:-1], qd, label='qD - ' + label, color=color, linestyle='dotted')
|
||||
ax1.plot(times[:-1], qd, label='qD - ' + label,
|
||||
color=color, linestyle='dotted')
|
||||
ax1.spines["right"].set_linestyle("--")
|
||||
ax1.spines["right"].set_linestyle((0,(1,5)))
|
||||
ax1.spines["right"].set_linestyle((0, (1, 5)))
|
||||
ax1.set_ylabel('Mean cumulative dose (virions)', fontsize=14)
|
||||
ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.2)
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
def compare_viruses_vr():
|
||||
|
||||
# Represented as tuples of three numbers on the interval [0, 1] (e.g. (1, 0, 0)) (R, G, B)
|
||||
colors = [(0., 0.5, 0.5), (0, 0, 0.5), (0.5, 0, 0)]
|
||||
colors_violin=['lightsteelblue', 'wheat', 'darkseagreen']
|
||||
colors_violin = ['lightsteelblue', 'wheat', 'darkseagreen']
|
||||
colors_violin_lines = ['royalblue', 'orange', 'forestgreen']
|
||||
|
||||
# The colors of the borders surrounding the violin plots
|
||||
border_colors = [(0, 0, 0), (0, 0, 0), (0, 0, 0)]
|
||||
|
||||
|
||||
whisker_width = 0.8
|
||||
positions = [1, 2, 3, 12]
|
||||
|
||||
infected_populations = [infected_model('No mask', 'Light activity', activity).build_model(size=SAMPLE_SIZE) for activity in ('Breathing', 'Talking', 'Shouting')]
|
||||
|
||||
vrs = [np.log10(pop.emission_rate_when_present()) for pop in infected_populations]
|
||||
exposure_models = [exposure_module('Light activity', expiration, 'No mask').build_model(
|
||||
size=SAMPLE_SIZE) for expiration in ('Breathing', 'Talking', 'Shouting')]
|
||||
|
||||
vrs = [np.log10(emission_rate_when_present(pop))
|
||||
for pop in exposure_models]
|
||||
|
||||
fig, ax = plt.subplots()
|
||||
ax.set_xlim((0, 11))
|
||||
|
||||
parts = ax.violinplot(vrs, quantiles=[(0.05, 0.95) for _ in vrs], showextrema=False)
|
||||
parts = ax.violinplot(
|
||||
vrs, quantiles=[(0.05, 0.95) for _ in vrs], showextrema=False)
|
||||
means = [np.log10(np.mean(10 ** vr)) for vr in vrs]
|
||||
ax.hlines(y=means,
|
||||
xmin=[pos - whisker_width / 2 for pos in positions[:3]],
|
||||
xmax=[pos + whisker_width / 2 for pos in positions[:3]],
|
||||
colors=colors_violin_lines,
|
||||
alpha=0.8)
|
||||
|
||||
xmin=[pos - whisker_width / 2 for pos in positions[:3]],
|
||||
xmax=[pos + whisker_width / 2 for pos in positions[:3]],
|
||||
colors=colors_violin_lines,
|
||||
alpha=0.8)
|
||||
|
||||
for pc, color, bc in zip(parts['bodies'], colors_violin, border_colors):
|
||||
pc.set_facecolor(color)
|
||||
pc.set_edgecolor(bc)
|
||||
pc.set_alpha(0.5)
|
||||
parts['cquantiles'].set_color([c for c in colors_violin_lines[:3] for _ in range(2)])
|
||||
parts['cquantiles'].set_color(
|
||||
[c for c in colors_violin_lines[:3] for _ in range(2)])
|
||||
parts['cquantiles'].set_alpha(1)
|
||||
|
||||
positions=np.linspace(4.5, 11.5, 20)
|
||||
positions = np.linspace(4.5, 11.5, 20)
|
||||
|
||||
######### SARS-CoV-2 #########
|
||||
lower_bound = [418, 216, 216, 518, 648, 878, 893, 1670, 1872, 1915, 2002, 2002, 2189, 3341, 9835, 13968, 60667]
|
||||
higher_bound = [4176, 2160, 2160, 5184, 6480, 8784, 8928, 16704, 18720, 19152, 20016, 20016, 21888, 33408, 98352, 139680, 606672]
|
||||
|
||||
lower_bound = [418, 216, 216, 518, 648, 878, 893, 1670,
|
||||
1872, 1915, 2002, 2002, 2189, 3341, 9835, 13968, 60667]
|
||||
higher_bound = [4176, 2160, 2160, 5184, 6480, 8784, 8928, 16704,
|
||||
18720, 19152, 20016, 20016, 21888, 33408, 98352, 139680, 606672]
|
||||
|
||||
for i in range(len(lower_bound)):
|
||||
data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
|
||||
ax.boxplot(np.log10(data), positions=[positions[i]], medianprops=dict(color=colors[0]+ (0.5,)), whiskerprops=dict(color=colors[0]+ (0.5,)), boxprops=dict(color=colors[0]+ (0.5,)))
|
||||
ax.boxplot(np.log10(data), positions=[positions[i]], medianprops=dict(
|
||||
color=colors[0] + (0.5,)), whiskerprops=dict(color=colors[0] + (0.5,)), boxprops=dict(color=colors[0] + (0.5,)))
|
||||
|
||||
######### Measles #########
|
||||
lower_bound = [259, 8640, 39816, 124416]
|
||||
higher_bound = [2592, 86400, 398160, 1244160]
|
||||
|
||||
|
||||
for i in range(len(lower_bound)):
|
||||
data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
|
||||
ax.boxplot(np.log10(data), positions=[positions[i+5]], medianprops=dict(color=colors[1]+ (0.5,)), whiskerprops=dict(color=colors[1]+ (0.5,)), boxprops=dict(color=colors[1]+ (0.5,)))
|
||||
ax.boxplot(np.log10(data), positions=[positions[i+5]], medianprops=dict(color=colors[1] + (
|
||||
0.5,)), whiskerprops=dict(color=colors[1] + (0.5,)), boxprops=dict(color=colors[1] + (0.5,)))
|
||||
|
||||
######### Influenza #########
|
||||
lower_bound = [2, 114, 1138]
|
||||
higher_bound = [16, 1145, 11376]
|
||||
|
||||
|
||||
for i in range(len(lower_bound)):
|
||||
data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
|
||||
ax.boxplot(np.log10(data), positions=[positions[i+12]], medianprops=dict(color=colors[2]+ (0.5,)), whiskerprops=dict(color=colors[2]+ (0.5,)), boxprops=dict(color=colors[2]+ (0.5,)))
|
||||
ax.boxplot(np.log10(data), positions=[positions[i+12]], medianprops=dict(color=colors[2] + (
|
||||
0.5,)), whiskerprops=dict(color=colors[2] + (0.5,)), boxprops=dict(color=colors[2] + (0.5,)))
|
||||
|
||||
######### Rhinovirus #########
|
||||
lower_bound = [45]
|
||||
higher_bound = [446]
|
||||
|
||||
|
||||
for i in range(len(lower_bound)):
|
||||
data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
|
||||
ax.boxplot(np.log10(data), positions=[positions[i+8]], medianprops=dict(color=(0.5, 0.5, 0.5, 0.5, )), whiskerprops=dict(color=(0.5, 0.5, 0.5, 0.5,)), boxprops=dict(color=(0.5, 0.5, 0.5, 0.5,)))
|
||||
|
||||
handles = [patches.Patch(edgecolor=c, facecolor='none', label=l) for c, l in zip([p + (0.5,) for p in [(0., 0.5, 0.5), (0, 0, 0.5), (0.5, 0, 0), (0.5, 0.5, 0.5)]], ('SARS-CoV-2', 'Measles', 'Influenza', 'Rhinovirus'))]
|
||||
ax.boxplot(np.log10(data), positions=[positions[i+8]], medianprops=dict(color=(0.5, 0.5, 0.5, 0.5, )),
|
||||
whiskerprops=dict(color=(0.5, 0.5, 0.5, 0.5,)), boxprops=dict(color=(0.5, 0.5, 0.5, 0.5,)))
|
||||
|
||||
handles = [patches.Patch(edgecolor=c, facecolor='none', label=l) for c, l in zip(
|
||||
[p + (0.5,) for p in [(0., 0.5, 0.5), (0, 0, 0.5), (0.5, 0, 0), (0.5, 0.5, 0.5)]], ('SARS-CoV-2', 'Measles', 'Influenza', 'Rhinovirus'))]
|
||||
boxplot_legend = plt.legend(handles=handles, loc='lower right')
|
||||
|
||||
ax.annotate("Bus ride", xy=(6, np.log10(5500)), color='k', fontsize=8,
|
||||
|
|
@ -650,17 +579,21 @@ def compare_viruses_vr():
|
|||
xytext=(-50, 40), textcoords='offset points',
|
||||
arrowprops=dict(arrowstyle="->",
|
||||
connectionstyle="arc3,rad=-0.2", color='lightgrey'))
|
||||
|
||||
handles = [patches.Patch(color=c, label=l) for c, l in zip([p for p in colors_violin], ('Breathing', 'Speaking', 'Shouting'))]
|
||||
|
||||
handles = [patches.Patch(color=c, label=l) for c, l in zip(
|
||||
[p for p in colors_violin], ('Breathing', 'Speaking', 'Shouting'))]
|
||||
plt.legend(handles=handles, loc='lower left', bbox_to_anchor=(0.12, 0.))
|
||||
plt.gca().add_artist(boxplot_legend)
|
||||
|
||||
ax.hlines(y=[-2, 0, 2, 4, 6], xmin=ax.get_xlim()[0], xmax=ax.get_xlim()[1], colors=(0.8, 0.8, 0.8), linestyles='--', alpha=0.3)
|
||||
ax.vlines(x=4, ymin=ax.get_ylim()[0], ymax=ax.get_ylim()[1], colors=(0.8, 0.8, 0.8))
|
||||
ax.hlines(y=[-2, 0, 2, 4, 6], xmin=ax.get_xlim()[0], xmax=ax.get_xlim()
|
||||
[1], colors=(0.8, 0.8, 0.8), linestyles='--', alpha=0.3)
|
||||
ax.vlines(x=4, ymin=ax.get_ylim()[0], ymax=ax.get_ylim()[
|
||||
1], colors=(0.8, 0.8, 0.8))
|
||||
ax.set_yticks([i for i in range(-4, 7, 2)])
|
||||
ax.set_yticklabels(['$10^{' + str(i) + '}$' for i in range(-4, 7, 2)])
|
||||
ax.set_xticks([2, 7])
|
||||
ax.set_xticklabels(['SARS-CoV-2\n(model)', 'Literature Data\n(recorded outbreaks) '])
|
||||
ax.set_xticklabels(
|
||||
['SARS-CoV-2\n(model)', 'Literature Data\n(recorded outbreaks) '])
|
||||
ax.set_ylabel('Emission rate (virions h$^{-1}$)', fontsize=12)
|
||||
|
||||
plt.tight_layout()
|
||||
|
|
@ -668,6 +601,7 @@ def compare_viruses_vr():
|
|||
|
||||
######### Auxiliar functions #########
|
||||
|
||||
|
||||
def get_enclosure_points(x_coordinates, y_coordinates):
|
||||
df = pd.DataFrame({'x': x_coordinates, 'y': y_coordinates})
|
||||
|
||||
|
|
@ -799,6 +733,3 @@ def print_er_info(er: np.array, log_er: np.array):
|
|||
print(f"ER_{quantile} = {np.quantile(er, quantile)}")
|
||||
|
||||
return
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue