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 
# Conflicts:
#	cara/mc-output.py
This commit is contained in:
Andrejh 2021-02-24 16:28:59 +01:00
commit 25b5a16668
2 changed files with 35 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

19
cara/mc-output.py
View file

@ -1,7 +1,8 @@
from cara.montecarlo import * from cara.montecarlo import *
from cara.model_scenarios import * from cara.model_scenarios import *
# compare_concentration_curves([classroom_model, classroom_model_with_hepa], ['Just window', 'Window and HEPA']) classroom_model.infection_probability()
compare_concentration_curves([classroom_model, classroom_model_with_hepa], ['Just window', 'Window and HEPA'])
#print(np.mean(chorale_model.infection_probability())) #print(np.mean(chorale_model.infection_probability()))
#print(np.mean(chorale_model.infection_probability())+np.std(chorale_model.infection_probability())) #print(np.mean(chorale_model.infection_probability())+np.std(chorale_model.infection_probability()))
@ -11,15 +12,13 @@ from cara.model_scenarios import *
#print(np.mean(classroom_model_with_hepa.infection_probability())) # #print(np.mean(classroom_model_with_hepa.infection_probability()))
# composite_plot_pi_vs_viral_load([shared_office_worst_model[1], shared_office_model[1],shared_office_better_model[1]],
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', 'Baseline +\nHEPA filter'],
labels=['No mask &\nwindows closed', 'Baseline', 'Baseline +\nHEPA filter'], # colors=['tomato', '#1f77b4', 'limegreen'],
colors=['tomato', '#1f77b4', 'limegreen'], # title='$P(I|qID)$ vs $vl$ - Shared office scenario',
title='', # vl_points=200)
#title='$P(I|qID=60)$ vs $vl$ - Shared office scenario', #plot_pi_vs_viral_load([shared_office_model[1]], labels=['Baseline, qID=60', 'HEPA, qID=60', 'No mask + windows closed, qID=60'],title='$P(I|qID)$ - Shared office scenario')
vl_points=200)
#plot_pi_vs_viral_load([shared_office_model[1]], labels=['Baseline'],title='')
#generate_cdf_curves_vs_qr(masked=False,qid=1000) #generate_cdf_curves_vs_qr(masked=False,qid=1000)

33
cara/montecarlo.py
View file

@ -1037,16 +1037,35 @@ def compare_concentration_curves(exp_models: typing.List[MCExposureModel], label
times = np.arange(start, stop, TIME_STEP) times = np.arange(start, stop, TIME_STEP)
concentrations = [[np.mean(model.concentration_model.concentration(t)) for t in times] for model in exp_models] concentrations = [[np.median(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): for c, label, color in zip(concentrations, labels, colors):
plt.plot(times, c, label=label, color=color) ax.plot(times, c, label=label, color=color)
plt.legend() ax.legend(loc='upper left')
ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2)
factors = [0.6 * model.exposed.activity.inhalation_rate * (1 - model.exposed.mask.η_inhale) for model in exp_models]
qds = [[np.trapz(c[:i + 1], times[:i + 1]) * factor for i in range(len(times))]
for c, factor in zip(concentrations, factors)]
plt.suptitle("SUPTITLE HERE")
plt.xlabel("XLABEL HERE")
plt.ylabel("YLABEL HERE")
if show_qd: if show_qd:
plt.twinx()
for c, label, color in zip(concentrations, labels, colors):
plt.plot(times, np.cumsum(c), label='qD - ' + label, color=color, linestyle='dotted')
plt.legend() ax1 = ax.twinx()
for qd, label, color in zip(qds, labels, colors):
ax1.plot(times, qd, label='qD - ' + label, color=color, linestyle='dotted')
ax1.set_ylabel('qD LABEL HERE')
ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.2)
ax2 = ax.twinx()
ax2.spines["right"].set_position(("axes", 1.2))
ax2.set_ylabel('RNA LABEL HERE')
ax2.set_ylim(tuple(y * exp_models[0].concentration_model.virus.qID for y in ax1.get_ylim()))
ax1.legend(loc='upper right')
plt.tight_layout()
plt.show() plt.show()