From cbc5ffb203c283369360806997f6a130f2b2a4e9 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 18 Aug 2021 14:23:46 +0200 Subject: [PATCH] added methods to different scenarios --- cara/model_scenarios.py | 267 ++++++++++++++++++++++++++++++++++++++++ cara/plot_output.py | 192 +---------------------------- 2 files changed, 271 insertions(+), 188 deletions(-) create mode 100644 cara/model_scenarios.py diff --git a/cara/model_scenarios.py b/cara/model_scenarios.py new file mode 100644 index 00000000..c05a91cb --- /dev/null +++ b/cara/model_scenarios.py @@ -0,0 +1,267 @@ +from cara import models +from cara.monte_carlo.data import activity_distributions +from tqdm import tqdm +import cara.monte_carlo as mc +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial import ConvexHull +import pandas as pd +import matplotlib.lines as mlines +from matplotlib.patches import Rectangle + +######### Plot material ######### +fig = plt.figure() +ax = fig.add_subplot(1, 1, 1) + +SAMPLE_SIZE = 50000 + +er_means = [] +er_medians = [] +lower_percentiles = [] +upper_percentiles = [] + +######### Scatter points ######### +############# Coleman ############# + +############# Coleman - Breathing ############# +coleman_etal_vl_breathing = [np.log10(821065925.4), np.log10(1382131207), np.log10(81801735.96), np.log10( + 487760677.4), np.log10(2326593535), np.log10(1488879159), np.log10(884480386.5)] +coleman_etal_er_breathing = [127, 455.2, 281.8, 884.2, 448.4, 1100.6, 621] +############# Coleman - Talking ############# +coleman_etal_vl_talking = [np.log10(70492378.55), np.log10(7565486.029), np.log10(7101877592), np.log10(1382131207), np.log10(821065925.4), + np.log10(1382131207), np.log10(81801735.96), np.log10( + 487760677.4), np.log10(2326593535), np.log10(1488879159), + np.log10(884480386.5)] +coleman_etal_er_talking = [417, 234.5, 79.9, 908.2, 310.9, + 4336, 733, 1356.5, 1373.3, 477.9, 2428.7] +############# Milton et al ############# +milton_vl = [np.log10(8.30E+04), np.log10(4.20E+05), np.log10(1.80E+06)] +milton_er = [22, 220, 1120] # removed first and last due to its dimensions +############# Milton et al ############# +yann_vl = [np.log10(7.86E+07), np.log10(2.23E+09), np.log10(1.51E+10)] +yann_er = [8396.78166, 45324.55964, 400054.0827] + +######### Standard exposure models ########### + + +def exposure_model_from_vl_talking(viral_loads): + + for vl in tqdm(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['Talking'], + ), + ), + 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) + emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present() + 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)) + + build_plot(viral_loads, er_means, + lower_percentiles, upper_percentiles, coleman_etal_vl_talking, coleman_etal_er_talking, milton_vl, milton_er, yann_vl, yann_er) + + ############ Plot ############ + plt.title('Exhaled virions while talking for 1h', + 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 + + +def exposure_model_from_vl_breathing(viral_loads): + + for vl in tqdm(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) + emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present() + 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)) + + build_plot(viral_loads, er_means, + lower_percentiles, upper_percentiles, coleman_etal_vl_breathing, coleman_etal_er_breathing, milton_vl, milton_er, yann_vl, yann_er) + + ############ Plot ############ + plt.title('Exhaled virions while breathing for 1h', + 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): + df = pd.DataFrame({'x': x_coordinates, 'y': y_coordinates}) + + points = df[['x', 'y']].values + # get convex hull + hull = ConvexHull(points) + # get x and y coordinates + # repeat last point to close the polygon + x_hull = np.append(points[hull.vertices, 0], + points[hull.vertices, 0][0]) + y_hull = np.append(points[hull.vertices, 1], + points[hull.vertices, 1][0]) + return x_hull, y_hull + + +def build_legend(fig): + 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)) + + +def build_plot(viral_loads, er_means, lower_percentiles, upper_percentiles, coleman_etal_vl, coleman_etal_er, milton_vl, milton_er, yann_vl, yann_er, ): + ax.plot(viral_loads, er_means) + ax.fill_between(viral_loads, lower_percentiles, + upper_percentiles, alpha=0.2) + ax.set_yscale('log') + + ############# Coleman ############# + plt.scatter(coleman_etal_vl, coleman_etal_er, marker='x') + x_hull, y_hull = get_enclosure_points( + coleman_etal_vl, coleman_etal_er) + # plot shape + plt.fill(x_hull, y_hull, '--', c='orange', alpha=0.2) + + ############# Markers ############# + markers = [5, 'd', 4] + + ############# Milton et al ############# + for index, m in enumerate(markers): + plt.scatter(milton_vl[index], milton_er[index], + marker=m, color='red') + x_hull, y_hull = get_enclosure_points(milton_vl, milton_er) + # plot shape + plt.fill(x_hull, y_hull, '--', c='red', alpha=0.2) + + ############# Yan et al ############# + plt.scatter(yann_vl[0], yann_er[0], marker=markers[0], color='green') + plt.scatter(yann_vl[1], yann_er[1], marker=markers[1], color='green', s=50) + plt.scatter(yann_vl[2], yann_er[2], marker=markers[2], color='green') + + x_hull, y_hull = get_enclosure_points(yann_vl, yann_er) + # plot shape + plt.fill(x_hull, y_hull, '--', c='green', alpha=0.2) + + ############ Legend ############ + build_legend(fig) + + return plt + + +# # Milton + # boxes = [ + # { + # 'label': "Milton data", + # 'whislo': 0, # Bottom whisker position + # 'q1': 22, # First quartile (25th percentile) + # 'med': 220, # Median (50th percentile) + # 'q3': 1120, # Third quartile (75th percentile) + # 'whishi': 260000, # Top whisker position + # 'fliers': [] # Outliers + # } + # ] + # # `box plot aligned with the viral load value of 5.62325 + # ax.bxp(boxes, showfliers=False, positions=[5.62324929]) + + # # Yan + + # boxes = [ + # { + # 'whislo': 1424.81, # Bottom whisker position + # 'q1': 8396.78, # First quartile (25th percentile) + # 'med': 45324.6, # Median (50th percentile) + # 'q3': 400054, # Third quartile (75th percentile) + # 'whishi': 88616200, # Top whisker position + # 'fliers': [] # Outliers + # } + # ] + # ax.bxp(boxes, showfliers=False, positions=[9.34786]) + # box plot aligned with the viral load value of 9.34786 diff --git a/cara/plot_output.py b/cara/plot_output.py index a39140d4..e7f58869 100644 --- a/cara/plot_output.py +++ b/cara/plot_output.py @@ -1,80 +1,11 @@ -from dataclasses import field +from cara.model_scenarios import exposure_model_from_vl_breathing, exposure_model_from_vl_talking import numpy as np -import matplotlib.pyplot as plt -import matplotlib.lines as mlines -from matplotlib.patches import Rectangle -import pandas as pd import csv -import cara.monte_carlo as mc -from cara import models -from cara.monte_carlo.data import activity_distributions -from tqdm import tqdm -from scipy.spatial import ConvexHull +viral_loads = np.linspace(2, 12, 600) - -def get_enclosure_points(x_coordinates, y_coordinates): - df = pd.DataFrame({'x': x_coordinates, 'y': y_coordinates}) - - points = df[['x', 'y']].values - # get convex hull - hull = ConvexHull(points) - # get x and y coordinates - # repeat last point to close the polygon - x_hull = np.append(points[hull.vertices, 0], - points[hull.vertices, 0][0]) - y_hull = np.append(points[hull.vertices, 1], - points[hull.vertices, 1][0]) - return x_hull, y_hull - - -SAMPLE_SIZE = 50000 - -fig = plt.figure() -ax = fig.add_subplot(1, 1, 1) - -points = 600 -viral_loads = np.linspace(2, 12, points) - -er_means = [] -er_medians = [] -lower_percentiles = [] -upper_percentiles = [] - -for vl in tqdm(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'], - expiration=models.Expiration.types['Talking'], - ), - ), - 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) - emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present() - 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)) +#er_means = exposure_model_from_vl_talking(viral_loads) +er_means = exposure_model_from_vl_breathing(viral_loads) with open('data.csv', 'w', newline='') as csvfile: fieldnames = ['viral load', 'emission rate'] @@ -83,118 +14,3 @@ with open('data.csv', 'w', newline='') as csvfile: for i, vl in enumerate(viral_loads): thewriter.writerow( {'viral load': 10**vl, 'emission rate': er_means[i]}) - -ax.plot(viral_loads, er_means) -ax.fill_between(viral_loads, lower_percentiles, - upper_percentiles, alpha=0.2) -ax.set_yscale('log') - -############# Coleman ############# -#Breathing -coleman_etal_vl = [np.log10(821065925.4), np.log10(1382131207), np.log10(81801735.96), np.log10( - 487760677.4), np.log10(2326593535), np.log10(1488879159), np.log10(884480386.5)] -coleman_etal_er = [127, 455.2, 281.8, 884.2, 448.4, 1100.6, 621] -#Talking -coleman_etal_vl = [np.log10(70492378.55), np.log10(7565486.029), np.log10(7101877592), np.log10(1382131207), np.log10(821065925.4), - np.log10(1382131207), np.log10(81801735.96), np.log10(487760677.4), np.log10(2326593535), np.log10(1488879159), - np.log10(884480386.5)] -coleman_etal_er = [417, 234.5, 79.9, 908.2, 310.9, 4336, 733, 1356.5, 1373.3, 477.9, 2428.7] - -plt.scatter(coleman_etal_vl, coleman_etal_er, marker='x') -x_hull, y_hull = get_enclosure_points(coleman_etal_vl, coleman_etal_er) -# plot shape -plt.fill(x_hull, y_hull, '--', c='orange', alpha=0.2) - -############# Markers ############# -markers = [5, 'd', 4] - -############# Milton et al ############# -milton_vl = [np.log10(8.30E+04), np.log10(4.20E+05), np.log10(1.80E+06)] -milton_er = [22, 220, 1120] # removed first and last due to its dimensions -plt.scatter(milton_vl[0], milton_er[0], marker=markers[0], color='red') -plt.scatter(milton_vl[1], milton_er[1], marker=markers[1], color='red', s=50) -plt.scatter(milton_vl[2], milton_er[2], marker=markers[2], color='red') -x_hull, y_hull = get_enclosure_points(milton_vl, milton_er) -# plot shape -plt.fill(x_hull, y_hull, '--', c='red', alpha=0.2) - -############# Yan et al ############# -# removed first and last due to its dimensions -yan_vl = [np.log10(7.86E+07), np.log10(2.23E+09), np.log10(1.51E+10)] -yan_er = [8396.78166, 45324.55964, 400054.0827] -plt.scatter(yan_vl[0], yan_er[0], marker=markers[0], color='green') -plt.scatter(yan_vl[1], yan_er[1], marker=markers[1], color='green', s=50) -plt.scatter(yan_vl[2], yan_er[2], marker=markers[2], color='green') - -x_hull, y_hull = get_enclosure_points(yan_vl, yan_er) -# plot shape -plt.fill(x_hull, y_hull, '--', c='green', alpha=0.2) - -# # Milton -# boxes = [ -# { -# 'label': "Milton data", -# 'whislo': 0, # Bottom whisker position -# 'q1': 22, # First quartile (25th percentile) -# 'med': 220, # Median (50th percentile) -# 'q3': 1120, # Third quartile (75th percentile) -# 'whishi': 260000, # Top whisker position -# 'fliers': [] # Outliers -# } -# ] -# # `box plot aligned with the viral load value of 5.62325 -# ax.bxp(boxes, showfliers=False, positions=[5.62324929]) - -# # Yan - -# boxes = [ -# { -# 'whislo': 1424.81, # Bottom whisker position -# 'q1': 8396.78, # First quartile (25th percentile) -# 'med': 45324.6, # Median (50th percentile) -# 'q3': 400054, # Third quartile (75th percentile) -# 'whishi': 88616200, # Top whisker position -# 'fliers': [] # Outliers -# } -# ] -# ax.bxp(boxes, showfliers=False, positions=[9.34786]) -# box plot aligned with the viral load value of 9.34786 - -############ 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 1h', 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()