boxplot format

# Conflicts:
#	cara/model_scenarios.py

.gitignore

@@ -4,6 +4,8 @@ __pycache__
 *.DS_Store
 *.pyc
 
+data.csv
+
 # Editor stuff
 *.swp
 .idea

cara/apps/calculator/model_generator.py

@@ -246,7 +246,7 @@ class FormData:
                 ventilation=self.ventilation(),
                 infected=self.infected_population(),
             ),
-            exposed=self.exposed_population()
+            exposed=self.exposed_population(),
         )
 
     def build_model(self, sample_size=_DEFAULT_MC_SAMPLE_SIZE) -> models.ExposureModel:

cara/apps/calculator/report_generator.py

@@ -339,3 +339,4 @@ class ReportGenerator:
     def render(self, context: dict) -> str:
         template = self._template_environment().get_template("calculator.report.html.j2")
         return template.render(**context)
+

cara/model_scenarios_paper.py

@@ -0,0 +1,687 @@
+from cara import models
+from cara.monte_carlo.data import activity_distributions, symptomatic_vl_frequencies, infectious_virus_distribution, infectious_dose_distribution
+import cara.monte_carlo as mc
+import numpy as np
+
+######### Scatter points (data taken: copies per hour) #########
+
+############# 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 = [1668, 938, 319.6, 3632.8, 1243.6,
+                           17344, 2932, 5426, 5493.2, 1911.6, 9714.8]
+
+############# 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]
+############# 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 ###########
+
+######### Breathing model ###########
+def breathing_exposure(activity: str, mask: str):
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types[mask],
+                activity=activity_distributions[activity],
+                expiration=models.Expiration.types['Breathing'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types[activity],
+            mask=models.Mask.types[mask],
+        ),
+    )
+    return exposure_mc
+
+######### Speaking model ###########
+def speaking_exposure(activity: str, mask: str):
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types[mask],
+                activity=activity_distributions[activity],
+                expiration=models.Expiration.types['Talking'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types[activity],
+            mask=models.Mask.types[mask],
+        ),
+    )
+    return exposure_mc
+
+######### Shouting model ###########
+def shouting_exposure(activity: str, mask: str):
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types[mask],
+                activity=activity_distributions[activity],
+                expiration=models.Expiration.types['Shouting'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types[activity],
+            mask=models.Mask.types[mask],
+        ),
+    )
+    return exposure_mc
+
+######### Breathing model for specific viral load ###########
+def breathing_exposure_vl(vl):
+    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=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                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"],
+        ),
+    )
+    return exposure_mc
+
+######### Talking model for specific viral load ###########
+def talking_exposure_vl(vl):
+    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=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                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"],
+        ),
+    )
+    return exposure_mc
+
+######### Shouting model for specific viral load ###########
+def shouting_exposure_vl(vl):
+    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=mc.Virus(
+                    viral_load_in_sputum=10**vl,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types['No mask'],
+                activity=activity_distributions['Light activity'],
+                expiration=models.Expiration.types['Shouting'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Light activity'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+######### Used for CDF Models ###########
+######### Breathing Models #########
+def breathing_seated_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                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"],
+        ),
+    )
+    return exposure_mc
+
+def breathing_light_activity_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Light activity'],
+                expiration=models.Expiration.types['Breathing'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Light activity'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+def breathing_heavy_exercise_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Heavy exercise'],
+                expiration=models.Expiration.types['Breathing'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Heavy exercise'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+######### Speaking Models #########
+def speaking_seated_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                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"],
+        ),
+    )
+    return exposure_mc
+
+def speaking_light_activity_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Light activity'],
+                expiration=models.Expiration.types['Talking'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Light activity'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+def speaking_heavy_exercise_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Heavy exercise'],
+                expiration=models.Expiration.types['Talking'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Heavy exercise'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+######### Shouting Models #########
+def shouting_seated_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Seated'],
+                expiration=models.Expiration.types['Shouting'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Seated'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+def shouting_light_activity_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Light activity'],
+                expiration=models.Expiration.types['Shouting'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Light activity'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+def shouting_heavy_exercise_exposure():
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Heavy exercise'],
+                expiration=models.Expiration.types['Shouting'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types['Heavy exercise'],
+            mask=models.Mask.types["No mask"],
+        ),
+    )
+    return exposure_mc
+
+########## Concentration curves ###########
+def office_model_no_mask_windows_closed():
+    office_model_no_vent = mc.ExposureModel(
+        concentration_model=mc.ConcentrationModel(
+            room=models.Room(volume=160, humidity=0.3),
+            ventilation=models.MultipleVentilation(
+                (models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.0), 
+                models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25))),
+            infected=mc.InfectedPopulation(
+                number=1,
+                presence=models.SpecificInterval(present_times = ((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
+                virus=mc.SARSCoV2(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Seated'],
+                expiration=models.MultipleExpiration(
+                    expirations = (models.Expiration.types['Talking'],
+                                models.Expiration.types['Breathing']),
+                    weights=(1, 2)
+                )
+            )
+        ),
+        exposed=models.Population(
+            number=18,
+            presence=models.SpecificInterval(present_times = ((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
+            activity=models.Activity.types['Seated'],
+            mask=models.Mask.types['No mask']
+        )
+    )
+    return office_model_no_vent
+
+def office_model_no_mask_windows_open_breaks():
+    office_model_no_vent = mc.ExposureModel(
+        concentration_model=mc.ConcentrationModel(
+            room=models.Room(volume=160, humidity=0.3),
+            ventilation = models.MultipleVentilation(
+                ventilations=(
+                    models.SlidingWindow(
+                        active=models.SpecificInterval(present_times=((1.5, 2), (3.5, 4.5), (6, 6.5))),
+                        inside_temp=models.PiecewiseConstant((0, 24), (295,)),
+                        outside_temp=models.PiecewiseConstant((0, 24), (291,)),
+                        window_height=1.6, 
+                        opening_length=0.6,
+                    ),
+                    models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25),
+                )  
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
+                virus=mc.SARSCoV2(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Seated'],
+                expiration=models.MultipleExpiration(
+                    expirations = (models.Expiration.types['Talking'],
+                                models.Expiration.types['Breathing']),
+                    weights=(1, 2)
+                )
+            )
+        ),
+        exposed=models.Population(
+            number=18,
+            presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
+            activity=models.Activity.types['Seated'],
+            mask=models.Mask.types['No mask']
+        )
+    )
+    return office_model_no_vent
+
+def office_model_no_mask_windows_open_alltimes():
+    office_model_no_vent = mc.ExposureModel(
+        concentration_model=mc.ConcentrationModel(
+            room=models.Room(volume=160, humidity=0.3),
+            ventilation=models.MultipleVentilation(
+                ventilations=(
+                    models.SlidingWindow(
+                        active=models.PeriodicInterval(period=120, duration=120),
+                        inside_temp=models.PiecewiseConstant((0, 24), (295,)),
+                        outside_temp=models.PiecewiseConstant((0, 24), (291,)),
+                        window_height=1.6, opening_length=0.6,
+                    ),
+                    models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25),
+                )
+            ),
+            infected=mc.InfectedPopulation(
+                number=1,
+                presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
+                virus=mc.SARSCoV2(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                mask=models.Mask.types["No mask"],
+                activity=activity_distributions['Seated'],
+                expiration=models.MultipleExpiration(
+                    expirations = (models.Expiration.types['Talking'],
+                                models.Expiration.types['Breathing']),
+                    weights=(1, 2)
+                )
+            )
+        ),
+        exposed=models.Population(
+            number=18,
+            presence=models.SpecificInterval(present_times=((0, 1.5), (2, 3.5), (4.5, 6), (6.5, 8))),
+            activity=models.Activity.types['Seated'],
+            mask=models.Mask.types['No mask']
+        )
+    )
+    return office_model_no_vent
+
+
+######### Standard exposure models ###########
+
+######### Breathing model ###########
+def breathing_exposure(activity: str, mask: str):
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types[mask],
+                activity=activity_distributions[activity],
+                expiration=models.Expiration.types['Breathing'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types[activity],
+            mask=models.Mask.types[mask],
+        ),
+    )
+    return exposure_mc
+
+######### Speaking model ###########
+def speaking_exposure(activity: str, mask: str):
+    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=mc.Virus(
+                    viral_load_in_sputum=symptomatic_vl_frequencies,
+                    infectious_dose=infectious_dose_distribution,
+                    viable_to_RNA=infectious_virus_distribution,
+                ),
+                presence=mc.SpecificInterval(((0, 2),)),
+                mask=models.Mask.types[mask],
+                activity=activity_distributions[activity],
+                expiration=models.Expiration.types['Talking'],
+            ),
+        ),
+        exposed=mc.Population(
+            number=14,
+            presence=mc.SpecificInterval(((0, 2),)),
+            activity=models.Activity.types[activity],
+            mask=models.Mask.types[mask],
+        ),
+    )
+    return exposure_mc
+
+######### Infected Population model ###########
+def infected_model(mask: str, activity: str, expiratory_activity: str):
+    infected=mc.InfectedPopulation(
+        number=1,
+        virus=mc.Virus(
+            viral_load_in_sputum=symptomatic_vl_frequencies,
+            infectious_dose=infectious_dose_distribution,
+            viable_to_RNA=infectious_virus_distribution,
+        ),
+        presence=mc.SpecificInterval(((0, 2),)),
+        mask=models.Mask.types[mask],
+        activity=activity_distributions[activity],
+        expiration=models.Expiration.types[expiratory_activity])
+
+    return infected

cara/models.py

@@ -50,7 +50,6 @@ from .utils import method_cache
 
 from .dataclass_utils import nested_replace
 
-
 # Define types for items supporting vectorisation. In the future this may be replaced
 # by ``np.ndarray[<type>]`` once/if that syntax is supported. Note that vectorization
 # implies 1d arrays: multi-dimensional arrays are not supported.
@@ -429,6 +428,9 @@ class Virus:
     #: Dose to initiate infection, in RNA copies
     infectious_dose: _VectorisedFloat
 
+    #: viable-to-RNA virus ratio as a function of the viral load
+    viable_to_RNA: _VectorisedFloat
+
     #: Pre-populated examples of Viruses.
     types: typing.ClassVar[typing.Dict[str, "Virus"]]
 
@@ -465,19 +467,23 @@ Virus.types = {
         # as per https://www.dhs.gov/publication/st-master-question-list-covid-19
         # 50 comes from Buonanno et al.
         infectious_dose=50.,
+        viable_to_RNA = 0.5,
     ),
     'SARS_CoV_2_B117': SARSCoV2(
         # also called VOC-202012/01
         viral_load_in_sputum=1e9,
         infectious_dose=30.,
+        viable_to_RNA = 0.5,
     ),
     'SARS_CoV_2_P1': SARSCoV2(
         viral_load_in_sputum=1e9,
         infectious_dose=1/0.045,
+        viable_to_RNA = 0.5,
     ),
     'SARS_CoV_2_B16172': SARSCoV2(
         viral_load_in_sputum=1e9,
         infectious_dose=30/1.6,
+        viable_to_RNA = 0.5,
     ),
 }
 
@@ -539,7 +545,7 @@ class _ExpirationBase:
     #: Pre-populated examples of Expirations.
     types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
 
-    def aerosols(self, mask: Mask):
+    def aerosols(self, mask: Mask, cn_B: float, cn_L: float):
         """
         total volume of aerosols expired per volume of air (mL/cm^3).
         """
@@ -561,19 +567,19 @@ class Expiration(_ExpirationBase):
     BLO_factors: typing.Tuple[float, float, float]
 
     @cached()
-    def aerosols(self, mask: Mask):
+    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) -> float:
+        def _Lmode(d: float, cn_L: float) -> float:
             # L-mode (see ref. above).
-            return ( (1 / d) * (1.0 / (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:
@@ -582,8 +588,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) +
+            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))
 
@@ -608,9 +614,9 @@ class MultipleExpiration(_ExpirationBase):
             raise ValueError("expirations and weigths should contain the"
                              "same number of elements")
 
-    def aerosols(self, mask: Mask):
+    def aerosols(self, mask: Mask, cn_B: float, cn_L: float):
         return np.array([
-            weight * expiration.aerosols(mask) / sum(self.weights)
+            weight * expiration.aerosols(mask, cn_B, cn_L) / sum(self.weights)
             for weight,expiration in zip(self.weights,self.expirations)
         ]).sum(axis=0)
 
@@ -676,8 +682,11 @@ class InfectedPopulation(Population):
     #: The type of expiration that is being emitted whilst doing the activity.
     expiration: _ExpirationBase
 
-    @method_cache
-    def emission_rate_when_present(self) -> _VectorisedFloat:
+    #: The percentage of host immunity
+    host_immunity: float = 0.
+
+   
+    def emission_rate_when_present(self, cn_B: float = 0.06, cn_L: float = 0.2) -> _VectorisedFloat:
         """
         The emission rate if the infected population is present.
 
@@ -687,7 +696,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)
+        aerosols = self.expiration.aerosols(self.mask, cn_B, cn_L)
 
         ER = (self.virus.viral_load_in_sputum *
               self.activity.exhalation_rate *
@@ -719,7 +728,7 @@ class InfectedPopulation(Population):
         # with a declaration of state change time, as is the case for things
         # like Ventilation.
 
-        return self.emission_rate_when_present()
+        return self.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
 
 
 @dataclass(frozen=True)
@@ -906,7 +915,33 @@ class ExposureModel:
     repeats: int = 1
 
     #: The fraction of viruses actually deposited in the respiratory tract
-    fraction_deposited: _VectorisedFloat = 0.6
+    """
+    To be updated in the future.
+    The diameter value of 1.37 is the correspondent to a fraction_deposited of 0.6
+    """
+    d = 1.37 
+    IF = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8)))))
+    DF = IF * (0.0587 + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d)))))
+    fraction_deposited: _VectorisedFloat = DF
+
+    def _normed_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
+        """The number of virions per meter^3 between any two times, normalized 
+        by the emission rate of the infected population"""
+        for start, stop in self.exposed.presence.boundaries():
+            if start > time2:
+                normed_exposure = 0.
+                break
+            elif time2 <= stop:
+                normed_exposure = self.concentration_model.normed_integrated_concentration(time1, time2)
+                break
+            else:
+                normed_exposure = self.concentration_model.normed_integrated_concentration(time1, time2)
+        return normed_exposure
+
+    def exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
+        """The number of virions per meter^3 between any two times."""
+        return (self._normed_exposure_between_bounds(time1, time2) * 
+                self.concentration_model.infected.emission_rate_when_present())
 
     def _normed_exposure(self) -> _VectorisedFloat:
         """
@@ -928,10 +963,11 @@ class ExposureModel:
     def infection_probability(self) -> _VectorisedFloat:
         exposure = self.exposure()
 
+        # Dose
         inf_aero = (
             self.exposed.activity.inhalation_rate *
             (1 - self.exposed.mask.inhale_efficiency()) *
-            exposure * self.fraction_deposited
+            exposure * self.fraction_deposited * (self.concentration_model.infected.virus.viable_to_RNA * (1 - self.concentration_model.infected.host_immunity))
         )
 
         # Probability of infection.

cara/monte_carlo/data.py

@@ -38,24 +38,33 @@ symptomatic_vl_frequencies = LogCustomKernel(
     kernel_bandwidth=0.1
 )
 
+# From https://doi.org/10.1093/cid/ciaa1579
+infectious_virus_distribution = Uniform(0.15, 0.45)
+
+# From discussion with virologists
+infectious_dose_distribution = Uniform(10., 100.)
 
 # From CERN-OPEN-2021-04 and refererences therein
 virus_distributions = {
     'SARS_CoV_2': mc.SARSCoV2(
                 viral_load_in_sputum=symptomatic_vl_frequencies,
-                infectious_dose=100,
+                infectious_dose=infectious_dose_distribution,
+                viable_to_RNA=infectious_virus_distribution,
                 ),
     'SARS_CoV_2_B117': mc.SARSCoV2(
                 viral_load_in_sputum=symptomatic_vl_frequencies,
-                infectious_dose=60,
+                infectious_dose=infectious_dose_distribution,
+                viable_to_RNA=infectious_virus_distribution,
                 ),
     'SARS_CoV_2_P1': mc.SARSCoV2(
                 viral_load_in_sputum=symptomatic_vl_frequencies,
-                infectious_dose=100/2.25,
+                infectious_dose=infectious_dose_distribution,
+                viable_to_RNA=infectious_virus_distribution,
                 ),
     'SARS_CoV_2_B16172': mc.SARSCoV2(
                 viral_load_in_sputum=symptomatic_vl_frequencies,
-                infectious_dose=60/1.6,
+                infectious_dose=infectious_dose_distribution,
+                viable_to_RNA=infectious_virus_distribution,
                 ),
 }
 
@@ -68,3 +77,5 @@ mask_distributions = {
     'Type I': mc.Mask(Uniform(0.25, 0.80)),
     'FFP2': mc.Mask(Uniform(0.83, 0.91)),
 }
+
+

cara/plot_output.py

@@ -0,0 +1,55 @@
+""" Title: COVID Airborne Risk Assessment
+Author: <author(s) names>
+Date: <date>
+Code version: <code version>
+Availability: <where it's located> """
+
+from cara.models import ExposureModel, InfectedPopulation
+from cara import model_scenarios_paper
+from cara.results_paper import *
+from cara.test_plots import *
+from cara.monte_carlo.data import symptomatic_vl_frequencies
+from itertools import product
+from dataclasses import dataclass
+
+# Exhaled virions while talking, seated #
+#print('\n<<<<<<<<<<< Vlout for Talking, seated >>>>>>>>>>>')
+#exposure_model_from_vl_talking()
+
+# Exhaled virions while breathing, seated #
+#print('\n<<<<<<<<<<< Vlout for Breathing, seated >>>>>>>>>>>')
+#exposure_model_from_vl_breathing()
+
+# Exhaled virions while breathing, light activity #
+#print('\n<<<<<<<<<<< Vlout for Shouting, light activity >>>>>>>>>>>')
+#exposure_model_from_vl_shouting()
+
+# Exhaled virions while talking according to BLO model, seated #
+#print('\n<<<<<<<<<<< Vlout for Talking, seated with chosen Cn,L >>>>>>>>>>>')
+#exposure_model_from_vl_talking_cn()
+
+# Exhaled virions while breathing according to BLO model, seated #
+#print('\n<<<<<<<<<<< Vlout for Breathing, seated with chosen Cn,B >>>>>>>>>>>')
+#exposure_model_from_vl_breathing_cn()
+#print('\n')
+
+############ Plots with viral loads and emission rates + statistical data ############
+#present_vl_er_histograms(activity='Seated', mask='No mask')
+#present_vl_er_histograms(activity='Light activity', mask='No mask')
+#present_vl_er_histograms(activity='Heavy exercise', mask='No mask')
+
+############ CDFs for comparing the QR-Values in different scenarios ############
+#generate_cdf_curves()
+
+############ Deposition Fraction Graph ############
+#print('\n<<<<<<<<<<< Deposition Fraction for Breathing, seated >>>>>>>>>>>') 
+#calculate_deposition_factor()
+
+############ Comparison between concentration curves ############ 
+#compare_concentration_curves()
+
+############ Emission Rate Violin plot ############ 
+compare_viruses_vr()
+
+############ Used for testing ############
+#exposure_model_from_vl_talking_new_points()

cara/results_paper.py

@@ -0,0 +1,822 @@
+from numpy.core.function_base import linspace
+from cara import models
+from cara.monte_carlo.data import activity_distributions
+from tqdm import tqdm
+from matplotlib.patches import Rectangle, Patch
+from scipy.spatial import ConvexHull
+from model_scenarios_paper import *
+import cara.monte_carlo as mc
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+import matplotlib.lines as mlines
+import matplotlib.patches as patches
+import matplotlib as mpl
+from random import randrange
+
+
+######### Plot material #########
+SAMPLE_SIZE = 250000
+viral_loads = np.linspace(2, 12, 600)
+
+############# Markers (for legend) #############
+markers = [5, 'd', 4]
+
+""" Exhaled virions from exposure models """
+
+######### Breathing #########
+
+
+def exposure_model_from_vl_breathing():
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+
+    er_means = []
+    er_means_1h = []
+    er_medians = []
+    lower_percentiles = []
+    upper_percentiles = []
+
+    for vl in tqdm(viral_loads):
+        exposure_mc = breathing_exposure_vl(vl)
+        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=0.06, cn_L=0.2) / 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))
+        emission_rate_1h = exposure_model.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
+        er_means_1h.append(np.mean(emission_rate_1h))
+
+    # 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)
+    ax.fill_between(viral_loads, lower_percentiles,
+                    upper_percentiles, alpha=0.2)
+    ax.set_yscale('log')
+
+    ratio = np.mean(10**viral_loads / er_means)
+    ratio_1h = np.mean(10**viral_loads / er_means_1h)
+    print('Mean swab-to-aersol vl ratio in 30min:')
+    print(format(ratio, "5.1e"))
+    print('Mean swab-to-aersol vl ratio emission rate per hour:')
+    print(format(ratio_1h, "5.1e"))
+
+    ############# Coleman #############
+    scatter_coleman_data(coleman_etal_vl_breathing,
+                         coleman_etal_er_breathing_2)
+
+    ############# Milton et al #############
+    scatter_milton_data(milton_vl, milton_er_2)
+
+    ############# Yan et al #############
+    scatter_yann_data(yann_vl, yann_er_2)
+
+    ############ Legend ############
+    build_breathing_legend(fig)
+
+    ############ Plot ############
+    plt.title('',
+              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()
+
+
+""" Variation according to the BLO model """
+############ Breathing ############
+
+
+def exposure_model_from_vl_breathing_cn():
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+
+    n_lines = 30
+    cns = np.linspace(0.01, 0.5, n_lines)
+
+    cmap = define_colormap(cns)
+
+    for cn in tqdm(cns):
+        er_means = np.array([])
+        for vl in viral_loads:
+            exposure_mc = breathing_exposure_vl(vl)
+            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=0.2) / 2
+            er_means = np.append(er_means, np.mean(emission_rate))
+
+        # 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, alpha=0.75), linewidth=0.5)
+
+    # The dashed line for the chosen Cn,B
+    er_means = np.array([])
+    for vl in viral_loads:
+        exposure_mc = breathing_exposure_vl(vl)
+        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=0.06, cn_L=0.2) / 2
+        er_means = np.append(er_means, np.mean(emission_rate))
+    ax.plot(viral_loads, er_means, color=cmap.to_rgba(
+        cn, alpha=0.75), linewidth=1, ls='--')
+    plt.text(viral_loads[int(len(viral_loads)*0.9)], 10**4.2,
+             r"$\mathbf{c_{n,B}=0.06}$", color=cmap.to_rgba(cn), fontsize=12)
+
+    cmap = fig.colorbar(cmap, ticks=[0.01, 0.1, 0.25, 0.5])
+    cmap.set_label(label='Particle emission concentration, ${c_{n,B}}$', fontsize=12)
+    ax.set_yscale('log')
+
+    ############# Coleman #############
+    scatter_coleman_data(coleman_etal_vl_breathing,
+                         coleman_etal_er_breathing_2)
+
+    ############# Milton et al #############
+    scatter_milton_data(milton_vl, milton_er_2)
+
+    ############# Yan et al #############
+    scatter_yann_data(yann_vl, yann_er_2)
+
+    ############ Legend ############
+    build_breathing_legend(fig)
+
+    ############ Plot ############
+    plt.title('',
+              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()
+
+############ Talking ############
+
+
+def exposure_model_from_vl_talking():
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+
+    er_means = []
+    er_means_1h = []
+    er_medians = []
+    lower_percentiles = []
+    upper_percentiles = []
+
+    for vl in tqdm(viral_loads):
+        exposure_mc = talking_exposure_vl(vl)
+        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(cn_B=0.06, cn_L=0.2)/4
+        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))
+        emission_rate_1h = exposure_model.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
+        er_means_1h.append(np.mean(emission_rate_1h))
+
+    # divide by 4 to have in 15min (quarter of an hour)
+    coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
+
+    ax.plot(viral_loads, er_means)
+    ax.fill_between(viral_loads, lower_percentiles,
+                    upper_percentiles, alpha=0.2)
+    ax.set_yscale('log')
+
+    ratio = np.mean(10**viral_loads / er_means)
+    ratio_1h = np.mean(10**viral_loads / er_means_1h)
+    print('Mean swab-to-aersol vl ratio in 30min:')
+    print(format(ratio, "5.1e"))
+    print('Mean swab-to-aersol vl ratio emission rate per hour:')
+    print(format(ratio_1h, "5.1e"))
+
+    ############# Coleman #############
+    scatter_coleman_data(coleman_etal_vl_talking, coleman_etal_er_talking_2)
+
+    ############ Legend ############
+    build_talking_legend(fig)
+
+    ############ Plot ############
+    plt.title('',
+              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()
+
+def exposure_model_from_vl_talking_cn():
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+
+    n_lines = 30
+    cns = np.linspace(0.01, 2, n_lines)
+    cmap = define_colormap(cns)
+
+    for cn in tqdm(cns):
+        er_means = np.array([])
+        for vl in viral_loads:
+            exposure_mc = talking_exposure_vl(vl)
+            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(
+                cn_B=0.1, cn_L=cn) / 4
+            er_means = np.append(er_means, np.mean(emission_rate))
+
+        # divide by 4 to have in 15min (quarter of an hour)
+        coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
+        ax.plot(viral_loads, er_means, color=cmap.to_rgba(
+            cn, alpha=0.75), linewidth=0.5)
+
+    # The dashed line for the chosen Cn,L
+    er_means = np.array([])
+    for vl in viral_loads:
+        exposure_mc = talking_exposure_vl(vl)
+        exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+        # divide by 4 to have in 15min
+        emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
+            cn_B=0.06, cn_L=0.2) / 4
+        er_means = np.append(er_means, np.mean(emission_rate))
+    ax.plot(viral_loads, er_means, color=cmap.to_rgba(
+        cn, alpha=0.75), linewidth=1, ls='--')
+    plt.text(viral_loads[int(len(viral_loads)*0.93)], 10**5.5,
+             r"$\mathbf{c_{n,L}=0.2}$", color=cmap.to_rgba(cn), fontsize=12)
+
+    cmap = fig.colorbar(cmap, ticks=[0.01, 0.5, 1.0, 2.0])
+    cmap.set_label(label='Particle emission concentration, ${c_{n,L}}$', fontsize=12)
+    ax.set_yscale('log')
+
+    ############# Coleman #############
+    scatter_coleman_data(coleman_etal_vl_talking, coleman_etal_er_talking_2)
+
+    ############ Legend ############
+    build_talking_legend(fig)
+
+    ############ Plot ############
+    plt.title('',
+              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()
+
+####### Shouting ########
+
+
+def exposure_model_from_vl_shouting():
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+
+    er_means_1h = []
+
+    for vl in tqdm(viral_loads):
+        exposure_mc = shouting_exposure_vl(vl)
+        exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+        emission_rate_1h = exposure_model.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)
+        er_means_1h.append(np.mean(emission_rate_1h))
+
+    ax.plot(viral_loads, er_means_1h)
+    ax.set_yscale('log')
+
+    ratio_1h = np.mean(10**viral_loads / er_means_1h)
+    print('Mean swab-to-aersol vl ratio emission rate per hour:')
+    print(format(ratio_1h, "5.1e"))
+
+    ############ Plot ############
+    plt.title('',
+              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()
+
+
+############ Plots with viral loads and emission rates ############
+############ Statistical Data ############
+
+
+def get_statistical_data(activity: str, mask: str):
+    ############ Breathing model ############
+    exposure_mc = breathing_exposure(activity, mask)
+    exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+    emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
+        cn_B=0.06, cn_L=0.2)
+    breathing_er = [np.log10(er) for er in emission_rate]
+    print('\n<<<<<<<<<<< Breathing model statistics >>>>>>>>>>>')
+    print_er_info(emission_rate, breathing_er)
+
+    ############ Speaking model ############
+    exposure_mc = speaking_exposure(activity, mask)
+    exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+    emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
+        cn_B=0.06, cn_L=0.2)
+    speaking_er = [np.log10(er) for er in emission_rate]
+    print('\n<<<<<<<<<<< Speaking model statistics >>>>>>>>>>>')
+    print_er_info(emission_rate, speaking_er)
+
+    ############ Shouting model ############
+    exposure_mc = shouting_exposure(activity, mask)
+    exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)
+    emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(
+        cn_B=0.06, cn_L=0.2)
+    shouting_er = [np.log10(er) for er in emission_rate]
+    print('\n<<<<<<<<<<< Shouting model statistics >>>>>>>>>>>')
+    print_er_info(emission_rate, shouting_er)
+
+    viral_load_in_sputum = exposure_model.concentration_model.infected.virus.viral_load_in_sputum
+
+    return viral_load_in_sputum, breathing_er, speaking_er, shouting_er
+
+
+def present_vl_er_histograms(activity: str, mask: str):
+    viral_load_in_sputum, breathing_er, speaking_er, shouting_er = get_statistical_data(activity, mask)
+
+    fig, axs = plt.subplots(1, 2, sharex=False, sharey=False)
+    fig.set_figheight(5)
+    fig.set_figwidth(10)
+    plt.tight_layout()
+    plt.subplots_adjust(wspace=0.2)
+    plt.subplots_adjust(top=0.9)
+    plt.subplots_adjust(bottom=0.15)
+
+    viral_loads = [np.log10(vl) for vl in viral_load_in_sputum]
+
+    axs[0].hist(viral_loads, bins = 300, color='lightgrey')
+    axs[0].set_xlabel('vl$_{\mathrm{in}}$ (log$_{10}$ RNA copies mL$^{-1}$)')
+
+    mean = np.mean(viral_loads)
+    axs[0].vlines(x=(mean), ymin=0, ymax=axs[0].get_ylim()[1], colors=('grey'), linestyles=('dashed'))
+
+    breathing_mean_er = np.mean(breathing_er)
+    speaking_mean_er = np.mean(speaking_er)
+    shouting_mean_er = np.mean(shouting_er)
+
+    axs[1].hist(breathing_er, bins = 300, color='lightsteelblue')
+    axs[1].hist(speaking_er, bins = 300, color='wheat')
+    axs[1].hist(shouting_er, bins = 300, color='darkseagreen')
+    axs[1].set_xlabel('vR (log$_{10}$ virions h$^{-1}$)')
+
+    axs[1].vlines(x=(breathing_mean_er, speaking_mean_er, shouting_mean_er), ymin=0, ymax=axs[1].get_ylim()[1], colors=('cornflowerblue', 'goldenrod', 'olivedrab'), alpha=(0.75, 0.75, 0.75), linestyles=('dashed', 'dashed', 'dashed'))
+
+    labels = [Patch([], [], color=color, label=label)
+             for color, label in zip(['lightgrey', 'lightsteelblue', 'wheat', 'darkseagreen'],
+                                            ['Viral Load', 'Breathing', 'Speaking', 'Shouting'])]
+    labels.append(mlines.Line2D([], [], color='black',
+                            marker='', linestyle='dashed', label='Mean'))
+
+
+    for x in (0, 1):
+        axs[x].set_yticklabels([])
+        axs[x].set_yticks([])
+
+    plt.legend(handles=labels, loc='upper left', bbox_to_anchor=(1, 1))
+    plt.tight_layout()
+    plt.show()
+
+############ CDFs for comparing the QR-Values in different scenarios ############
+
+
+def generate_cdf_curves():
+    fig, axs = plt.subplots(3, 1, sharex='all')
+
+    ############ Breathing models ############
+    br_seated = breathing_seated_exposure()
+    br_seated_model = br_seated.build_model(size=SAMPLE_SIZE)
+
+    br_light_activity = breathing_light_activity_exposure()
+    br_light_activity_model = br_light_activity.build_model(size=SAMPLE_SIZE)
+
+    br_heavy_exercise = breathing_heavy_exercise_exposure()
+    br_heavy_exercise_model = br_heavy_exercise.build_model(size=SAMPLE_SIZE)
+
+    ############ Speaking models ############
+    sp_seated = speaking_seated_exposure()
+    sp_seated_model = sp_seated.build_model(size=SAMPLE_SIZE)
+
+    sp_light_activity = speaking_light_activity_exposure()
+    sp_light_activity_model = sp_light_activity.build_model(size=SAMPLE_SIZE)
+
+    sp_heavy_exercise = speaking_heavy_exercise_exposure()
+    sp_heavy_exercise_model = sp_heavy_exercise.build_model(size=SAMPLE_SIZE)
+
+    ############ Shouting models ############
+    sh_seated = shouting_seated_exposure()
+    sh_seated_model = sh_seated.build_model(size=SAMPLE_SIZE)
+
+    sh_light_activity = shouting_light_activity_exposure()
+    sh_light_activity_model = sh_light_activity.build_model(size=SAMPLE_SIZE)
+
+    sh_heavy_exercise = shouting_heavy_exercise_exposure()
+    sh_heavy_exercise_model = sh_heavy_exercise.build_model(size=SAMPLE_SIZE)
+
+    er_values = [np.log10(scenario.concentration_model.infected.emission_rate_when_present(cn_B=0.06, cn_L=0.2)) for scenario in (br_seated_model, br_light_activity_model,
+                                                                                                                                  br_heavy_exercise_model, sp_seated_model, sp_light_activity_model, sp_heavy_exercise_model, sh_seated_model, sh_light_activity_model, sh_heavy_exercise_model)]
+    left = min(np.min(ers) for ers in er_values)
+    right = max(np.max(ers) for ers in er_values)
+
+    # Colors can be changed here
+    colors_breathing = ['lightsteelblue', 'cornflowerblue', 'royalblue']
+    colors_speaking = ['wheat', 'tan', 'orange']
+    colors_shouting = ['palegreen', 'darkseagreen', 'forestgreen']
+    colors = [colors_breathing, colors_speaking, colors_shouting]
+
+    breathing_rates = ['Seated', 'Light activity', 'Heavy activity']
+    activities = ['Breathing', 'Speaking', 'Shouting']
+    lines_breathing = [mlines.Line2D([], [], color=color, markersize=15, label=label)
+                       for color, label in zip(colors_breathing, breathing_rates)]
+    lines_speaking = [mlines.Line2D([], [], color=color, markersize=15, label=label)
+                      for color, label in zip(colors_speaking, breathing_rates)]
+    lines_shouting = [mlines.Line2D([], [], color=color, markersize=15, label=label)
+                      for color, label in zip(colors_shouting, breathing_rates)]
+    lines = [lines_breathing, lines_speaking, lines_shouting]
+
+    samples: int = 50000
+    for i in range(3):
+        axs[i].hist(er_values[3 * i:3 * (i + 1)], bins=2000,
+                    histtype='step', cumulative=True, range=(-7, 6), color=colors[i])
+        axs[i].set_xlim(-6, 6)
+        axs[i].set_yticks([0, samples / 2, samples])
+        axs[i].set_yticklabels(['0.0', '0.5', '1.0'])
+        axs[i].yaxis.set_label_position("right")
+        axs[i].set_ylabel(activities[i], fontsize=14)
+        axs[i].grid(linestyle='--')
+        axs[i].legend(handles=lines[i], loc='upper left')
+
+    plt.xlabel('$\mathrm{vR}$', fontsize=16)
+    tick_positions = np.arange(-6, 6, 2)
+    plt.xticks(ticks=tick_positions, labels=[
+               '$\;10^{' + str(i) + '}$' for i in tick_positions])
+
+    fig.text(0.02, 0.5, 'Cumulative Distribution Function',
+             va='center', rotation='vertical', fontsize=14)
+    fig.set_figheight(8)
+    fig.set_figwidth(5)
+    plt.show()
+
+############ Deposition Fraction Graph ############# 
+def calculate_deposition_factor():
+
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+    
+    diameters = np.linspace(0.001, 0.01, 1000) #particle diameter (μm)
+    diameters = np.append(diameters, np.linspace(0.01, 0.1, 1000))
+    diameters = np.append(diameters, np.linspace(0.1, 1., 1000))
+    diameters = np.append(diameters, np.linspace(1., 10., 1000))
+    diameters = np.append(diameters, np.linspace(10, 100, 1000))
+
+    fractions_et = []
+    fractions_tb = []
+    fractions_al = []
+    fractions_df = []
+    for d in diameters:
+        
+        IF = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8)))))
+        DF_et = IF * (
+            (1 / ( 1 + np.exp(6.84 + 1.183 * np.log(d)))) + (1 / (1 + np.exp(0.924 - 1.885 * np.log(d))))
+        )
+        fractions_et.append(DF_et)
+
+        DF_tb = (0.00352/d) * (np.exp(-0.234*((np.log(d) + 3.40)**2)) + (63.9 * np.exp(-0.819*((np.log(d) - 1.61)**2))))
+        fractions_tb.append(DF_tb)
+
+        DF_al = (0.0155/d) * (np.exp(-0.416*((np.log(d) + 2.84)**2)) + (19.11*np.exp(-0.482 * ((np.log(d) - 1.362)**2))))
+        fractions_al.append(DF_al)
+
+        DF = IF * (0.0587 + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d)))))
+        fractions_df.append(DF)
+
+    ax.plot(diameters, fractions_df, label = 'Total Deposition', color='k')
+    ax.plot(diameters, fractions_et, label = 'Extrathoracic', ls='-.', lw=0.9, color='grey')
+    ax.plot(diameters, fractions_tb, label = 'Tracheobronchial', ls='--', lw=0.9, color='darkgray')
+    ax.plot(diameters, fractions_al, label = 'Alveolar', ls=(0,(1,1)), lw=0.9, color='darkgray')
+    
+    ax.grid(linestyle='--')
+    ax.set_xscale('log')
+    ax.margins(x=0, y=0)
+    plt.legend(bbox_to_anchor=(1.04,1), loc="upper left")
+
+    y_ticks = [0., 0.2, 0.4, 0.6, 0.8, 1]
+    x_ticks = [0.001, 0.01, 0.1, 1, 10, 100]
+    plt.yticks(ticks=y_ticks, labels=[
+               str(i) for i in y_ticks])
+    plt.xticks(ticks=x_ticks, labels=[
+               str(i) for i in x_ticks])
+    plt.xlabel('Particle diameter (μm)', fontsize=14)
+    plt.ylabel('Deposition fraction\nf$_{dep}$', fontsize=14)
+    
+    fig.set_figwidth(10)
+    plt.tight_layout()
+    plt.show()
+
+############ Compare concentration curves ############
+def compare_concentration_curves():
+
+    exp_models=[office_model_no_mask_windows_closed().build_model(size=SAMPLE_SIZE),
+                office_model_no_mask_windows_open_breaks().build_model(size=SAMPLE_SIZE),
+                office_model_no_mask_windows_open_alltimes().build_model(size=SAMPLE_SIZE)]
+
+    labels=['Windows closed', 'Window open during breaks', 'Window open at all times']
+    colors=['tomato', 'lightskyblue', 'limegreen', '#1f77b4', 'seagreen', 'lightskyblue', 'deepskyblue']
+
+    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)
+
+    TIMESTEP = 0.01
+    times=np.arange(start, stop, TIMESTEP)
+    
+    concentrations = [[np.mean(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):
+        ax.plot(times, c, label=label, color=color)
+    
+    ax.legend(loc='upper left')
+    ax.set_ylim(ax.get_ylim()[0], ax.get_ylim()[1] * 1.2)
+    ax.spines["right"].set_visible(False)
+
+    cumulative_doses = [np.cumsum([
+        np.array(model.exposure_between_bounds(float(time1), float(time2))).mean()
+        for time1, time2 in zip(times[:-1], times[1:])
+    ]) for model in exp_models]
+
+    plt.xlabel("Exposure time ($h$)", fontsize=14)
+    plt.ylabel("Mean concentration (virions m$^{-3}$)", fontsize=14)
+
+    ax1 = ax.twinx()
+    for qd, label, color in zip(cumulative_doses, labels, colors):
+        ax1.plot(times[:-1], qd, label='qD - ' + label, color=color, linestyle='dotted')
+    ax1.spines["right"].set_linestyle("--")
+    ax1.spines["right"].set_linestyle((0,(1,5)))
+    ax1.set_ylabel('Mean cumulative dose (virions)', fontsize=14)
+    ax1.set_ylim(ax1.get_ylim()[0], ax1.get_ylim()[1] * 1.2)
+
+    plt.show()
+
+def compare_viruses_vr():
+
+    # Represented as tuples of three numbers on the interval [0, 1] (e.g. (1, 0, 0)) (R, G, B)
+    colors = [(0., 0.5, 0.5), (0, 0, 0.5), (0.5, 0, 0), (0., 0.78, 0.)]
+    colors_violin=['lightsteelblue', 'wheat', 'darkseagreen']
+    colors_violin_lines = ['royalblue', 'orange', 'forestgreen']
+
+    # The colors of the borders surrounding the violin plots
+    border_colors = [(0, 0, 0), (0, 0, 0), (0, 0, 0)]
+    
+    whisker_width = 0.8
+    positions = [1, 2, 3, 12]
+
+    infected_populations = [infected_model('No mask', 'Light activity', activity).build_model(size=SAMPLE_SIZE) for activity in ('Breathing', 'Talking', 'Shouting')]
+    
+    vrs = [np.log10(pop.emission_rate_when_present(cn_B=0.06, cn_L=0.2)) for pop in infected_populations]
+
+    fig, ax = plt.subplots()
+    ax.set_xlim((0, 11))
+
+    parts = ax.violinplot(vrs, quantiles=[(0.05, 0.95) for _ in vrs], showextrema=False)
+    means = [np.log10(np.mean(10 ** vr)) for vr in vrs]
+    ax.hlines(y=means,
+                xmin=[pos - whisker_width / 2 for pos in positions[:3]],
+                xmax=[pos + whisker_width / 2 for pos in positions[:3]],
+                colors=colors_violin_lines,
+                alpha=0.8)
+
+
+    for pc, color, bc in zip(parts['bodies'], colors_violin, border_colors):
+        pc.set_facecolor(color)
+        pc.set_edgecolor(bc)
+        pc.set_alpha(0.5)
+    parts['cquantiles'].set_color([c for c in colors_violin_lines[:3] for _ in range(2)])
+    parts['cquantiles'].set_alpha(1)
+
+    positions=np.linspace(4.5, 11.5, 20)
+
+    ######### SARS-CoV #########
+    lower_bound = [418]
+    higher_bound = [4176]
+
+    for i in range(len(lower_bound)):
+        data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
+        ax.boxplot(np.log10(data), positions=[positions[i]], medianprops=dict(color=colors[3] + (0.5,)),
+                   whiskerprops=dict(color=colors[3] + (0.5,)), boxprops=dict(color=colors[3] + (0.5,)))
+
+    ######### SARS-CoV-2 #########
+    lower_bound = [216, 216, 518, 648, 878, 893, 1670, 1872, 1915, 2002, 2002, 2189, 3341, 9835, 13968, 60667]
+    higher_bound = [2160, 2160, 5184, 6480, 8784, 8928, 16704, 18720, 19152, 20016, 20016, 21888, 33408, 98352, 139680, 606672]
+    
+    for i in range(len(lower_bound)):
+        data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
+        ax.boxplot(np.log10(data), positions=[positions[i+1]], medianprops=dict(color=colors[0]+ (0.5,)),
+                   whiskerprops=dict(color=colors[0]+ (0.5,)), boxprops=dict(color=colors[0]+ (0.5,)))
+
+    ######### Measles #########
+    lower_bound = [259, 8640, 39816, 124416]
+    higher_bound = [2592, 86400, 398160, 1244160]
+    
+    for i in range(len(lower_bound)):
+        data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
+        ax.boxplot(np.log10(data), positions=[positions[i+5]], medianprops=dict(color=colors[1]+ (0.5,)),
+                   whiskerprops=dict(color=colors[1]+ (0.5,)), boxprops=dict(color=colors[1]+ (0.5,)))
+
+    ######### Influenza #########
+    lower_bound = [2, 114, 1138]
+    higher_bound = [16, 1145, 11376]
+    
+    for i in range(len(lower_bound)):
+        data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
+        ax.boxplot(np.log10(data), positions=[positions[i+12]], medianprops=dict(color=colors[2]+ (0.5,)),
+                   whiskerprops=dict(color=colors[2]+ (0.5,)), boxprops=dict(color=colors[2]+ (0.5,)))
+
+    ######### Rhinovirus #########
+    lower_bound = [45]
+    higher_bound = [446]
+    
+    for i in range(len(lower_bound)):
+        data = np.random.uniform(lower_bound[i], higher_bound[i], size=200000)
+        ax.boxplot(np.log10(data), positions=[positions[i+8]], medianprops=dict(color=(0.5, 0.5, 0.5, 0.5, )),
+                   whiskerprops=dict(color=(0.5, 0.5, 0.5, 0.5,)), boxprops=dict(color=(0.5, 0.5, 0.5, 0.5,)))
+    
+    handles = [patches.Patch(edgecolor=c, facecolor='none', label=l)
+               for c, l in zip([p + (0.5,) for p in [(0., 0.78, 0.), (0., 0.5, 0.5), (0, 0, 0.5), (0.5, 0, 0), (0.5, 0.5, 0.5)]],
+                               ('SARS-CoV', 'SARS-CoV-2', 'Measles', 'Influenza', 'Rhinovirus'))]
+    boxplot_legend = plt.legend(handles=handles, loc='lower right')
+
+    ax.annotate("Bus ride", xy=(6, np.log10(5500)), color='k', fontsize=8,
+                xycoords='data',
+                xytext=(-50, 50), textcoords='offset points',
+                arrowprops=dict(arrowstyle="->",
+                                connectionstyle="arc3,rad=-0.2", color='lightgrey'))
+
+    ax.annotate("S V Chorale", xy=(10, np.log10(110000)), color='k', fontsize=8,
+                xycoords='data',
+                xytext=(-50, 40), textcoords='offset points',
+                arrowprops=dict(arrowstyle="->",
+                                connectionstyle="arc3,rad=-0.2", color='lightgrey'))
+    
+    handles = [patches.Patch(color=c, label=l) for c, l in zip([p  for p in colors_violin], ('Breathing', 'Speaking', 'Shouting'))]
+    plt.legend(handles=handles, loc='lower left', bbox_to_anchor=(0.12, 0.))
+    plt.gca().add_artist(boxplot_legend)
+
+    ax.hlines(y=[-2, 0, 2, 4, 6], xmin=ax.get_xlim()[0], xmax=ax.get_xlim()[1], colors=(0.8, 0.8, 0.8), linestyles='--', alpha=0.3)
+    ax.vlines(x=4, ymin=ax.get_ylim()[0], ymax=ax.get_ylim()[1], colors=(0.8, 0.8, 0.8))
+    ax.set_yticks([i for i in range(-4, 7, 2)])
+    ax.set_yticklabels(['$10^{' + str(i) + '}$' for i in range(-4, 7, 2)])
+    ax.set_xticks([2, 7])
+    ax.set_xticklabels(['SARS-CoV-2\n(model)', 'Literature Data\n(recorded outbreaks) '])
+    ax.set_ylabel('Emission rate (virions h$^{-1}$)', fontsize=12)
+
+    plt.tight_layout()
+    plt.show()
+
+######### 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 define_colormap(cns):
+    min_val, max_val = 0.25, 0.85
+    n = 10
+    orig_cmap = plt.cm.Blues
+    colors = orig_cmap(np.linspace(min_val, max_val, n))
+
+    norm = mpl.colors.Normalize(vmin=cns.min(), vmax=cns.max())
+    cmap = mpl.cm.ScalarMappable(
+        norm=norm, cmap=mpl.colors.LinearSegmentedColormap.from_list("mycmap", colors))
+    cmap.set_array([])
+
+    return cmap
+
+
+def scatter_coleman_data(coleman_etal_vl_breathing, coleman_etal_er_breathing):
+    plt.scatter(coleman_etal_vl_breathing,
+                coleman_etal_er_breathing, marker='x', c='orange')
+    x_hull, y_hull = get_enclosure_points(
+        coleman_etal_vl_breathing, coleman_etal_er_breathing)
+    # plot shape
+    plt.fill(x_hull, y_hull, '--', c='orange', alpha=0.2)
+
+
+def scatter_milton_data(milton_vl, milton_er):
+    try:
+        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)
+    except:
+        print("No data for Milton et al")
+
+
+def scatter_yann_data(yann_vl, yann_er):
+    try:
+        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)
+    except:
+        print("No data for Yan et al")
+
+
+def build_talking_legend(fig):
+    result_from_model = mlines.Line2D(
+        [], [], color='blue', marker='_', linestyle='None')
+    coleman = mlines.Line2D([], [], color='orange',
+                            marker='x', linestyle='None')
+
+    title_proxy = Rectangle((0, 0), 0, 0, color='w')
+    titles = ["$\\bf{CARA \, \\it{(SARS-CoV-2)}:}$",
+              "$\\bf{Coleman \, et \, al. \, \\it{(SARS-CoV-2)}:}$"]
+    leg = plt.legend([title_proxy, result_from_model, title_proxy, coleman],
+                     [titles[0], "Results from model", titles[1], "Dataset"])
+
+    # 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_breathing_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], "Results 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 print_er_info(er: np.array, log_er: np.array):
+    """
+    Prints statistical parameters of a given distribution of ER-values
+    """
+    print(f"MEAN of ER = {np.mean(er)}\n"
+          f"MEAN of log ER = {np.mean(log_er)}\n"
+          f"SD of ER = {np.std(er)}\n"
+          f"SD of log ER = {np.std(log_er)}\n"
+          f"Median of ER = {np.quantile(er, 0.5)}\n")
+
+    print(f"Percentiles of ER:")
+    for quantile in (0.01, 0.05, 0.25, 0.50, 0.55, 0.65, 0.75, 0.95, 0.99):
+        print(f"ER_{quantile} = {np.quantile(er, quantile)}")
+
+    return
+
+
+

cara/test_plots.py

@@ -0,0 +1,86 @@
+from numpy.core.function_base import linspace
+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
+import matplotlib as mpl
+from model_scenarios_paper import *
+
+# Used for testing
+
+######### Plot material #########
+
+SAMPLE_SIZE = 50000
+viral_loads = np.linspace(2, 12, 600)
+
+er_means = []
+er_medians = []
+lower_percentiles = []
+upper_percentiles = []
+def exposure_model_from_vl_talking_new_points():
+    fig = plt.figure()
+    ax = fig.add_subplot(1, 1, 1)
+    
+    for vl in tqdm(viral_loads):
+        exposure_mc = talking_exposure_vl(vl)
+        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(
+            1.0)/4
+        er_means.append((10**vl) / 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 4 to have in 15min (quarter of an hour)
+    coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking]
+
+    ax.plot(viral_loads, er_means)
+    ax.set_yscale('log')
+
+    new_datapoints = [
+        10**(a) / b for a, b in zip(coleman_etal_vl_talking, coleman_etal_er_talking_2)]
+    print(new_datapoints)
+
+    ############# Coleman #############
+    plt.scatter(coleman_etal_vl_talking, new_datapoints, marker='x')
+
+    ############# Markers #############
+    markers = [5, 'd', 4]
+
+    ############ Legend ############
+    result_from_model = mlines.Line2D(
+        [], [], color='blue', marker='_', linestyle='None')
+    coleman = mlines.Line2D([], [], color='orange',
+                            marker='x', linestyle='None')
+
+    title_proxy = Rectangle((0, 0), 0, 0, color='w')
+    titles = ["$\\bf{CARA \, \\it{(SARS-CoV-2)}:}$",
+              "$\\bf{Coleman \, et \, al. \, \\it{(SARS-CoV-2)}:}$"]
+    leg = plt.legend([title_proxy, result_from_model, title_proxy, coleman],
+                     [titles[0], "Result from model", titles[1], "Dataset"])
+
+    # 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 talking for 15min',
+              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.ylim([10**0, 10**10])
+    plt.show()
+    

requirements.txt

@@ -84,3 +84,4 @@ webencodings==0.5.1
 websocket-client==1.1.0
 wheel==0.36.2
 widgetsnbextension==3.5.1
+pandas==1.3.2

setup.cfg

@@ -26,6 +26,12 @@ ignore_missing_imports = True
 [mypy-scipy.*]
 ignore_missing_imports = True
 
+[mypy-tqdm.*]
+ignore_missing_imports = True
+
+[mypy-pandas.*]
+ignore_missing_imports = True
+
 [mypy-timezonefinder.*]
 ignore_missing_imports = True