Changed colormap color scale and plot linewidth
This commit is contained in:
Luis Aleixo
commit
5c07820a9e
3 changed files with 155 additions and 92 deletions
|
|
@ -156,7 +156,7 @@ def exposure_model_from_vl_talking(viral_loads):
|
|||
)
|
||||
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
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(1.0)/4
|
||||
er_means.append(np.mean(emission_rate))
|
||||
er_medians.append(np.median(emission_rate))
|
||||
lower_percentiles.append(np.quantile(emission_rate, 0.01))
|
||||
|
|
@ -214,7 +214,7 @@ def exposure_model_from_vl_talking(viral_loads):
|
|||
def exposure_model_from_vl_talking_cn(viral_loads):
|
||||
|
||||
n_lines = 5
|
||||
cns = np.linspace(0., .2, n_lines)
|
||||
cns = np.linspace(0.1, 1, n_lines)
|
||||
norm = mpl.colors.Normalize(vmin=cns.min(), vmax=cns.max())
|
||||
cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.jet)
|
||||
cmap.set_array([])
|
||||
|
|
@ -264,8 +264,8 @@ def exposure_model_from_vl_talking_cn(viral_loads):
|
|||
|
||||
|
||||
ax.plot(viral_loads, er_means, color=cmap.to_rgba(cn))
|
||||
ax.fill_between(viral_loads, lower_percentiles,
|
||||
upper_percentiles, alpha=0.2, color=cmap.to_rgba(cn))
|
||||
#ax.fill_between(viral_loads, lower_percentiles,
|
||||
# upper_percentiles, alpha=0.2, color=cmap.to_rgba(cn))
|
||||
|
||||
fig.colorbar(cmap, ticks=cns)
|
||||
ax.set_yscale('log')
|
||||
|
|
@ -341,7 +341,7 @@ def exposure_model_from_vl_breathing(viral_loads):
|
|||
)
|
||||
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
|
||||
emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(1.0) / 2
|
||||
er_means.append(np.mean(emission_rate))
|
||||
er_medians.append(np.median(emission_rate))
|
||||
lower_percentiles.append(np.quantile(emission_rate, 0.01))
|
||||
|
|
@ -436,6 +436,144 @@ def exposure_model_from_vl_breathing(viral_loads):
|
|||
return er_means
|
||||
|
||||
|
||||
def exposure_model_from_vl_breathing_cn(viral_loads):
|
||||
n_lines = 5
|
||||
cns = np.linspace(0.01, 0.5, n_lines)
|
||||
norm = mpl.colors.Normalize(vmin=cns.min(), vmax=cns.max())
|
||||
cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.gray)
|
||||
cmap.set_array([])
|
||||
|
||||
for cn in tqdm(cns):
|
||||
er_means = []
|
||||
er_medians = []
|
||||
lower_percentiles = []
|
||||
upper_percentiles = []
|
||||
for vl in viral_loads:
|
||||
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=models.Virus(
|
||||
viral_load_in_sputum=10**vl,
|
||||
infectious_dose=50.,
|
||||
),
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
mask=models.Mask.types["No mask"],
|
||||
activity=activity_distributions['Seated'],
|
||||
expiration=models.Expiration.types['Breathing'],
|
||||
),
|
||||
),
|
||||
exposed=mc.Population(
|
||||
number=14,
|
||||
presence=mc.SpecificInterval(((0, 2),)),
|
||||
activity=models.Activity.types['Seated'],
|
||||
mask=models.Mask.types["No mask"],
|
||||
),
|
||||
)
|
||||
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 = 1.0) / 2
|
||||
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))
|
||||
|
||||
# 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]
|
||||
ax.plot(viral_loads, er_means, color=cmap.to_rgba(cn), linewidth=1)
|
||||
|
||||
#ax.fill_between(viral_loads, lower_percentiles,
|
||||
# upper_percentiles, alpha=0.2)
|
||||
|
||||
fig.colorbar(cmap, ticks=cns)
|
||||
ax.set_yscale('log')
|
||||
|
||||
############# Coleman #############
|
||||
plt.scatter(coleman_etal_vl_breathing,
|
||||
coleman_etal_er_breathing_2, marker='x')
|
||||
x_hull, y_hull = get_enclosure_points(
|
||||
coleman_etal_vl_breathing, coleman_etal_er_breathing_2)
|
||||
# plot shape
|
||||
plt.fill(x_hull, y_hull, '--', c='orange', alpha=0.2)
|
||||
|
||||
############# Markers #############
|
||||
markers = [5, 'd', 4]
|
||||
|
||||
############# Milton et al #############
|
||||
try:
|
||||
for index, m in enumerate(markers):
|
||||
plt.scatter(milton_vl[index], milton_er_2[index],
|
||||
marker=m, color='red')
|
||||
x_hull, y_hull = get_enclosure_points(milton_vl, milton_er_2)
|
||||
# plot shape
|
||||
plt.fill(x_hull, y_hull, '--', c='red', alpha=0.2)
|
||||
except:
|
||||
print("No data for Milton et al")
|
||||
|
||||
############# Yan et al #############
|
||||
try:
|
||||
plt.scatter(yann_vl[0], yann_er_2[0], marker=markers[0], color='green')
|
||||
plt.scatter(yann_vl[1], yann_er_2[1],
|
||||
marker=markers[1], color='green', s=50)
|
||||
plt.scatter(yann_vl[2], yann_er_2[2], marker=markers[2], color='green')
|
||||
|
||||
x_hull, y_hull = get_enclosure_points(yann_vl, yann_er_2)
|
||||
# plot shape
|
||||
plt.fill(x_hull, y_hull, '--', c='green', alpha=0.2)
|
||||
except:
|
||||
print("No data for Yan et al")
|
||||
|
||||
############ Legend ############
|
||||
result_from_model = mlines.Line2D(
|
||||
[], [], color='blue', marker='_', linestyle='None')
|
||||
coleman = mlines.Line2D([], [], color='orange',
|
||||
marker='x', linestyle='None')
|
||||
milton_mean = mlines.Line2D(
|
||||
[], [], color='red', marker='d', linestyle='None') # mean
|
||||
milton_25 = mlines.Line2D(
|
||||
[], [], color='red', marker=5, linestyle='None') # 25
|
||||
milton_75 = mlines.Line2D(
|
||||
[], [], color='red', marker=4, linestyle='None') # 75
|
||||
yann_mean = mlines.Line2D([], [], color='green',
|
||||
marker='d', linestyle='None') # mean
|
||||
yann_25 = mlines.Line2D([], [], color='green',
|
||||
marker=5, linestyle='None') # 25
|
||||
yann_75 = mlines.Line2D([], [], color='green',
|
||||
marker=4, linestyle='None') # 75
|
||||
|
||||
title_proxy = Rectangle((0, 0), 0, 0, color='w')
|
||||
titles = ["$\\bf{CARA \, \\it{(SARS-CoV-2)}:}$", "$\\bf{Coleman \, et \, al. \, \\it{(SARS-CoV-2)}:}$",
|
||||
"$\\bf{Milton \, et \, al. \,\\it{(Influenza)}:}$", "$\\bf{Yann \, et \, al. \,\\it{(Influenza)}:}$"]
|
||||
leg = plt.legend([title_proxy, result_from_model, title_proxy, coleman, title_proxy, milton_mean, milton_25, milton_75, title_proxy, yann_mean, yann_25, yann_75],
|
||||
[titles[0], "Result from model", titles[1], "Dataset", titles[2], "Mean", "25th per.", "75th per.", titles[3], "Mean", "25th per.", "75th per."])
|
||||
|
||||
# Move titles to the left
|
||||
for item, label in zip(leg.legendHandles, leg.texts):
|
||||
if label._text in titles:
|
||||
width = item.get_window_extent(fig.canvas.get_renderer()).width
|
||||
label.set_ha('left')
|
||||
label.set_position((-3*width, 0))
|
||||
|
||||
############ Plot ############
|
||||
plt.title('Exhaled virions while breathing for 30 min',
|
||||
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()
|
||||
|
||||
return er_means
|
||||
|
||||
|
||||
######### Auxiliar functions #########
|
||||
|
||||
def get_enclosure_points(x_coordinates, y_coordinates):
|
||||
|
|
|
|||
|
|
@ -555,19 +555,19 @@ class Expiration(_ExpirationBase):
|
|||
BLO_factors: typing.Tuple[float, float, float]
|
||||
|
||||
@cached()
|
||||
def aerosols(self, mask: Mask, cn: float):
|
||||
def aerosols(self, mask: Mask, cn_B: float, cn_L: float):
|
||||
""" Result is in mL.cm^-3 """
|
||||
def volume(d):
|
||||
return (np.pi * d**3) / 6.
|
||||
|
||||
def _Bmode(d: float) -> float:
|
||||
def _Bmode(d: float, cn_B: float) -> float:
|
||||
# B-mode (see ref. above).
|
||||
return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) *
|
||||
return ( (1 / d) * (cn_B / (np.sqrt(2 * np.pi) * 0.262364)) *
|
||||
np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2)))
|
||||
|
||||
def _Lmode(d: float, cn: float) -> float:
|
||||
def _Lmode(d: float, cn_L: float) -> float:
|
||||
# L-mode (see ref. above).
|
||||
return ( (1 / d) * (cn / (np.sqrt(2 * np.pi) * 0.506818)) *
|
||||
return ( (1 / d) * (cn_L / (np.sqrt(2 * np.pi) * 0.506818)) *
|
||||
np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2)))
|
||||
|
||||
def _Omode(d: float) -> float:
|
||||
|
|
@ -576,8 +576,8 @@ class Expiration(_ExpirationBase):
|
|||
np.exp(-1 * (np.log(d) - 4.97673) ** 2 / (2 * 0.585005 ** 2)))
|
||||
|
||||
def integrand(d: float) -> float:
|
||||
return (self.BLO_factors[0] * _Bmode(d) +
|
||||
self.BLO_factors[1] * _Lmode(d, cn) +
|
||||
return (self.BLO_factors[0] * _Bmode(d, cn_B) +
|
||||
self.BLO_factors[1] * _Lmode(d, cn_L) +
|
||||
self.BLO_factors[2] * _Omode(d)
|
||||
) * volume(d) * (1 - mask.exhale_efficiency(d))
|
||||
|
||||
|
|
@ -670,7 +670,7 @@ class InfectedPopulation(Population):
|
|||
#: The type of expiration that is being emitted whilst doing the activity.
|
||||
expiration: _ExpirationBase
|
||||
|
||||
def emission_rate_when_present(self, cn: float) -> _VectorisedFloat:
|
||||
def emission_rate_when_present(self, cn_B: float, cn_L: float) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate if the infected population is present.
|
||||
|
||||
|
|
@ -680,7 +680,7 @@ class InfectedPopulation(Population):
|
|||
# Emission Rate (virions / h)
|
||||
# Note on units: exhalation rate is in m^3/h, aerosols in mL/cm^3
|
||||
# and viral load in virus/mL -> 1e6 conversion factor
|
||||
aerosols = self.expiration.aerosols(self.mask, cn)
|
||||
aerosols = self.expiration.aerosols(self.mask, cn_B, cn_L)
|
||||
|
||||
ER = (self.virus.viral_load_in_sputum *
|
||||
self.activity.exhalation_rate *
|
||||
|
|
|
|||
|
|
@ -1,13 +1,14 @@
|
|||
from cara.model_scenarios import exposure_model_from_vl_breathing, exposure_model_from_vl_talking, exposure_model_from_vl_talking_cn, exposure_model_from_vl_talking_new_points
|
||||
from cara.model_scenarios import *
|
||||
import numpy as np
|
||||
import csv
|
||||
|
||||
viral_loads = np.linspace(2, 12, 600)
|
||||
|
||||
#er_means = exposure_model_from_vl_talking(viral_loads)
|
||||
#er_means = exposure_model_from_vl_talking_new_points(viral_loads)
|
||||
#er_means = exposure_model_from_vl_breathing(viral_loads)
|
||||
#er_means = exposure_model_from_vl_talking_new_points(viral_loads)
|
||||
#er_means = exposure_model_from_vl_talking_cn(viral_loads)
|
||||
er_means = exposure_model_from_vl_breathing_cn(viral_loads)
|
||||
|
||||
# with open('data.csv', 'w', newline='') as csvfile:
|
||||
# fieldnames = ['viral load', 'emission rate']
|
||||
|
|
@ -16,79 +17,3 @@ viral_loads = np.linspace(2, 12, 600)
|
|||
# for i, vl in enumerate(viral_loads):
|
||||
# thewriter.writerow(
|
||||
# {'viral load': 10**vl, 'emission rate': er_means[i]})
|
||||
|
||||
|
||||
def fit_function_to_data_points():
|
||||
|
||||
rna_copies = np.array([4.01624,
|
||||
4.38393,
|
||||
4.65486,
|
||||
4.99213,
|
||||
5.35982,
|
||||
5.66392,
|
||||
5.90444,
|
||||
6.11178,
|
||||
6.30254,
|
||||
6.47947,
|
||||
6.67023,
|
||||
6.83057,
|
||||
6.97433,
|
||||
7.13467,
|
||||
7.24802,
|
||||
7.4056,
|
||||
7.59912,
|
||||
7.80646,
|
||||
7.9834,
|
||||
8.11057,
|
||||
8.23774,
|
||||
8.41467,
|
||||
8.55843,
|
||||
8.74918,
|
||||
8.97311,
|
||||
9.23022,
|
||||
9.43756,
|
||||
9.74166,
|
||||
10.06235,
|
||||
10.34987,
|
||||
10.59038,
|
||||
10.76455,
|
||||
10.92489])
|
||||
|
||||
r_inf = np.array([0.004036804,
|
||||
0.003189439,
|
||||
0.003189439,
|
||||
0.007288068,
|
||||
0.013790595,
|
||||
0.022835218,
|
||||
0.03258901,
|
||||
0.043190166,
|
||||
0.05392952,
|
||||
0.066083573,
|
||||
0.080077827,
|
||||
0.093079245,
|
||||
0.108626396,
|
||||
0.121773284,
|
||||
0.135622068,
|
||||
0.149616322,
|
||||
0.171666,
|
||||
0.192864676,
|
||||
0.212510456,
|
||||
0.228057606,
|
||||
0.238658763,
|
||||
0.254205913,
|
||||
0.268905699,
|
||||
0.284452849,
|
||||
0.3,
|
||||
0.315547151,
|
||||
0.326995672,
|
||||
0.338444194,
|
||||
0.346641452,
|
||||
0.353143979,
|
||||
0.356391606,
|
||||
0.357242608,
|
||||
0.357242608])
|
||||
|
||||
result = np.polyfit(rna_copies, r_inf, 1)
|
||||
print(result)
|
||||
|
||||
fit_function_to_data_points()
|
||||
|
|
|
|||
Loading…
Reference in a new issue