diff --git a/cara/short_range_plots/model_scenarios.py b/cara/short_range_plots/model_scenarios.py index dbe0ddd7..dee31201 100644 --- a/cara/short_range_plots/model_scenarios.py +++ b/cara/short_range_plots/model_scenarios.py @@ -62,11 +62,11 @@ def exposure_module_without_short_range(activity: str, expiration: str, mask: st room=models.Room(volume=100, humidity=0.5), ventilation=models.AirChange( active=models.SpecificInterval(((0, 24),)), - air_exch=0.25, + air_exch=500/100, ), infected=mc.InfectedPopulation( number=1, - virus=virus_distributions['SARS_CoV_2'], + virus=virus_distributions['SARS_CoV_2_OMICRON'], presence=mc.SpecificInterval(((8.5, 12.5),(13.5, 17.5),)), mask=exposure_mask, activity=activity_distributions[activity], @@ -89,11 +89,7 @@ def exposure_module_without_short_range(activity: str, expiration: str, mask: st ) return exposure_mc -def exposure_module_with_short_range(activity: str, - expiration: str, - mask: str, - sr_presence: list, - sr_activities: list): +def exposure_module_with_short_range(activity: str, expiration: str, mask: str, sr_presence: list, sr_activities: list): if mask == 'No mask': exposure_mask = models.Mask.types['No mask'] else: @@ -104,11 +100,11 @@ def exposure_module_with_short_range(activity: str, room=models.Room(volume=100, humidity=0.5), ventilation=models.AirChange( active=models.SpecificInterval(((0, 24),)), - air_exch=0.25, + air_exch=500/100, ), infected=mc.InfectedPopulation( number=1, - virus=virus_distributions['SARS_CoV_2'], + virus=virus_distributions['SARS_CoV_2_OMICRON'], presence=mc.SpecificInterval(((8.5, 12.5),(13.5, 17.5),)), mask=exposure_mask, activity=activity_distributions[activity], @@ -132,12 +128,7 @@ def exposure_module_with_short_range(activity: str, ) return exposure_mc - -def baseline_model(activity: str, - expiration: str, - mask: str, - sr_presence: list, - sr_activities: list): +def exposure_module_with_short_range_outdoors(activity: str, expiration: str, mask: str, sr_presence: list, sr_activities: list): if mask == 'No mask': exposure_mask = models.Mask.types['No mask'] else: @@ -145,7 +136,47 @@ def baseline_model(activity: str, exposure_mc = mc.ExposureModel( concentration_model=mc.ConcentrationModel( - room=models.Room(volume=10, humidity=0.5), + room=models.Room(volume=100000, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0, 24),)), + air_exch=1000000, + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_OMICRON'], + presence=mc.SpecificInterval(((8.5, 12.5),)), + mask=exposure_mask, + activity=activity_distributions[activity], + expiration=build_expiration(expiration), + host_immunity=0., + ), + ), + short_range=mc.ShortRangeModel( + presence=[models.SpecificInterval(interval) for interval in sr_presence], + expirations=[short_range_expiration_distributions[activity] for activity in sr_activities], + dilutions=dilution_factor(activities=sr_activities, + distance=np.random.uniform(0.5, 1.5, 250000)), + ), + exposed=mc.Population( + number=14, + presence=mc.SpecificInterval(((8.5, 12.5),)), + activity=activity_distributions[activity], + mask=exposure_mask, + host_immunity=0., + ), + ) + return exposure_mc + + +def baseline_model(activity: str, expiration: str, mask: str, sr_presence: list, sr_activities: list): + if mask == 'No mask': + exposure_mask = models.Mask.types['No mask'] + else: + exposure_mask = mask_distributions[mask] + + exposure_mc = mc.ExposureModel( + concentration_model=mc.ConcentrationModel( + room=models.Room(volume=100, humidity=0.5), ventilation=models.AirChange( active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25, @@ -153,7 +184,7 @@ def baseline_model(activity: str, infected=mc.InfectedPopulation( number=1, virus=virus_distributions['SARS_CoV_2_OMICRON'], - presence=models.SpecificInterval(((8.5, 12.5),)), + presence=models.SpecificInterval(((8.5, 9.5),)), mask=exposure_mask, activity=activity_distributions[activity], expiration=build_expiration(expiration), @@ -168,7 +199,7 @@ def baseline_model(activity: str, ), exposed=mc.Population( number=3, - presence=models.SpecificInterval(((8.5, 12.5),)), + presence=models.SpecificInterval(((8.5, 9.5),)), activity=activity_distributions[activity], mask=exposure_mask, host_immunity=0., diff --git a/cara/short_range_plots/results.py b/cara/short_range_plots/results.py index 658ab820..955f7d8a 100644 --- a/cara/short_range_plots/results.py +++ b/cara/short_range_plots/results.py @@ -13,48 +13,70 @@ from dataclasses import dataclass from cara.monte_carlo.data import symptomatic_vl_frequencies # print('\n<<<<<<<<<<< Peak viral concentration without short range for baseline scenarios >>>>>>>>>>>') -# concentration_curve(models = [exposure_module_without_short_range(activity='Seated', expiration='Breathing', mask='No mask')], -# labels = ['Baseline'], +# concentration_curve(models = [exposure_module_without_short_range(activity='Light activity', expiration={"Speaking": 1, "Breathing": 2}, mask='No mask')], +# labels = ['Background (long-range) concentration'], +# labelsDose = ['Dose (long-range)'], # colors = ['royalblue'], -# ) +# linestyles = ['--'], +# thickness = [2]) # print('\n<<<<<<<<<<< Peak viral concentration with short range interactions for baseline scenarios >>>>>>>>>>>') # concentration_curve(models=[exposure_module_with_short_range( -# activity='Seated', -# expiration={"Speaking": 1, "Breathing": 2}, +# activity='Light activity', +# expiration={"Speaking": 1, "Breathing": 2}, # mask='No mask', -# sr_presence=[(10.5, 11.0), (15.5, 16.0)], -# sr_activities=['Breathing', 'Shouting']), +# sr_presence=[(10.5, 11.0)], +# sr_activities=['Breathing']), # exposure_module_without_short_range( -# activity='Seated', -# expiration={"Speaking": 1, "Breathing": 2}, -# mask='No mask', -# )], -# labels = ['Mean concentration with short range interactions', 'Mean concentration without short range interactions'], -# colors = ['royalblue', 'darkviolet'], -# linestyles = ['-', '-'], -# thickness = [2, 2.5]) +# activity='Light activity', +# expiration={"Speaking": 1, "Breathing": 2}, +# mask='No mask',) +# ], +# labels = ['Concentration with short range interactions', 'Background (long-range) concentration'], +# labelsDose = ['Dose (full)', 'Dose (long-range)'], +# colors = ['salmon', 'royalblue'], +# linestyles = ['-', '--'], +# thickness = [2, 2]) -# print('\n<<<<<<<<<<< Probability of infections vs exposure time >>>>>>>>>>>') -# Always assume 1h for the short range interactions. -# Always assume that in each model there is only ONE short range interaction. -# plot_pi_vs_exposure_time(exp_models = [ -# baseline_model( -# activity='Seated', -# expiration={"Speaking": 2, "Breathing": 1}, -# mask='No mask', -# sr_presence=[(10.0, 11.0)], -# sr_activities=['Breathing']), -# baseline_model( -# activity='Seated', -# expiration={"Shouting": 2, "Breathing": 1}, -# mask='No mask', -# sr_presence=[(10.0, 11.0)], -# sr_activities=['Shouting'])], -# labels = ['Baseline model breathing', 'Baseline model shouting'], -# colors=['royalblue', 'darkviolet'], -# linestyles=['solid', 'solid'], -# points=20, -# time_in_minutes=True, -# normalize_y_axis=True) +# print('\n<<<<<<<<<<< Peak viral concentration with short range interactions for baseline scenarios >>>>>>>>>>>') +# concentration_curve(models=[exposure_module_with_short_range( +# activity='Light activity', +# expiration={"Speaking": 1, "Breathing": 2}, +# mask='No mask', +# sr_presence=[(10.5, 11.0), (15.0, 16.0)], +# sr_activities=['Breathing', 'Speaking']), +# exposure_module_without_short_range( +# activity='Light activity', +# expiration={"Speaking": 1, "Breathing": 2}, +# mask='No mask',) +# ], +# labels = ['Concentration with short range interactions', 'Background (long-range) concentration'], +# labelsDose = ['Dose (full)', 'Dose (long-range)'], +# colors = ['salmon', 'royalblue'], +# linestyles = ['-', '--'], +# thickness = [2, 2]) + + +print('\n<<<<<<<<<<< Dose vs SR exposure time >>>>>>>>>>>') +#Always assume 1h for the short range interactions. +#Always assume that in each model there is only ONE short range interaction. +plot_vD_vs_exposure_time(exp_models = [ + baseline_model( + activity='Light activity', + expiration={"Speaking": 2, "Breathing": 1}, + mask='No mask', + sr_presence=[(8.5, 9.5)], + sr_activities=['Breathing']), + baseline_model( + activity='Light activity', + expiration={"Speaking": 2, "Breathing": 1}, + mask='No mask', + sr_presence=[(8.5, 9.5)], + sr_activities=['Speaking'])], + labels = ['Breathing', 'Speaking'], + colors=['royalblue', 'darkviolet'], + linestyles=['solid', 'solid'], + points=20, + time_in_minutes=True, + normalize_y_axis=True) diff --git a/cara/short_range_plots/scripts.py b/cara/short_range_plots/scripts.py index de456b0b..e6f7d6de 100644 --- a/cara/short_range_plots/scripts.py +++ b/cara/short_range_plots/scripts.py @@ -23,7 +23,7 @@ from mpl_toolkits.axes_grid1.inset_locator import mark_inset np.random.seed(2000) SAMPLE_SIZE = 250000 TIMESTEP = 0.01 -viral_loads = np.linspace(2, 12, 600) +#viral_loads = np.linspace(2, 12, 600) _VectorisedFloat = typing.Union[float, np.ndarray] @@ -59,7 +59,7 @@ def previous_deposited_exposure_between_bounds(model: ExposureModel, time1: floa (1 - model.exposed.mask.inhale_efficiency())) -def concentration_curve(models, labels, colors, linestyles, thickness): +def concentration_curve(models, labels, labelsDose, colors, linestyles, thickness): exp_models = [model.build_model(size=SAMPLE_SIZE) for model in models] @@ -73,12 +73,13 @@ def concentration_curve(models, labels, colors, linestyles, thickness): concentrations = [[np.mean(model.concentration( t)) for t in times] for model in tqdm(exp_models)] - fig, ax = plt.subplots() + fig, ax = plt.subplots(figsize=(8,6)) for c, color, linestyle, width in zip(concentrations, colors, linestyles, thickness): ax.plot(times, c, color=color, ls=linestyle, lw=width) - ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2) + #ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2) + ax.set_ylim(ax.get_ylim()[0], 90) ax.spines["right"].set_visible(False) cumulative_doses = [np.cumsum([ @@ -129,10 +130,11 @@ def concentration_curve(models, labels, colors, linestyles, thickness): 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.3) + ax1.set_ylabel('Mean cumulative dose\n(infectious virus)', fontsize=14) + #ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.3) + ax1.set_ylim(ax1.get_ylim()[0], 40) - complete_labels = labels + ['vD - ' + label for label in labels] + complete_labels = labels + [label for label in labelsDose] complete_colors = colors + [color for color in colors] complete_linestyles = linestyles + ['dotted' for linestyle in linestyles] @@ -146,21 +148,17 @@ def concentration_curve(models, labels, colors, linestyles, thickness): f"5th per = {quantile_05[i][-1]}\n" f"95th per = {quantile_95[i][-1]}\n") - plt.legend(handles=labels_legend, loc='upper right', bbox_to_anchor=(1, 1)) + plt.legend(handles=labels_legend, loc='upper left') plt.show() -def plot_pi_vs_exposure_time(exp_models: typing.List[mc.ExposureModel], - labels, - colors, - linestyles, - points: int = 20, time_in_minutes: bool = False, normalize_y_axis: bool = False) -> None: +def plot_vD_vs_exposure_time(exp_models: typing.List[mc.ExposureModel], labels, colors, linestyles, points: int = 20, time_in_minutes: bool = False, normalize_y_axis: bool = False) -> None: TIMESTEP = 0.001 concentration_models = [model.concentration_model for model in exp_models] exposed_models = [model.exposed for model in exp_models] - pis: typing.List[typing.List[float]] = [[] for _ in exp_models] + vDs: typing.List[typing.List[float]] = [[] for _ in exp_models] presence_intervals = [model.short_range.presence[0].boundaries() for model in exp_models] start, final = presence_intervals[0] @@ -177,20 +175,20 @@ def plot_pi_vs_exposure_time(exp_models: typing.List[mc.ExposureModel], ) for cm, exposed, em in zip(concentration_models, exposed_models, exp_models)] for i, m in enumerate(current_models): - pis[i].append(np.mean(m.build_model(SAMPLE_SIZE).infection_probability() / 100)) + vDs[i].append(np.mean(m.build_model(SAMPLE_SIZE).deposited_exposure())) times = np.linspace(0, 60, points) - for i, pi in enumerate(pis): - plt.plot(times, pi, color=colors[i], label=labels[i]) + for i, vD in enumerate(vDs): + plt.plot(times, vD, color=colors[i], label=labels[i]) # plt.xlim((0, 60)) - if normalize_y_axis: - plt.ylim((0, 1)) + #if normalize_y_axis: + # plt.ylim((0, 1)) for m in exp_models: - print(np.mean(m.build_model(SAMPLE_SIZE).infection_probability() / 100)) + print(np.mean(m.build_model(SAMPLE_SIZE).deposited_exposure())) - plt.xlabel(f'Exposure time (m)', fontsize=12) - plt.ylabel('Probability of infection\n$P(\,\mathtt{I}\,)$', fontsize=12) + plt.xlabel(f'Duration of close-proximity encounter (min)', fontsize=12) + plt.ylabel('Mean cumulative dose\n(infectious virus)', fontsize=12) plt.legend() plt.show()