pingu_98/cara
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge remote-tracking branch 'origin/feature/mc' into feature/mc

Browse source
This commit is contained in:
markus 2021-03-22 17:13:33 +01:00
commit 741baa04be
6 changed files with 92 additions and 37 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

9
cara/mc-output-publication.py Normal file
View file

@ -0,0 +1,9 @@
""" Title: COVID Airborne Risk Assessment
Author: <author(s) names>
Date: <date>
Code version: <code version>
Availability: <where it's located> """
from cara.montecarlo import *
from cara.model_scenarios import *

26
cara/mc-output.py
View file

@ -7,7 +7,7 @@ from cara.model_scenarios import *
#print(np.mean(classroom_model[0].infection_probability()))
#print(np.mean(classroom_model[1].infection_probability()))
#print(np.mean(chorale_model.infection_probability())+np.std(chorale_model.infection_probability()))
#print(np.quantile(chorale_model.infection_probability(),0.8))
#print(np.quantile(chorale_model[0].infection_probability(),0.8))
#print(np.quantile(chorale_model.infection_probability(),0.90))
#print(np.quantile(chorale_model.infection_probability(),0.1))
@ -20,20 +20,20 @@ from cara.model_scenarios import *
# time_in_minutes=True,
# normalize_y_axis=True)
compare_viruses_qr(violins=True)
#compare_viruses_qr(violins=True)
# print_qd_info(large_population_baselines[0])
#print_qd_info(waiting_room_model_full_winter[1])
#print(np.mean(shared_office_model[1].infection_probability()))
#composite_plot_pi_vs_viral_load([shared_office_worst_model[1], shared_office_model[1], shared_office_better_model[1]],
# labels=['No mask &\nwindows closed', 'Baseline:Windows\nopen (10min/2h)', 'Windows open\n(10min/2h) +\nHEPA filter'],
# colors=['tomato', '#1f77b4', 'limegreen'],
# title='Shared office scenario',
# title='',
# vl_points=200)
#composite_plot_pi_vs_viral_load([classroom_model_no_vent[1], classroom_model[1], classroom_model_with_hepa[1], classroom_model_full_open_multi[1], classroom_model_full_open_multi_masks[1]],
# labels=['Windows closed', 'Baseline:Windows\nopen (10min/2h)', 'Windows open\n(10min/2h) + HEPA', 'Multiple windows open', 'Multiple windows open\n+masks'],
# colors=['tomato','#1f77b4', 'dodgerblue', 'seagreen', 'limegreen'],
# title='Classroom scenario',
# title='',
# vl_points=200)
#composite_plot_pi_vs_viral_load([ski_cabin_model_60[1], ski_cabin_model_30[1], ski_cabin_model_baseline_20[1], ski_cabin_model_10[1]],
@ -49,16 +49,16 @@ compare_viruses_qr(violins=True)
# vl_points=200)
#compare_concentration_curves([classroom_model_no_vent[1], classroom_model[1], classroom_model_with_hepa[1], classroom_model_full_open_multi[1]],
# labels=['Windows closed', 'Baseline:(windows 10min/2h)', 'Baseline:(windows 10min/2h) + HEPA', 'Multiple windows open'],
# labels=['Windows closed', 'Baseline: Windows open (10min/2h)', 'Windows open (10min/2h) + HEPA', 'Multiple windows open'],
# colors=['tomato','#1f77b4', 'seagreen', 'limegreen'],
# title='Classroom scenario'
# title=''
# )
#
#compare_concentration_curves([waiting_room_model[1], waiting_room_model_full_summer[1],
# waiting_room_model_full_winter[1], waiting_room_model_periodic_winter[1]],
# labels=['Baseline:(windows closed)', 'Windows open (summer)', 'Windows open (winter)', 'Windows open 5min/20min (winter)'],
# labels=['Baseline: Windows closed', 'Windows open (summer)', 'Windows open (winter)', 'Windows open 5min/20min (winter)'],
# colors=['#1f77b4', 'darkorange', 'deepskyblue', 'lightskyblue'],
# title='Waiting room scenario'
# title=''
# )
#plot_pi_vs_viral_load([shared_office_model[1]], labels=['Baseline'],title='')
@ -91,3 +91,9 @@ compare_viruses_qr(violins=True)
# plt.title(f'Probability of infection in baseline case - {"English" if model.concentration_model.infected.qid == 60 else "Original"} variant')
# plt.yticks([], [])
# plt.show()
#composite_plot_pi_vs_viral_load([shared_office_worst_model[1], shared_office_worst_model_infiltration[1]],
# labels=['No mask &\nwindows closed', 'Baseline:Windows\nopen (10min/2h)', 'Windows open\n(10min/2h) +\nHEPA filter'],
# colors=['tomato', '#1f77b4', 'limegreen'],
# title='Shared office scenario',
# vl_points=200)

28
cara/model_scenarios.py
View file

@ -129,7 +129,7 @@ classroom_model = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -156,7 +156,7 @@ classroom_model_full_open = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -183,7 +183,7 @@ classroom_model_summer = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -210,7 +210,7 @@ classroom_model_full_open_summer = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -227,17 +227,15 @@ classroom_model_full_open_summer = [MCExposureModel(
classroom_model_no_vent = [MCExposureModel(
concentration_model=MCConcentrationModel(
room=models.Room(volume=160),
ventilation=models.SlidingWindow(
ventilation=models.AirChange(
active=models.PeriodicInterval(period=120, duration=120),
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
outside_temp=models.PiecewiseConstant((0, 24), (291,)),
window_height=1.6, opening_length=0.,
air_exch=0.25,
),
infected=MCInfectedPopulation(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -264,7 +262,7 @@ classroom_model_teacher_mask_full_open = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=True,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -297,7 +295,7 @@ classroom_model_with_hepa = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -330,7 +328,7 @@ classroom_model_full_open_multi = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -363,7 +361,7 @@ classroom_model_full_open_multi_masks = [MCExposureModel(
number=1,
presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))),
masked=True,
virus=MCVirus(halflife=1.1, qID=qid),
virus=MCVirus(halflife=3.8, qID=qid),
expiratory_activity=2,
samples=200000,
breathing_category=3,
@ -439,7 +437,7 @@ shared_office_worst_model = [MCExposureModel(
room=models.Room(volume=50),
ventilation=models.AirChange(
active=models.PeriodicInterval(period=120, duration=120),
air_exch=0.,
air_exch=0.25,
),
infected=MCInfectedPopulation(
number=1,
@ -788,7 +786,7 @@ waiting_room_model = [MCExposureModel(
room=models.Room(volume=100),
ventilation=models.AirChange(
active=models.PeriodicInterval(period=120, duration=120),
air_exch=0.,
air_exch=0.25,
),
infected=MCInfectedPopulation(
number=1,

31
cara/model_scenarios_publication.py Normal file
View file

@ -0,0 +1,31 @@
from cara import models
from cara.montecarlo import *
######### Standard exposure models ###########
exposure_models = [MCExposureModel(
concentration_model=MCConcentrationModel(
room=models.Room(volume=45),
ventilation=models.SlidingWindow(
active=models.PeriodicInterval(period=120, duration=120),
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
outside_temp=models.PiecewiseConstant((0, 24), (283,)),
window_height=1.6, opening_length=0.2,
),
infected=MCInfectedPopulation(
number=1,
presence=models.SpecificInterval(((0, 4), (5, 9))),
masked=False,
virus=MCVirus(halflife=1.1, qID=100),
expiratory_activity=a,
samples=2000000,
breathing_category=3,
expiratory_activity_weights=(0.7, 0.3, 0)
)
),
exposed=models.Population(
number=2,
presence=models.SpecificInterval(((0, 4), (5, 9))),
activity=models.Activity.types['Seated'],
mask=models.Mask.types['Type I']
)
) for a in (1, 2, 3)]

2
cara/models.py
View file

@ -236,7 +236,7 @@ class WindowOpening(Ventilation):
def air_exchange(self, room: Room, time: float) -> float:
# If the window is shut, no air is being exchanged.
if not self.active.triggered(time):
return 0.
return 0.25
# Reminder, no dependence on time in the resulting calculation.
inside_temp = self.inside_temp.value(time)

33
cara/montecarlo.py
View file

@ -11,7 +11,7 @@ import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.lines as mlines
from sklearn.neighbors import KernelDensity
TIME_STEP = 0.01
TIME_STEP = 0.001
USE_SCOEH = False
@ -472,11 +472,22 @@ def print_qr_info(log_qr: np.ndarray) -> None:
print(f"MEAN of log_10(qR) = {np.mean(log_qr)}\n"
f"SD of log_10(qR) = {np.std(log_qr)}\n"
f"MEAN of qR = {np.mean(qr_values)}\n"
f"SD of qR = {np.std(qr_values)}\n")
f"SD of qR = {np.std(qr_values)}\n"
f"Median of qR = {np.quantile(qr_values, 0.5)}\n")
for quantile in (0.01, 0.05, 0.25, 0.50, 0.75, 0.95, 0.99):
print(f"qR_{quantile} = {np.quantile(qr_values, quantile)}")
def present_qR_values(model: MCConcentrationModel) -> None:
"""
Displays a handful of key parameters and results of qR from a given MCConcentrationModel
:param model: The MCConcentrationModel representing the scenario to be presented
:return: Nothing, parameters are printed
"""
qRs = (np.log10(model.infected.emission_rate_when_present()))
print_qr_info(np.asarray(qRs))
def present_model(model: MCConcentrationModel, bins: int = 200,
title: str = 'Summary of $qR$ model parameters') -> None:
@ -734,7 +745,7 @@ def composite_plot_pi_vs_viral_load(baselines: typing.List[MCExposureModel], lab
axs[1, 0].set_xticklabels(['$10^{' + str(i) + '}$' for i in range(2, 13, 2)])
axs[1, 0].set_xlim(2, 12)
axs[1, 0].set_xlabel('Viral load (RNA copies mL$^{-1}$)', fontsize=12)
axs[0, 0].set_ylabel('Probability of infection\n$P(I|qID=60)$', fontsize=12)
axs[0, 0].set_ylabel('Probability of infection\n$P(\,\mathtt{I}\,|\,\mathrm{vl}\,)$', fontsize=12)
plt.suptitle(title, fontsize=12)
axs[0, 0].text(11, -0.01, '$(i)$')
@ -873,8 +884,8 @@ def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid:
"""
fig, axs = plt.subplots(3, 1, sharex='all')
plt.suptitle("$F(qR|qID=$" + str(qid) + "$)$",fontsize=14, y=0.93)
fig.text(0.02, 0.5, 'Cumulative Distribution Function', va='center', rotation='vertical',fontsize=14)
plt.suptitle("$F(\mathrm{qR}|\mathrm{qID}=$" + str(qid) + "$)$",fontsize=16, y=0.93)
fig.text(0.02, 0.5, '', va='center', rotation='vertical',fontsize=16)
scenarios = [MCInfectedPopulation(
number=1,
@ -899,17 +910,17 @@ def generate_cdf_curves_vs_qr(masked: bool = False, samples: int = 200000, qid:
for i in range(3):
axs[i].hist(qr_values[3 * i:3 * (i + 1)], bins=2000, histtype='step',
color=colors, cumulative=True, range=(left, right))
axs[i].set_xlim(left, right)
color=colors, cumulative=True, range=(-7, 6))
axs[i].set_xlim(-6, 6)
axs[i].set_yticks([0, samples / 2, samples])
axs[i].set_yticklabels(['0.0', '0.5', '1.0'])
axs[i].yaxis.set_label_position("right")
axs[i].set_ylabel(activities[i],fontsize=12)
axs[i].set_ylabel(activities[i], fontsize=14)
axs[i].grid(linestyle='--')
axs[0].legend(handles=lines, loc='upper left')
plt.xlabel('qR (q h$^{-1}$)', fontsize=12)
tick_positions = np.arange(int(np.ceil(left)), int(np.ceil(right)), 2)
plt.xlabel('$\mathrm{qR}$ (q h$^{-1}$)', fontsize=16)
tick_positions = np.arange(-6, 6, 2)
plt.xticks(ticks=tick_positions, labels=['$\;10^{' + str(i) + '}$' for i in tick_positions])
fig.set_figheight(8)
@ -1304,6 +1315,6 @@ def plot_pi_vs_exposure_time(exp_models: typing.List[MCExposureModel], labels: t
plt.title('')
plt.xlabel(f'Travel time ({"min" if time_in_minutes else "h"})', fontsize=12)
plt.ylabel('Probability of infection\n$P(I|qID=60)$', fontsize=12)
plt.ylabel('Probability of infection\n$P(\,\mathtt{I}\,|\,\mathrm{vl}\,)$', fontsize=12)
plt.legend()
plt.show()