Removed script files

# 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/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) +
+            return (self.BLO_factors[0] * _Bmode(d, cn_B) +
-                    self.BLO_factors[1] * _Lmode(d) +
+                    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
+    #: The percentage of host immunity
-    def emission_rate_when_present(self) -> _VectorisedFloat:
+    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,
                 ),
 }
 
@@ -67,4 +76,6 @@ virus_distributions = {
 mask_distributions = {
     'Type I': mc.Mask(Uniform(0.25, 0.80)),
     'FFP2': mc.Mask(Uniform(0.83, 0.91)),
-}
+}
+
+

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