Added method to generate the BLO curve

cara/models.py

@@ -1097,7 +1097,7 @@ class ShortRangeModel:
             # Verifies if the given time falls within a short range interaction
             if start < time <= finish:
                 dilution = self.dilutions[index]
-                jet_origin_concentration = concentration_model.infected.expiration.jet_origin_concentration()
+                jet_origin_concentration = self.expirations[index].jet_origin_concentration()
                 # Long range concentration normalized by the virus viral load
                 long_range_normed_concentration = concentration_model.concentration(time) / concentration_model.virus.viral_load_in_sputum
                 

cara/short_range_plots/model_scenarios.py

@@ -0,0 +1,208 @@
+""" Title: CARA - COVID Airborne Risk Assessment
+Author: A. Henriques et al
+Date: 18/02/2021
+Code version: 4.0.0
+"""
+
+from cara import models, data
+from cara.monte_carlo.data import activity_distributions, short_range_expiration_distributions, mask_distributions, virus_distributions, dilution_factor
+import cara.monte_carlo as mc
+import numpy as np
+from cara.monte_carlo.sampleable import CustomKernel
+from cara.monte_carlo.data import BLOmodel
+import typing
+from cara.apps.calculator.model_generator import build_expiration
+
+_VectorisedFloat = typing.Union[float, np.ndarray]
+
+scenario_activity_and_expiration = {
+    'office': (
+        'Seated',
+        # Mostly silent in the office, but 1/3rd of time speaking.
+        {'Speaking': 1, 'Breathing': 2}
+    ),
+    'controlroom-day': (
+        'Seated',
+        # Daytime control room shift, 50% speaking.
+        {'Speaking': 1, 'Breathing': 1}
+    ),
+    'controlroom-night': (
+        'Seated',
+        # Nightshift control room, 10% speaking.
+        {'Speaking': 1, 'Breathing': 9}
+    ),
+    # 'meeting': (
+    #     'Seated',
+    #     # Conversation of N people is approximately 1/N% of the time speaking.
+    #     {'Speaking': 1, 'Breathing': self.total_people - 1}
+    # ),
+    'callcentre': ('Seated', 'Speaking'),
+    'library': ('Seated', 'Breathing'),
+    'training': ('Standing', 'Speaking'),
+    'lab': (
+        'Light activity',
+        #Model 1/2 of time spent speaking in a lab.
+        {'Speaking': 1, 'Breathing': 1}),
+    'workshop': (
+        'Moderate activity',
+        #Model 1/2 of time spent speaking in a workshop.
+        {'Speaking': 1, 'Breathing': 1}),
+    'gym':('Heavy exercise', 'Breathing'),
+    }
+
+######### Standard exposure models ###########
+def exposure_module_without_short_range(activity: str, expiration: str, mask: str):
+    if mask == 'No mask':
+        exposure_mask = models.Mask.types['No mask']
+    else:
+        exposure_mask = mask_distributions[mask]
+
+    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=500/100,
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                virus=virus_distributions['SARS_CoV_2_OMICRON'],
+                presence=mc.SpecificInterval(((8.5, 12.5),(13.5, 17.5),)),
+                mask=exposure_mask,
+                activity=activity_distributions[activity],
+                expiration=build_expiration(expiration),
+                host_immunity=0.,
+            ),
+        ),
+        short_range=mc.ShortRangeModel(
+            presence=[],
+            expirations=[],
+            dilutions=[]
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((8.5, 12.5),(13.5, 17.5),)),
+            activity=activity_distributions[activity],
+            mask=exposure_mask,
+            host_immunity=0.,
+        ),
+    )
+    return exposure_mc
+
+def exposure_module_with_short_range(activity: str, expiration: str, mask: str, sr_presence: list, sr_activities: list):
+    if mask == 'No mask':
+        exposure_mask = models.Mask.types['No mask']
+    else:
+        exposure_mask = mask_distributions[mask]
+
+    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=500/100,
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                virus=virus_distributions['SARS_CoV_2_OMICRON'],
+                presence=mc.SpecificInterval(((8.5, 12.5),(13.5, 17.5),)),
+                mask=exposure_mask,
+                activity=activity_distributions[activity],
+                expiration=build_expiration(expiration),
+                host_immunity=0.,
+            ),
+        ),
+        short_range=mc.ShortRangeModel(
+            presence=[models.SpecificInterval(interval) for interval in sr_presence],
+            expirations=[short_range_expiration_distributions[activity] for activity in sr_activities],
+            dilutions=dilution_factor(activities=sr_activities,
+                        distance=np.random.uniform(0.5, 1.5, 250000)),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((8.5, 12.5),(13.5, 17.5),)),
+            activity=activity_distributions[activity],
+            mask=exposure_mask,
+            host_immunity=0.,
+        ),
+    )
+    return exposure_mc
+
+def exposure_module_with_short_range_outdoors(activity: str, expiration: str, mask: str, sr_presence: list, sr_activities: list):
+    if mask == 'No mask':
+        exposure_mask = models.Mask.types['No mask']
+    else:
+        exposure_mask = mask_distributions[mask]
+
+    exposure_mc = mc.ExposureModel(
+        concentration_model=mc.ConcentrationModel(
+            room=models.Room(volume=100000, humidity=0.5),
+            ventilation=models.AirChange(
+                active=models.SpecificInterval(((0, 24),)),
+                air_exch=1000000,
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                virus=virus_distributions['SARS_CoV_2_OMICRON'],
+                presence=mc.SpecificInterval(((8.5, 12.5),)),
+                mask=exposure_mask,
+                activity=activity_distributions[activity],
+                expiration=build_expiration(expiration),
+                host_immunity=0.,
+            ),
+        ),
+        short_range=mc.ShortRangeModel(
+            presence=[models.SpecificInterval(interval) for interval in sr_presence],
+            expirations=[short_range_expiration_distributions[activity] for activity in sr_activities],
+            dilutions=dilution_factor(activities=sr_activities,
+                        distance=np.random.uniform(0.5, 1.5, 250000)),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((8.5, 12.5),)),
+            activity=activity_distributions[activity],
+            mask=exposure_mask,
+            host_immunity=0.,
+        ),
+    )
+    return exposure_mc
+
+
+def baseline_model(activity: str, expiration: str, mask: str, sr_presence: list, sr_activities: list):
+    if mask == 'No mask':
+        exposure_mask = models.Mask.types['No mask']
+    else:
+        exposure_mask = mask_distributions[mask]
+
+    exposure_mc = mc.ExposureModel(
+        concentration_model=mc.ConcentrationModel(
+            room=models.Room(volume=100, humidity=0.5),
+            ventilation=models.AirChange(
+                active=models.PeriodicInterval(period=120, duration=120),
+                air_exch=0.25,
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                virus=virus_distributions['SARS_CoV_2_OMICRON'],
+                presence=models.SpecificInterval(((8.5, 9.5),)),
+                mask=exposure_mask,
+                activity=activity_distributions[activity],
+                expiration=build_expiration(expiration),
+                host_immunity=0.,
+            ),
+        ),
+        short_range=mc.ShortRangeModel(
+            presence=[models.SpecificInterval(interval) for interval in sr_presence],
+            expirations=[short_range_expiration_distributions[activity] for activity in sr_activities],
+            dilutions=dilution_factor(activities=sr_activities,
+                        distance=np.random.uniform(0.5, 1.5, 250000)),
+        ),
+        exposed=mc.Population(
+            number=3,
+            presence=models.SpecificInterval(((8.5, 9.5),)),
+            activity=activity_distributions[activity],
+            mask=exposure_mask,
+            host_immunity=0.,
+        ),
+    )
+    return exposure_mc

cara/short_range_plots/results.py

@@ -0,0 +1,85 @@
+""" Title: CARA - COVID Airborne Risk Assessment
+Author: A. Henriques et al
+Date: 18/02/2021
+Code version: 4.0.0
+"""
+
+from email.mime import base
+from cara.models import ExposureModel, InfectedPopulation
+from model_scenarios import *
+from scripts import *
+from itertools import product
+from dataclasses import dataclass
+from cara.monte_carlo.data import symptomatic_vl_frequencies
+
+# print('\n<<<<<<<<<<< Peak viral concentration without short range for baseline scenarios >>>>>>>>>>>')
+# concentration_curve(models = [exposure_module_without_short_range(activity='Light activity', expiration={"Speaking": 1, "Breathing": 2}, mask='No mask')],
+#                             labels = ['Background (long-range) concentration'],
+#                             labelsDose = ['Dose (long-range)'],
+#                             colors = ['royalblue'],
+#                             linestyles = ['--'],
+#                             thickness = [2])
+
+# print('\n<<<<<<<<<<< Peak viral concentration with short range interactions for baseline scenarios >>>>>>>>>>>')
+# concentration_curve(models=[exposure_module_with_short_range(
+#                                     activity='Light activity',
+#                                     expiration={"Speaking": 1, "Breathing": 2},
+#                                     mask='No mask',
+#                                     sr_presence=[(10.5, 11.0)],
+#                                     sr_activities=['Breathing']),
+#                             exposure_module_without_short_range(
+#                                 activity='Light activity',
+#                                 expiration={"Speaking": 1, "Breathing": 2},
+#                                 mask='No mask',)
+#                             ],
+#                             labels = ['Concentration with short range interactions', 'Background (long-range) concentration'],
+#                             labelsDose = ['Dose (full)', 'Dose (long-range)'],
+#                             colors = ['salmon', 'royalblue'],
+#                             linestyles = ['-', '--'],
+#                             thickness = [2, 2])
+
+# print('\n<<<<<<<<<<< Peak viral concentration with short range interactions for baseline scenarios >>>>>>>>>>>')
+# concentration_curve(models=[exposure_module_with_short_range(
+#                                     activity='Light activity',
+#                                     expiration={"Speaking": 1, "Breathing": 2},
+#                                     mask='No mask',
+#                                     sr_presence=[(10.5, 11.0)],
+#                                     sr_activities=['Speaking']),
+#                             exposure_module_without_short_range(
+#                                 activity='Light activity',
+#                                 expiration={"Speaking": 1, "Breathing": 2},
+#                                 mask='No mask',)
+#                             ],
+#                             labels = ['Concentration with short range interactions', 'Background (long-range) concentration'],
+#                             labelsDose = ['Dose (full)', 'Dose (long-range)'],
+#                             colors = ['salmon', 'royalblue'],
+#                             linestyles = ['-', '--'],
+#                             thickness = [2, 2])
+
+
+# print('\n<<<<<<<<<<< Dose vs SR exposure time >>>>>>>>>>>')
+#Always assume 1h for the short range interactions.
+#Always assume that in each model there is only ONE short range interaction.
+# plot_vD_vs_exposure_time(exp_models = [
+#                             baseline_model(
+#                                 activity='Light activity',
+#                                 expiration={"Speaking": 2, "Breathing": 1},
+#                                 mask='No mask',
+#                                 sr_presence=[(8.5, 9.5)],
+#                                 sr_activities=['Breathing']),
+#                             baseline_model(
+#                                 activity='Light activity',
+#                                 expiration={"Speaking": 2, "Breathing": 1},
+#                                 mask='No mask',
+#                                 sr_presence=[(8.5, 9.5)],
+#                                 sr_activities=['Speaking'])],
+#                          labels = ['Breathing', 'Speaking'],
+#                          colors=['royalblue', 'darkviolet'],
+#                          linestyles=['solid', 'solid'],
+#                          points=20,
+#                          time_in_minutes=True,
+#                          normalize_y_axis=True)
+
+print('\n<<<<<<<<<<< BLO curve >>>>>>>>>>>')
+generate_BLO_curve(activity='Speaking')
+    

cara/short_range_plots/scripts.py

@@ -0,0 +1,243 @@
+""" Title: CARA - COVID Airborne Risk Assessment
+Author: A. Henriques et al
+Date: 18/02/2021
+Code version: 4.0.0
+"""
+
+from cara.monte_carlo.data import short_range_expiration_distributions, expiration_distributions, expiration_BLO_factors
+from tqdm import tqdm
+from matplotlib.patches import Rectangle, Patch
+from scipy.spatial import ConvexHull
+from model_scenarios import *
+from cara.models import *
+import cara.monte_carlo as mc
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.lines as mlines
+import matplotlib.patches as patches
+import matplotlib as mpl
+from scipy.interpolate import make_interp_spline
+from mpl_toolkits.axes_grid1.inset_locator import mark_inset
+
+
+######### Plot material #########
+np.random.seed(2000)
+SAMPLE_SIZE = 250000
+TIMESTEP = 0.01
+#viral_loads = np.linspace(2, 12, 600)
+_VectorisedFloat = typing.Union[float, np.ndarray]
+
+
+def previous_deposited_exposure_between_bounds(model: ExposureModel, time1: float, time2: float) -> _VectorisedFloat:
+        """
+        The number of virus per m^3 deposited on the respiratory tract
+        between any two times.
+        """
+        emission_rate_per_aerosol = model.concentration_model.infected.emission_rate_per_aerosol_when_present()
+        aerosols = model.concentration_model.infected.aerosols()
+        fdep = model.long_range_fraction_deposited()
+        f_inf = model.concentration_model.infected.fraction_of_infectious_virus()
+
+        diameter = model.concentration_model.infected.particle.diameter
+
+        if not np.isscalar(diameter) and diameter is not None:
+            # we compute first the mean of all diameter-dependent quantities
+            # to perform properly the Monte-Carlo integration over
+            # particle diameters (doing things in another order would
+            # lead to wrong results).
+            dep_exposure_integrated = np.array(model._long_range_normed_exposure_between_bounds(time1, time2) *
+                                               aerosols *
+                                               fdep).mean()
+        else:
+            # in the case of a single diameter or no diameter defined,
+            # one should not take any mean at this stage.
+            dep_exposure_integrated = model._long_range_normed_exposure_between_bounds(time1, time2)*aerosols*fdep
+
+        # then we multiply by the diameter-independent quantity emission_rate_per_aerosol,
+        # and parameters of the vD equation (i.e. f_inf, BR_k and n_in).
+        return (dep_exposure_integrated * emission_rate_per_aerosol * 
+                f_inf * model.exposed.activity.inhalation_rate * 
+                (1 - model.exposed.mask.inhale_efficiency()))
+
+
+def concentration_curve(models, labels, labelsDose, colors, linestyles, thickness):
+
+    exp_models = [model.build_model(size=SAMPLE_SIZE) for model in models]
+
+    start = min(min(model.concentration_model.infected.presence.transition_times())
+                for model in exp_models)
+    stop = max(max(model.concentration_model.infected.presence.transition_times())
+               for model in exp_models)
+
+    times = np.arange(start, stop, TIMESTEP)
+
+    concentrations = [[np.mean(model.concentration(
+        t)) for t in times] for model in tqdm(exp_models)]
+    
+    fig, (ax, ax2) = plt.subplots(2, 1, sharex=True, figsize=(8,6))
+    
+    for c, color, linestyle, width in zip(concentrations, colors, linestyles, thickness):
+        ax.plot(times, c, color=color, ls=linestyle, lw=width)
+        ax2.plot(times, c, color=color, ls=linestyle, lw=width)
+
+    # zoom-in / limit the view to different portions of the data
+    ax.set_ylim(ax.get_ylim()[1]*0.8, ax.get_ylim()[1])  # outliers only
+    ax2.set_ylim(0, max(concentrations[1]))  # most of the data
+
+    ax.spines['bottom'].set_visible(False)
+    ax2.spines['top'].set_visible(False)
+    ax.xaxis.tick_top()
+    ax.tick_params(labeltop=False)  # don't put tick labels at the top
+    ax2.xaxis.tick_bottom()
+    #ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2)
+    # ax.set_ylim(ax.get_ylim()[0], 90)
+    ax.spines["right"].set_visible(False)
+
+    cumulative_doses = [np.cumsum([
+        np.array(model.deposited_exposure_between_bounds(
+            float(time1), float(time2))).mean()
+        for time1, time2 in tqdm(zip(times[:-1], times[1:]))
+    ]) for model in exp_models]
+
+    # quantile_05 = [np.cumsum([
+    #     np.quantile(np.array(model.deposited_exposure_between_bounds(
+    #         float(time1), float(time2))), 0.05)
+    #     for time1, time2 in tqdm(zip(times[:-1], times[1:]))
+    # ]) for model in exp_models]
+
+    # quantile_95 = [np.cumsum([
+    #     np.quantile(np.array(model.deposited_exposure_between_bounds(
+    #         float(time1), float(time2))), 0.95)
+    #     for time1, time2 in tqdm(zip(times[:-1], times[1:]))
+    # ]) for model in exp_models]
+
+    plt.xlabel("Time of day", fontsize=14)
+    plt.ylabel("Mean concentration\n(virions m$^{-3}$)", fontsize=14)
+
+    ax1 = ax.twinx()
+    ax11 = ax2.twinx()
+    for vd, color, width in tqdm(zip(cumulative_doses, colors, thickness)):
+        ax1.plot(times[:-1], vd,
+                 color=color, linestyle='dotted', lw=1)
+        ax11.plot(times[:-1], vd,
+                 color=color, linestyle='dotted', lw=1)
+        ax1.scatter([times[-1]], [vd[-1]], marker='.', color=color)
+        ax11.scatter([times[-1]], [vd[-1]], marker='.', color=color)
+
+    # zoom-in / limit the view to different portions of the data
+    ax1.set_ylim(ax1.get_ylim()[1]*0.8, ax1.get_ylim()[1])  # outliers only
+    ax11.set_ylim(0, max(cumulative_doses[1]))  # most of the data
+    ax11.spines["bottom"].set_visible(False)
+
+    # # Plot presence of exposed person
+    # for i, model in enumerate(models):
+    #     for i, (presence_start, presence_finish) in enumerate(model.exposed.presence.boundaries()):
+    #         plt.fill_between(
+    #             times, quantile_95[i], 0,
+    #             where=(np.array(times) > presence_start) & (np.array(times) < presence_finish),
+    #             color=color[i], alpha=0.1,
+    #         )
+
+    # # Plot short range interaction area
+    # for i, model in enumerate(models):
+    #     for presence in model.short_range.presence:
+    #         (presence_start, presence_finish) = presence.boundaries()
+    #         plt.fill_between(
+    #             times, quantile_95[i],
+    #             where=(np.array(times) > presence_start) & (np.array(times) < presence_finish),
+    #             color='#1f77b4', alpha=0.1,
+    #         )
+
+    
+    # ax1.spines["right"].set_linestyle((0, (1, 5)))
+    # ax1.set_ylabel('Mean cumulative dose\n(infectious virus)', fontsize=14)
+    #ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.3)
+    # ax1.set_ylim(ax1.get_ylim()[0], 40)
+    
+    complete_labels = labels + [label for label in labelsDose]
+    complete_colors = colors + [color for color in colors]
+    complete_linestyles = linestyles + ['dotted' for linestyle in linestyles]
+
+    labels_legend = [mlines.Line2D([], [], color=color, label=label, linestyle=linestyle)
+              for color, label, linestyle in zip(complete_colors, complete_labels, complete_linestyles)]
+
+    # for i in range(len(models)):
+    #     print('Scenario: ', labels[i])
+    #     print(
+    #         f"MEAN vD = {cumulative_doses[i][-1]}\n"
+    #         f"5th per = {quantile_05[i][-1]}\n"
+    #         f"95th per = {quantile_95[i][-1]}\n")
+    
+    ax.legend(handles=labels_legend, loc='upper right')
+
+    d = .015  # how big to make the diagonal lines in axes coordinates
+    # arguments to pass to plot, just so we don't keep repeating them
+    kwargs = dict(transform=ax.transAxes, color='k', clip_on=False)
+    ax.plot((-d, +d), (-d, +d), **kwargs)        # top-left diagonal
+    ax.plot((1 - d, 1 + d), (-d, +d), **kwargs)  # top-right diagonal
+
+    kwargs.update(transform=ax2.transAxes)  # switch to the bottom axes
+    ax2.plot((-d, +d), (1 - d, 1 + d), **kwargs)  # bottom-left diagonal
+    ax2.plot((1 - d, 1 + d), (1 - d, 1 + d), **kwargs)  # bottom-right diagonal
+
+    plt.show()
+
+def plot_vD_vs_exposure_time(exp_models: typing.List[mc.ExposureModel], labels, colors, linestyles, points: int = 20, time_in_minutes: bool = False, normalize_y_axis: bool = False) -> None:
+    
+    TIMESTEP = 0.001
+
+    concentration_models = [model.concentration_model for model in exp_models]
+    exposed_models = [model.exposed for model in exp_models]
+    
+    vDs: typing.List[typing.List[float]] = [[] for _ in exp_models]
+    
+    presence_intervals = [model.short_range.presence[0].boundaries() for model in exp_models]
+    start, final = presence_intervals[0]
+    times = np.linspace(start, final, points)
+    for finish in tqdm(times):
+        current_models = [mc.ExposureModel(
+            concentration_model=cm,
+            short_range=mc.ShortRangeModel(
+                presence=[models.SpecificInterval((start, finish), ),],
+                expirations=em.short_range.expirations,
+                dilutions=em.short_range.dilutions,
+            ),
+            exposed=exposed,
+        ) for cm, exposed, em in zip(concentration_models, exposed_models, exp_models)]
+
+        for i, m in enumerate(current_models):
+            vDs[i].append(np.mean(m.build_model(SAMPLE_SIZE).deposited_exposure()))
+
+    times = np.linspace(0, 60, points)
+    for i, vD in enumerate(vDs):
+        plt.plot(times, vD, color=colors[i], label=labels[i])
+
+    # plt.xlim((0, 60))
+    #if normalize_y_axis:
+    #    plt.ylim((0, 1))
+        
+    for m in exp_models:
+        print(np.mean(m.build_model(SAMPLE_SIZE).deposited_exposure()))
+
+    plt.xlabel(f'Duration of close-proximity encounter (min)', fontsize=12)
+    plt.ylabel('Mean cumulative dose\n(infectious virus)', fontsize=12)
+    plt.legend()
+    plt.show()
+
+def generate_BLO_curve(activity):
+    diameters = short_range_expiration_distributions[activity].diameter
+    BLOcurve = BLOmodel(expiration_BLO_factors[activity]).distribution(diameters
+            )/BLOmodel(expiration_BLO_factors[activity]).integrate(0.1,100)
+    # you have to divide by the integral to get the normalized distribution,
+    # so that you can compare with the normalized histogram
+
+    plt.xscale('log')
+    plt.yscale('log')
+
+    plt.plot(diameters,BLOcurve,'.',ms=5)
+    plt.hist(diameters,bins=100,range=(0,100),density=True)
+    plt.title(activity)
+
+    plt.xlabel("diameters (microns)")
+    plt.ylabel("Normalized distribution and histogram")
+    plt.show()