---
title: "What-if: calibrating exposure"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{What-if: calibrating exposure}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 4.5,
  out.width = "100%"
)
```

```{r setup}
library(ambre)
set.seed(2024)
```

```{r libs}
library(dplyr)
library(purrr)
```

Two numbers decide how much of the water each person actually meets:

- the **volume** swallowed or inhaled per exposure event (litres per event), and
- the **frequency** of those events (how many times it happen per year).

Together with the pathogen concentration, these drive the *dose*, and the dose
drives the risk. This vignette is about the exposure side of the model: where
those two numbers come from by default, how to read them, and how to **override**
them to ask "what if?" — including the two most common management measures,
personal protective equipment and reduced-contact irrigation.

If you are new to the pipeline, read `vignette("a-get-started", package = "ambre")`
first. The barrier-side levers (treatment plants, on-field practices, their
log-reductions) are a different story, told in
`vignette("b-initial-vs-new-scenario", package = "ambre")`.

## How exposure is set by default

You never set volume, frequency or concentration by hand in a normal run. When you call
`inflow_concentration()` (which `run_qmra_intial_situation()` and 
`run_qmra_supplementary_process()` call for you), it runs three updates in 
order:

```
inflow_concentration()
  ├─ update_pathogen()    # pathogen concentration
  ├─ update_frequency()   # number of events / year, per path
  └─ update_volume()      # litres per event, per path
```

```{r default-exposure, results = "hide"}
scenario <- create_scenario(
  system.file("input_1culture_2pop.xlsx", package = "ambre")
)
scenario_conc <- inflow_concentration(
  scenario = scenario,
  pathogenName = "Campylobacter jejuni"
)
```

`update_pathogen()` look up **inflow** database. 
`update_frequency()` and `update_volume()` look up the **exposure path** of each
scenario row and copy the path-specific numbers out of the shipped database
`config_ambre$path` into that row's embedded `config$exposure` table. So the 
model does **not** use one generic volume for every path — each crop × population
× path combination gets its own exposure profile.

You can see the result. The `config$exposure` table has three rows —
`number_of_repeatings` (the Monte-Carlo count), `number_of_exposures` (events per
year) and `volume_perEvent`:

```{r default-exposure-show}
# Row 1 of this two-row scenario: irrigation-staff droplet ingestion
scenario_conc$config[[1]]$exposure
```

Two quirks are worth knowing up front, and we
return to them at the end:

- `update_frequency()` writes only the path's **`max`** frequency, and it writes
  it into the `number_of_exposures` slot — the annual event count *is* the
  frequency.
- both updates use **positional** writes (`min[[3]]`, `max[[3]]`, `value[[2]]`)
  that assume the exposure table keeps its current row order.

## Inspect the base numbers

Before changing anything, look at what the model is assuming. 

Start from the PathogenName, you can see the default value of concentration for the simulated pathogen : 

```{r default-concentration-show}
# Row 1 of this two-row scenario: Campylobacter jejuni concentration
dplyr::filter(scenario_conc$config[[1]]$inflow, PathogenName == "Campylobacter jejuni")
```

Then start from a path
description and get its `PathID` with `query_exp_path()`:

```{r inspect-pathid}
pid <- query_exp_path(
  pathName = "Ingestion of water droplets during maintenance of the irrigation system"
)
pid
```

Then read the per-path volume and frequency straight from the database. Both
helpers accept a vector of IDs, so you can inspect a whole scenario at once:

```{r inspect-values}
scenario$PathID          # the two paths in this input file

query_volume(scenario$PathID)     # litres per event: min / max
query_frequency(scenario$PathID)  # events per year: min / max
```

These are exactly the numbers `update_volume()` and `update_frequency()` inject.
The volumes here are tiny — `0.001` litre is one millilitre of incidentally
swallowed droplets — the frequency `40`–`60` becomes `60` because only the
maximum is used and the concentration is between 100 and 500 log/L.

## Override for a what-if

To explore sensitivity you replace those numbers with your own. Two helpers do
this, and both expect **one value per scenario row**.

`update_volume_with_desired_value()` takes a `volume` argument that is a
data frame (or list) with a `min`, a `max` and a `type` column. Use `min == max` for a
fixed volume, or `min < max` for a distribution spread according to the `type` specified:

```{r override-volume}
scenario_ppe <- update_volume_with_desired_value(
  scenario = scenario,
  volume = data.frame(
    min = c(0.0005, 0.0005),   # 0.5 mL per event, one entry per row
    max = c(0.0005, 0.0005),
    type = c("triangle", "triangle")
  )
)

scenario_ppe$config[[1]]$exposure
```

`update_frequency_with_desired_value()` takes a plain numeric **vector** of 
**integer**, one number of events per year per row:

```{r override-frequency}
scenario_night <- update_frequency_with_desired_value(
  scenario = scenario,
  frequency = c(30L, 30L)   # cap both paths at 30 events / year
)

scenario_night$config[[1]]$exposure
```

The `volume_perEvent` row now reads `0.0005`, and `number_of_exposures` reads
`30`: your what-if values have replaced the ones the database would have injected.


`update_concentration()` takes a **data.frame** with **PathogenName**, **min**, **max** and **type** of ditribution law values.

```{r override-concentration}
concentration_custom <- data.frame(PathogenName = c("Campylobacter jejuni"),
                                   min = c(1),
                                   max = c(2),
                                   type = c("uniform"))

scenario_pathogen <- update_pathogen(scenario = scenario, pathoName = concentration_custom$PathogenName)
scenario_low_pathogen <- update_concentration(scenario = scenario_pathogen ,
                     concentration = concentration_custom)

dplyr::filter(scenario_low_pathogen$config[[1]]$inflow, PathogenName == "Campylobacter jejuni")

```

## Modelling management measures

Read those two overrides as **interventions**. This is the exposure-side
complement to the barrier log-reductions in
`vignette("b-treatment-vs-multibarrier", package = "ambre")`: instead of removing
pathogens from the water, you change how much of the water each person meets.

- **Personal protective equipment (gloves, masks, hygiene rules)** lowers the
  **volume** ingested or inhaled per event. Model it by setting
  `volume_perEvent` below the default — a gloved worker who swallows a fraction
  of a millilitre instead of several.
- **Night-time irrigation, signage, fencing, keeping workers off wet fields**
  lowers the **frequency** of contact. Model it by setting `number_of_exposures`
  to a smaller annual count. Because the annual risk aggregates as
  $1 - \prod(1 - p)$ over events, fewer events mean lower annual risk even when
  each event is unchanged.

Both are assumptions *you* choose and should defend — they are not credited from
a barrier database. (The `barriere_path` / `barriere_specific` tables in
`config_ambre` that *would* encode route-specific exposure reductions are shipped
data, not yet wired into the engine.)

## A mini sensitivity study

Let us quantify how much the volume assumption matters. One catch first:
`run_qmra_initial_situation()` and `run_qmra_supplementary_process()`
re-derives the per-path volumes from the database on every
call (via `inflow_concentration()`), so it would overwrite any override you set.
To inject a what-if value you should use `run_qmra_custom` function:

Run it twice on the same file and pathogen — once with a pessimistic bare-hand
volume, once with a PPE volume ten times smaller:

```{r qmra-custom, results = "hide"}
sc <- create_scenario(system.file("input_1culture_2pop.xlsx", package = "ambre"))

regulation_reduction <- config_ambre$regulation$regulation_value |> 
  dplyr::filter(Country == "France") |>
  dplyr::select(-c(Concentration, Country, RegulationID))
regulation_concentration <- config_ambre$regulation$regulation_value |> 
  dplyr::filter(Country == "France") |>
  dplyr::select(-c(Country, RegulationID, Reduction))

concentration_custom <- data.frame(PathogenName = c("Rotavirus"),
                                   min = c(1000),
                                   max = c(2000),
                                   type = c("uniform"))

set.seed(2024)
bare <- run_qmra_custom(scenario = sc,
                  concentration = concentration_custom,
                  volume = data.frame(min=c(0.005,0.005), max = c(0.01, 0.01), type = c("triangle", "triangle")),
                  frequency = c(48L, 60L),
                  regulationLog = regulation_reduction,
                  regulationConcentration = regulation_concentration,
                  initialSituation = TRUE) # min 5 mL, max 10 mL

set.seed(2024)
ppe  <- run_qmra_custom(scenario = sc,
                  concentratio = concentration_custom,
                  volume = data.frame(min=c(0.0005,0.0005), max = c(0.001, 0.001), type = c("triangle", "triangle")),
                  frequency = c(48L, 60L),
                  regulationLog = regulation_reduction,
                  regulationConcentration = regulation_concentration,
                  initialSituation = TRUE) # min 0.5 mL min 0.1 mL

```

As the the other functions `run_qmra_*` this custom function return 2 plots and to formattable. See `d-interpreting-risk`for more detail on the output.

```{r qmra-custom-comparison, results = "hide"}
library(ggplot2)
cowplot::plot_grid(
      bare$dalys$Rotavirus +
        labs(subtitle = "min 5 mL, max 10 mL") +
        theme(plot.subtitle = element_text(hjust = 0.5)),
      ppe$dalys$Rotavirus +
        labs(subtitle = "min 0.5 mL min 0.1 mL") +
        theme(plot.subtitle = element_text(hjust = 0.5)),
      ncol = 2,
      align = "h"
    )
```


The risk tracks the volume almost proportionally: a tenfold cut in swallowed
volume gives roughly a tenfold cut in DALYs. More tellingly, the 95th percentile
crosses the line — the bare-hand upper tail sits **above** 1e-6 while the PPE tail
drops **below** it. Deciding whether a scenario meets the target can come down to
this single exposure assumption. For how to read these ranges against the
regulatory line, see
`vignette("d-interpreting-risk", package = "ambre")`; for why every number is a
distribution rather than a point, see
`vignette("g-monte-carlo-engine", package = "ambre")`.

## Technical notes and caveats

- **Positional writes are fragile.** `update_volume()` writes `config$exposure$min[[3]]`
  and `max[[3]]`; `update_frequency()` writes `config$exposure$value[[2]]`; the
  `_with_desired_value()` helpers do the same. The indices assume `volume_perEvent`
  is the third row and `number_of_exposures` the second. If the exposure table's
  row order ever changed, these would silently target the wrong parameter — so
  always re-inspect `config$exposure` after an override, as we did above.
- **Frequency uses the maximum only.** `update_frequency()` copies the path's
  `max` frequency (ignoring `min`) into `number_of_exposures`. The default run is
  therefore conservative on contact frequency. Your own
  `update_frequency_with_desired_value()` value is used exactly as given.
- **Volume units.** `volume_perEvent` multiplies the inflow concentration
  (organisms per litre) to yield organisms per event, so it is a volume in
  **litres**. The per-path values (roughly `1e-4` to `1.5e-2` L) are fractions of
  a millilitre up to a few millilitres of incidentally ingested water — or, for
  consumption paths, an equivalent mass of product treated as a volume. The CSVs
  carry no explicit unit column, so keep your overrides on the same litre basis.
- **Overrides and the one-call entry points.** Because `run_qmra_initial_situation()` and
  `run_qmra_supplementary_process()` rebuild the exposure profile from the database each time,
  the `_with_desired_value()` helpers only take effect in a hand-assembled chain
  like the one above. There is no exposure-override argument on the one-call
  functions.