diff --git a/cara/model_scenarios_paper.py b/cara/model_scenarios_paper.py index 89d65fe2..4e70e5f8 100644 --- a/cara/model_scenarios_paper.py +++ b/cara/model_scenarios_paper.py @@ -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 \ No newline at end of file diff --git a/cara/plot_output.py b/cara/plot_output.py index c63f27da..b0e82305 100644 --- a/cara/plot_output.py +++ b/cara/plot_output.py @@ -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() diff --git a/cara/results_paper.py b/cara/results_paper.py index 54c74bb3..ce0571fd 100644 --- a/cara/results_paper.py +++ b/cara/results_paper.py @@ -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 - - -