From 9d0c93947dd7a4588c8e2d71638b1a5a6d221b3b Mon Sep 17 00:00:00 2001 From: Andrejh Date: Thu, 2 Sep 2021 09:50:12 +0200 Subject: [PATCH 1/2] add 1.0 as cn default --- cara/model_scenarios.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/model_scenarios.py b/cara/model_scenarios.py index 3f142db3..b9a70885 100644 --- a/cara/model_scenarios.py +++ b/cara/model_scenarios.py @@ -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)) @@ -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)) From 32f9fa0d5093360e7b102cd21a2ef85cc0e1760b Mon Sep 17 00:00:00 2001 From: Andrejh Date: Thu, 2 Sep 2021 10:19:30 +0200 Subject: [PATCH 2/2] add cn_L, cn_B --- cara/model_scenarios.py | 144 +++++++++++++++++++++++++++++++++++++++- cara/models.py | 18 ++--- cara/plot_output.py | 5 +- 3 files changed, 153 insertions(+), 14 deletions(-) diff --git a/cara/model_scenarios.py b/cara/model_scenarios.py index b9a70885..1312115b 100644 --- a/cara/model_scenarios.py +++ b/cara/model_scenarios.py @@ -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') @@ -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.jet) + 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)) + + #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): diff --git a/cara/models.py b/cara/models.py index 35930a16..af986cc9 100644 --- a/cara/models.py +++ b/cara/models.py @@ -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 * diff --git a/cara/plot_output.py b/cara/plot_output.py index 104a4c5d..79dd9dc3 100644 --- a/cara/plot_output.py +++ b/cara/plot_output.py @@ -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']