UI adjustments on the short range interactions

Removed unused code

cara/apps/calculator/model_generator.py

@@ -13,8 +13,8 @@ from cara import data
 import cara.data.weather
 import cara.monte_carlo as mc
 from .. import calculator
-from cara.monte_carlo.data import activity_distributions, dilution_factor, virus_distributions, mask_distributions, initial_concentrations_mouth
-from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions
+from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions, dilution_factor
+from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions, short_range_expiration_distributions
 
 
 LOG = logging.getLogger(__name__)
@@ -234,7 +234,7 @@ class FormData:
             raise ValueError("mechanical_ventilation_type cannot be 'not-applicable' if "
                              "ventilation_type is 'mechanical_ventilation'")
 
-    def build_mc_model(self) -> mc.ExposureModel:
+    def build_mc_model(self) -> mc.SimulationModel:
         # Initializes room with volume either given directly or as product of area and height
         if self.volume_type == 'room_volume_explicit':
             volume = self.room_volume
@@ -253,28 +253,29 @@ class FormData:
             sr_presence=self.short_range_intervals()
             sr_activities=self.short_range_activities()
 
-        jet_origin_concentration = [initial_concentrations_mouth[activity] for activity in sr_activities]
-        
-        short_range = models.ShortRangeModel(
-            presence=sr_presence,
-            activities=sr_activities,
-            dilutions=dilution_factor(activities=sr_activities, distance=np.random.uniform(0.5, 1.5, 250000)),
-            jet_origin_concentrations=jet_origin_concentration,
-        )
+        short_range_expirations = [short_range_expiration_distributions[activity] for activity in sr_activities]
 
         # Initializes and returns a model with the attributes defined above
-        return mc.ExposureModel(
-            concentration_model=mc.ConcentrationModel(
-                room=room,
-                ventilation=self.ventilation(),
-                infected=self.infected_population(),
-                short_range=short_range,
-                evaporation_factor=0.3,
+        return mc.SimulationModel( 
+            mc.ExposureModel(
+                concentration_model=mc.ConcentrationModel(
+                    room=room,
+                    ventilation=self.ventilation(),
+                    infected=self.infected_population(),
+                    evaporation_factor=0.3,
+                ),
+                exposed=self.exposed_population(),
+            ),
+            mc.ShortRangeModel(
+                presence=sr_presence,
+                long_range_presence=self.long_range_intervals(),
+                activities=sr_activities,
+                expirations=short_range_expirations,
+                dilutions=dilution_factor(activities=sr_activities, distance=np.random.uniform(0.5, 1.5, 250000)),
             ),
-            exposed=self.exposed_population(),
         )
 
-    def build_model(self, sample_size=_DEFAULT_MC_SAMPLE_SIZE) -> models.ExposureModel:
+    def build_model(self, sample_size=_DEFAULT_MC_SAMPLE_SIZE) -> models.SimulationModel:
         return self.build_mc_model().build_model(size=sample_size)
 
     def tz_name_and_utc_offset(self) -> typing.Tuple[str, float]:
@@ -643,7 +644,7 @@ class FormData:
     def infected_present_interval(self) -> models.Interval:
         return self.present_interval(
             self.infected_start, self.infected_finish,
-            breaks=self.infected_lunch_break_times() + self.infected_coffee_break_times(),
+            breaks=self.infected_lunch_break_times() + self.infected_coffee_break_times()
         )
 
     def short_range_intervals(self) -> typing.List[models.SpecificInterval]:
@@ -657,10 +658,22 @@ class FormData:
         else:
             return []
 
+    def long_range_intervals(self) -> models.Interval:
+        short_range_intervals = []
+        for interval in self.short_range_interactions:
+            short_range_intervals.append(
+                (time_string_to_minutes(interval['start_time']), 
+                time_string_to_minutes(interval['start_time']) + float(interval['duration'])),
+            )
+        return self.present_interval(
+            self.exposed_start, self.exposed_finish,
+            breaks=self.infected_lunch_break_times() + self.infected_coffee_break_times() + tuple(short_range_intervals),
+        )
+
     def exposed_present_interval(self) -> models.Interval:
         return self.present_interval(
             self.exposed_start, self.exposed_finish,
-            breaks=self.exposed_lunch_break_times() + self.exposed_coffee_break_times(),
+            breaks=self.exposed_lunch_break_times() + self.exposed_coffee_break_times()
         )
 
     def short_range_activities(self) -> typing.List[str]:
@@ -678,12 +691,6 @@ def build_expiration(expiration_definition) -> models._ExpirationBase:
             for exp_type, weight in expiration_definition.items()
             ], axis=0)
         return expiration_distribution(tuple(BLO_factors))
-
-
-def build_concentration_mouth(short_range_def) -> models._ExpirationBase:
-    expiration_definition = short_range_def[1]['activity']
-    if isinstance(expiration_definition, str):
-        return initial_concentrations_mouth[expiration_definition]
     
 
 def baseline_raw_form_data():

cara/apps/calculator/report_generator.py

@@ -96,18 +96,29 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L
     return nice_times
 
 
-def calculate_report_data(model: models.ExposureModel):
-    times = interesting_times(model)
+def short_range_interesting_times(model: models.ExposureModel, times: typing.List[float]) -> typing.List[float]:
+    short_range_times : typing.List[float] = []
+    for period in model.concentration_model.infected.short_range_presence:
+        start, finish = tuple(period.boundaries())
+        short_range_times = short_range_times + [time for time in times if time >= start and time <= finish]
+    return short_range_times
+
+
+def calculate_report_data(model: models.SimulationModel):
+    times = interesting_times(model.exposure_model)
+    short_range_intervals = []
+    for interval in model.short_range.presence:
+        short_range_intervals.append(list(interval.boundaries()))
 
     concentrations = [
-        np.array(model.concentration_model.concentration(float(time))).mean()
+        np.array(model.concentration(float(time))).mean()
         for time in times
     ]
     highest_const = max(concentrations)
     prob = np.array(model.infection_probability()).mean()
-    er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean()
-    exposed_occupants = model.exposed.number
-    expected_new_cases = np.array(model.expected_new_cases()).mean()
+    er = np.array(model.exposure_model.concentration_model.infected.emission_rate_when_present()).mean()
+    exposed_occupants = model.exposure_model.exposed.number
+    expected_new_cases = np.array(model.exposure_model.expected_new_cases()).mean()
     cumulative_doses = np.cumsum([
         np.array(model.deposited_exposure_between_bounds(float(time1), float(time2))).mean()
         for time1, time2 in zip(times[:-1], times[1:])
@@ -115,7 +126,8 @@ def calculate_report_data(model: models.ExposureModel):
 
     return {
         "times": list(times),
-        "exposed_presence_intervals": [list(interval) for interval in model.exposed.presence.boundaries()],
+        "short_range_intervals": short_range_intervals,
+        "exposed_presence_intervals": [list(interval) for interval in model.exposure_model.exposed.presence.boundaries()],
         "cumulative_doses": list(cumulative_doses),
         "concentrations": concentrations,
         "highest_const": highest_const,
@@ -234,20 +246,20 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp
     return scenarios
 
 
-def scenario_statistics(mc_model: mc.ExposureModel, sample_times: np.ndarray):
+def scenario_statistics(mc_model: mc.SimulationModel, sample_times: np.ndarray):
     model = mc_model.build_model(size=_DEFAULT_MC_SAMPLE_SIZE)
     return {
-        'probability_of_infection': np.mean(model.infection_probability()),
-        'expected_new_cases': np.mean(model.expected_new_cases()),
+        'probability_of_infection': np.mean(model.exposure_model.infection_probability()),
+        'expected_new_cases': np.mean(model.exposure_model.expected_new_cases()),
         'concentrations': [
-            np.mean(model.concentration_model.concentration(time))
+            np.mean(model.exposure_model.concentration_model.concentration(time))
             for time in sample_times
         ],
     }
 
 
 def comparison_report(
-        scenarios: typing.Dict[str, mc.ExposureModel],
+        scenarios: typing.Dict[str, mc.SimulationModel],
         sample_times: typing.List[float],
         executor_factory: typing.Callable[[], concurrent.futures.Executor],
 ):
@@ -285,7 +297,7 @@ class ReportGenerator:
     def prepare_context(
             self,
             base_url: str,
-            model: models.ExposureModel,
+            model: models.SimulationModel,
             form: FormData,
             executor_factory: typing.Callable[[], concurrent.futures.Executor],
     ) -> dict:
@@ -293,12 +305,12 @@ class ReportGenerator:
         time = now.strftime("%Y-%m-%d %H:%M:%S UTC")
 
         context = {
-            'model': model,
+            'model': model.exposure_model,
             'form': form,
             'creation_date': time,
         }
 
-        scenario_sample_times = interesting_times(model)
+        scenario_sample_times = interesting_times(model.exposure_model)
 
         context.update(calculate_report_data(model))
         alternative_scenarios = manufacture_alternative_scenarios(form)

cara/apps/calculator/static/js/form.js

@@ -404,6 +404,10 @@ function validate_form(form) {
       short_range_interactions.push(JSON.stringify(obj));
   });
   $("input[type=text][name=short_range_interactions]").val('[' + short_range_interactions + ']');
+  if (short_range_interactions.length == 0) {
+    $("input[type=radio][id=short_range_no]").prop("checked", true);
+    on_short_range_option_change();
+  }
 
   if (submit) {
     $("#generate_report").prop("disabled", true);
@@ -859,7 +863,11 @@ $(document).ready(function () {
   //Short range modal - save button
   $("body").on("click", ".save_btn_frm_field", function() {
     var last_element = $(".form_field_outer").find(".form_field_outer_row").last().find(".short_range_option").prop("id");
-    if (!last_element) $('#short_range_dialog').modal('hide');
+    if (!last_element) {
+      $('#short_range_dialog').modal('hide');
+      $("input[type=radio][id=short_range_no]").prop("checked", true);
+      on_short_range_option_change();
+    } 
     else {
       let index = last_element.split("_").slice(-1)[0];
       let activity = validate_sr_parameter('#sr_activity_no_' + String(index)[0], "You must specify the activity type.");
@@ -875,13 +883,14 @@ $(document).ready(function () {
     }
   });
 
-  //Short range modal - close button
+  //Short range modal - reset button
   $("body").on("click", ".dismiss_btn_frm_field", function() {
     $(".form_field_outer_row").remove();
     $("#sr_interactions").text(0);
+    $('input[type=radio][id=short_range_no]').prop("checked", true);
+    on_short_range_option_change();
   });
 
-
 });
 
 

cara/apps/calculator/static/js/report.js

@@ -1,5 +1,5 @@
 /* Generate the concentration plot using d3 library. */
-function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses, exposed_presence_intervals) {
+function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses, exposed_presence_intervals, short_range_intervals) {
 
     var time_format = d3.timeFormat('%H:%M');
 
@@ -50,6 +50,16 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
             .attr('fill-opacity', '0.1');
     });
 
+    // Area representing the short range interaction(s).
+    var shortRangeArea = {};
+    var drawShortRangeArea = {};
+    short_range_intervals.forEach((b, index) => {
+        shortRangeArea[index] = d3.area();
+        drawShortRangeArea[index] = vis.append('svg:path')
+            .attr('fill', '#1f00b4')
+            .attr('fill-opacity', '0.1');
+    });
+
     // Plot tittle.
     var plotTitleEl = vis.append('svg:foreignObject')
         .attr("background-color", "transparent")
@@ -105,10 +115,16 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
 
     var legendAreaIcon = vis.append('rect')
         .attr('width', 20)
-        .attr('height', 20)
+        .attr('height', 15)
         .attr('fill', '#1f77b4')
         .attr('fill-opacity', '0.1');
 
+    var legendShortRangeAreaIcon = vis.append('rect')
+        .attr('width', 20)
+        .attr('height', 15)
+        .attr('fill', '#1f00b4')
+        .attr('fill-opacity', '0.1');
+
     var legendLineText = vis.append('text')
         .text('Mean concentration')
         .style('font-size', '15px')
@@ -124,10 +140,15 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
         .style('font-size', '15px')
         .attr('alignment-baseline', 'central');
 
+    var legendShortRangeText = vis.append('text')
+        .text('Short range interaction(s)')
+        .style('font-size', '15px')
+        .attr('alignment-baseline', 'central');
+
     // Legend bounding
     var legendBBox = vis.append('rect')
         .attr('width', 255)
-        .attr('height', 70)
+        .attr('height', 90)
         .attr('stroke', 'lightgrey')
         .attr('stroke-width', '2')
         .attr('rx', '5px')
@@ -217,7 +238,7 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
             .y(d => yRange(d.concentration));
         draw_line.attr("d", lineFunc(data_for_graphs.concentrations));
 
-        // Cumulative line
+        // Cumulative line.
         lineCumulative.defined(d => !isNaN(d.concentration))
             .x(d => xTimeRange(d.time))
             .y(d => yCumulativeRange(d.concentration));
@@ -234,6 +255,18 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
             })));
         });
 
+        // Short Range Area.
+        short_range_intervals.forEach((b, index) => {
+            shortRangeArea[index].x(d => xTimeRange(d.time))
+                .y0(graph_height - margins.bottom)
+                .y1(d => yRange(d.concentration));
+
+            drawShortRangeArea[index].attr('d', shortRangeArea[index](data_for_graphs.concentrations.filter(d => {
+                return d.time >= b[0] && d.time <= b[1]
+            })))
+        })
+
+
         // Title.
         plotTitleEl.attr('width', graph_width);
 
@@ -282,10 +315,15 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
                 .attr('y', margins.top + 2 * size);
             
             legendAreaIcon.attr('x', graph_width + size * 2.5)
-                .attr('y', margins.top + 2.5 * size);
+                .attr('y', margins.top + 2.6 * size);
             legendAreaText.attr('x', graph_width + 4 * size)
                 .attr('y', margins.top + 3 * size);
             
+            legendShortRangeAreaIcon.attr('x', graph_width + size * 2.5)
+                .attr('y', margins.top + 3.6 * size);
+            legendShortRangeText.attr('x', graph_width + 4 * size)
+                .attr('y', margins.top + 4 * size);
+            
             legendBBox.attr('x', graph_width * 1.07)
                 .attr('y', margins.top * 1.2);
         }
@@ -306,7 +344,12 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
             legendAreaIcon.attr('x', size * 0.50)
                 .attr('y', graph_height * 1.09 + size);
             legendAreaText.attr('x', 2 * size)
-                .attr('y', graph_height + 2.7 * size);
+                .attr('y', graph_height + 2.6 * size);
+
+            legendShortRangeAreaIcon.attr('x', size * 0.50)
+                .attr('y', graph_height * 1.175 + size);
+            legendShortRangeText.attr('x', 2 * size)
+                .attr('y', graph_height + 3.65 * size);
 
             legendBBox.attr('x', 1)
                 .attr('y', graph_height);

cara/apps/expert.py

@@ -501,7 +501,6 @@ baseline_model = models.ExposureModel(
             expiration=models.Expiration.types['Speaking'],
             host_immunity=0.,
         ),
-        short_range=[],
         evaporation_factor=0.3,
     ),
     exposed=models.Population(

cara/apps/templates/base/calculator.report.html.j2

@@ -96,7 +96,8 @@
 									var concentrations = {{ concentrations | JSONify }}
 									var cumulative_doses = {{ cumulative_doses | JSONify }}
 									var exposed_presence_intervals = {{ exposed_presence_intervals | JSONify }}
-									draw_concentration_plot("concentration_plot", times, concentrations, cumulative_doses, exposed_presence_intervals);
+									var short_range_intervals = {{ short_range_intervals | JSONify }}
+									draw_concentration_plot("concentration_plot", times, concentrations, cumulative_doses, exposed_presence_intervals, short_range_intervals);
 								</script>
               				</p>
 						</div>
@@ -309,6 +310,18 @@
 							Gym = For comparison only, all persons doing heavy physical exercise, breathing and not speaking.
 							{% endif %}
 							</p></li>
+							{% if form.short_range_option == "short_range_yes" %}
+							<li><p class="data_text">
+							Short range interactions: {{ form.short_range_interactions|length }}
+							</p></li>
+							<ul>
+							{% for interaction in form.short_range_interactions %}
+								<li>Expiratory activity {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.activity }} </li>
+								<li>Start time {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.start_time }} </li>
+								<li>Duration {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.duration }} {{ "minutes" if interaction.duration|int > 1 else "minute" }}</li>
+							{% endfor %}
+							</ul>
+							{% endif %}
 							<li><p class="data_text">Exposed occupant(s) activity time:</p></li>
 								<ul>
 									<li><p class="data_subtext">Start time: {{ form.exposed_start | minutes_to_time }}</p></li>

cara/models.py

@@ -667,6 +667,11 @@ class Expiration(_ExpirationBase):
         return self.cn * (volume(self.diameter) *
                 (1 - mask.exhale_efficiency(self.diameter))) * 1e-12
 
+    def jet_origin_concentration(self):
+        def volume(d):
+            return (np.pi * d**3) / 6.
+        return self.cn * volume(self.diameter)
+
 
 @dataclass(frozen=True)
 class MultipleExpiration(_ExpirationBase):
@@ -885,24 +890,22 @@ class InfectedPopulation(_PopulationWithVirus):
 @dataclass(frozen=True)
 class ShortRangeModel:
     #: Short range interactions
-    presence: typing.List[SpecificInterval]
+    presence: typing.List[Interval]
+
+    #: Long range presence intervals
+    long_range_presence: Interval
 
     #: The type of expiractory activities in the short range interactions
     activities: typing.List[str]
 
+    #: Expiration type
+    expirations: typing.List[Expiration]
+
     #: The dilution factors for each of the expiratory activity in the short range interactions
     dilutions: _VectorisedFloat
-    
-    #: The concentration on the jet origin for each of the expiratory activity in the short range interactions
-    jet_origin_concentrations: _VectorisedFloat
 
-    def short_range_concentration(self, time: float, background_concentration: _VectorisedFloat, viral_load: _VectorisedFloat) -> _VectorisedFloat:
-        for index, period in enumerate(self.presence):
-            start, finish = tuple(period.boundaries())
-            if time >= start and time <= finish:
-                dilution = self.dilutions[index]
-                jet_origin_concentration = self.jet_origin_concentrations[index] * 1e-6 * viral_load
-                return background_concentration + ((1/dilution)*(jet_origin_concentration - background_concentration))
+    def short_range_concentration(self, long_range_concentration: _VectorisedFloat, dilution: _VectorisedFloat, jet_origin_concentration: _VectorisedFloat) -> _VectorisedFloat: 
+        return long_range_concentration + ((1/dilution)*(jet_origin_concentration - long_range_concentration))
 
 
 @dataclass(frozen=True)
@@ -910,7 +913,6 @@ class ConcentrationModel:
     room: Room
     ventilation: _VentilationBase
     infected: _PopulationWithVirus
-    short_range: ShortRangeModel
 
     #: evaporation factor: the particles' diameter is multiplied by this
     # factor as soon as they are in the air (but AFTER going out of the,
@@ -1042,13 +1044,8 @@ class ConcentrationModel:
         Note that time is not vectorised. You can only pass a single float
         to this method.
         """
-        background_concentration = self._normed_concentration(time) * self.infected.emission_rate_when_present()
-        for period in self.short_range.presence:
-            start, finish = tuple(period.boundaries())
-            if time >= start and time <= finish:
-                return self.short_range.short_range_concentration(time, background_concentration, self.virus.viral_load_in_sputum)
-
-        return background_concentration
+        return (self._normed_concentration(time) * 
+            self.infected.emission_rate_when_present())
 
     @method_cache
     def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat:
@@ -1056,11 +1053,6 @@ class ConcentrationModel:
         Get the integrated concentration of viruses in the air  between the times start and stop,
         normalized by the emission rate.
         """
-        for period in self.short_range.presence:
-            time1, time2 = tuple(period.boundaries())
-            if (time1 <= start <= time2 and time1 <= stop <= time2): 
-                # Check if the given times are within the short range interactions
-                return scipy.integrate.quad_vec(lambda t: self.concentration(t), start, stop)[0]
 
         if stop <= self._first_presence_time():
             return 0.0
@@ -1226,3 +1218,111 @@ class ExposureModel:
         )
 
         return single_exposure_model.expected_new_cases()
+
+
+@dataclass(frozen=True)
+class SimulationModel:
+    exposure_model: ExposureModel
+    short_range: ShortRangeModel
+
+    def normed_integrated_concentration(self, time1, time2):
+        return scipy.integrate.quad_vec(lambda t: self._normed_concentration(t), time1, time2)[0]
+
+    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"""
+        exposure = 0.
+        for start, stop in self.exposure_model.exposed.presence.boundaries():
+            if stop < time1:
+                continue
+            elif start > time2:
+                break
+            elif start <= time1 and time2<= stop:
+                exposure += self.normed_integrated_concentration(time1, time2)
+            elif start <= time1 and stop < time2:
+                exposure += self.normed_integrated_concentration(time1, stop)
+            elif time1 < start and time2 <= stop:
+                exposure += self.normed_integrated_concentration(start, time2)
+            elif time1 <= start and stop < time2:
+                exposure += self.normed_integrated_concentration(start, stop)
+        return exposure
+
+    def _normed_concentration(self, time: float) -> _VectorisedFloat:
+        for index, period in enumerate(self.short_range.presence):
+            start, finish = tuple(period.boundaries())
+            if start <= time <= finish:
+                model = nested_replace(
+                    self, {'exposure_model.concentration_model.infected.expiration': self.short_range.expirations[index]}
+                )
+                dilution = self.short_range.dilutions[index]
+                jet_origin_concentration = (model.exposure_model.concentration_model.infected.expiration.jet_origin_concentration()
+                                            * 1e-6 * model.exposure_model.concentration_model.virus.viral_load_in_sputum)
+                long_range_normed_concentration=np.interp(model.short_range.expirations[index].particle.diameter, self.exposure_model.concentration_model.infected.particle.diameter, self.exposure_model.concentration_model._normed_concentration(time))
+                return self.short_range.short_range_concentration(long_range_normed_concentration, dilution, jet_origin_concentration) 
+        
+        return self.exposure_model.concentration_model._normed_concentration(time)
+
+    def concentration(self, time: float) -> _VectorisedFloat:
+        for index, period in enumerate(self.short_range.presence):
+            start, finish = tuple(period.boundaries())
+            if start <= time <= finish:
+                model = nested_replace(
+                    self, {'exposure_model.concentration_model.infected.expiration': self.short_range.expirations[index]}
+                )
+                return model._normed_concentration(time) * model.exposure_model.concentration_model.infected.emission_rate_when_present()
+        return self._normed_concentration(time) * self.exposure_model.concentration_model.infected.emission_rate_when_present()
+
+    def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
+        for index, period in enumerate(self.short_range.presence):
+            start, finish = tuple(period.boundaries())
+            if (start <= time1 <= finish and start <= time2 <= finish): # What if one is SR and another LR?
+                # Check if the given times are within the short range interactions
+                model = nested_replace(
+                    self, {'exposure_model.concentration_model.infected.expiration': self.short_range.expirations[index]}
+                )
+                emission_rate_per_aerosol = model.exposure_model.concentration_model.infected.emission_rate_per_aerosol_when_present()
+                aerosols = model.exposure_model.concentration_model.infected.aerosols()
+                fdep = model.exposure_model.fraction_deposited()
+                f_inf = model.exposure_model.concentration_model.infected.fraction_of_infectious_virus()
+
+                diameter = self.exposure_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(self._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 = self._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.exposure_model.exposed.activity.inhalation_rate * 
+                        (1 - self.exposure_model.exposed.mask.inhale_efficiency()))
+
+        return self.exposure_model.deposited_exposure_between_bounds(time1, time2)
+
+    def infection_probability(self):
+        dose = 0.0
+        for start, stop in self.short_range.long_range_presence.boundaries():
+            dose += self.exposure_model.deposited_exposure_between_bounds(start, stop)
+
+        for presence in self.short_range.presence:
+            start, stop = presence.boundaries()
+            dose += self.deposited_exposure_between_bounds(start, stop)
+
+        vD = dose * self.exposure_model.repeats
+
+        # oneoverln2 multiplied by ID_50 corresponds to ID_63.
+        infectious_dose = oneoverln2 * self.exposure_model.concentration_model.virus.infectious_dose
+
+        # Probability of infection.        
+        return (1 - np.exp(-((vD * (1 - self.exposure_model.exposed.host_immunity))/(infectious_dose * 
+                self.exposure_model.concentration_model.virus.transmissibility_factor)))) * 100
+        

cara/monte_carlo/data.py

@@ -161,7 +161,7 @@ mask_distributions = {
 }
 
 
-def expiration_distribution(BLO_factors):
+def expiration_distribution(BLO_factors, d_max=30.):
     """
     Returns an Expiration with an aerosol diameter distribution, defined
     by the BLO factors (a length-3 tuple).
@@ -170,10 +170,10 @@ def expiration_distribution(BLO_factors):
     an historical choice based on previous implementations of the model
     (it limits the influence of the O-mode).
     """
-    dscan = np.linspace(0.1, 30. ,3000)
+    dscan = np.linspace(0.1, d_max ,3000)
     return mc.Expiration(CustomKernel(dscan,
                 BLOmodel(BLO_factors).distribution(dscan),kernel_bandwidth=0.1),
-                cn=BLOmodel(BLO_factors).integrate(0.1, 30.))
+                cn=BLOmodel(BLO_factors).integrate(0.1, d_max))
 
 
 def dilution_factor(activities, distance, D=0.02):
@@ -199,11 +199,6 @@ def dilution_factor(activities, distance, D=0.02):
     return factors
 
 
-def initial_concentration_mouth(BLO_factors):
-    value, error = scipy.integrate.quad(lambda d: BLOmodel(BLO_factors).distribution(d) * BLOmodel(BLO_factors).volume(d), 0.1, 1000)
-    return value
-
-
 expiration_BLO_factors = {
     'Breathing': (1., 0., 0.),
     'Speaking':   (1., 1., 1.),
@@ -218,7 +213,7 @@ expiration_distributions = {
 }
 
 
-initial_concentrations_mouth = {
-    exp_type: initial_concentration_mouth(BLO_factors)
+short_range_expiration_distributions = {
+    exp_type: expiration_distribution(BLO_factors, d_max=1000).build_model(250000)
     for exp_type,BLO_factors in expiration_BLO_factors.items()
 }