---
title: "The Monte-Carlo engine"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{The Monte-Carlo engine}
  %\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)
```

What separates `ambre` from a spreadsheet is that almost nothing in it is a
single number. Concentrations, exposure volumes and log-reductions are all
**probability distributions**, and the package draws from them thousands of
times to report a *range* of risk rather than one fragile point estimate. This
vignette opens up that Monte-Carlo core: the two loops it runs, the distributions
it understands, and how to reproduce a result. It is the machinery behind the
boxplots in `vignette("a-get-started", package = "ambre")` and the spread you
interpret in `vignette("d-interpreting-risk", package = "ambre")`.

## Two dimensions: repeats and events

Every simulated quantity is drawn on a grid of two indices, both controlled by
`config_ambre$exposure`:

```{r exposure-control}
config_ambre$exposure[, c("name", "type", "value", "min", "max")]
```

- **`number_of_repeatings`** (default **1000**) is the number of Monte-Carlo
  **runs**. Each run is one plausible world; together they form the *uncertainty
  range*. This index is called `repeatID`.
- **`number_of_exposures`** (default **365**) is the number of exposure
  **events** in a year. This index is called `eventID`, and it is what the annual
  risk aggregation `1 - prod(1 - p)` runs over.

So a single draw of, say, exposure volume is really a `repeatID x eventID` table
-- 1000 x 365 values -- and the whole pipeline carries that shape through to the
final DALYs.

## The distribution catalog

The third row above, `volume_perEvent`, is a `triangle` with `min` 0.5 and `max`
3. That `type` / `value` / `min` / `max` / `mode` / `mean` / `sd` / `meanlog` /
`sdlog` block is the **shared vocabulary** of the whole database: the same nine
columns describe an inflow concentration, a treatment log-reduction and a cost.
`create_random_distribution()` reads that vocabulary. The supported `type` values
are:

| `type` | drawn with | parameters used |
|---|---|---|
| `value` | (constant) | `value` |
| `uniform` | `runif` | `min`, `max` |
| `log10_uniform` | `10^runif` | `min`, `max` (as log10) |
| `norm` | `rnorm` | `mean`, `sd` |
| `log10_norm` | `10^rnorm` | `mean`, `sd` (as log10) |
| `lognorm` | `rlnorm` | `meanlog`, `sdlog` |
| `triangle` | `EnvStats::rtri` | `min`, `max`, `mode` |

Call it directly to see what a draw looks like. It returns an `events` table
(one row per `repeatID x eventID`) and the `paras` actually used:

```{r crd}
draw <- create_random_distribution(
  type = "uniform",
  number_of_repeatings = 3,
  number_of_events = 5,
  min = 1, max = 10,
  debug = FALSE
)
head(draw$events)
draw$paras
```

In practice you rarely call it by hand. The pipeline calls
`generate_random_values()`, a thin wrapper that reads one row of a config table
and fills in sensible defaults before delegating. Feed it the real
`volume_perEvent` row and you get the exact draw the engine would use for
exposure volume:

```{r grv}
vol_row <- config_ambre$exposure[config_ambre$exposure$name == "volume_perEvent", ]
gv <- generate_random_values(
  vol_row,
  number_of_repeatings = 2,
  number_of_events = 4,
  debug = FALSE
)
summary(gv$events$volume_perEvent)
```

The distribution schema is described further in
`vignette("h-config-ambre", package = "ambre")`.

## From min/max to a spread

For the normal and log-normal families you often supply only a `min` and a
`max`, not a mean and standard deviation. `ambre` derives the missing parameters
so that a given fraction of the mass falls inside `[min, max]`, using
`get_percentile()` (a wrapper around the normal quantile) together with
`default_min()` / `default_max()`. For example, the z-score that places 90% of a
normal distribution within the min-max window is:

```{r percentile}
ambre:::get_percentile(0.9)
```

You do not normally call these yourself, but knowing they exist explains why a
`norm` row with only `min` and `max` still produces a sensible bell curve.

## A special case worth knowing

`generate_random_values()` has one silent behaviour. A `triangle` whose `min`
equals its `max` cannot be a triangle, so it is quietly converted to a
(degenerate) `uniform` -- effectively a constant:

```{r triangle-collapse}
degenerate <- data.frame(
  type = "triangle", value = NA,
  min = 2, max = 2, mode = NA,
  mean = NA, sd = NA, meanlog = NA, sdlog = NA
)
gv2 <- generate_random_values(degenerate, number_of_repeatings = 1,
                              number_of_events = 3, debug = FALSE)
unique(gv2$events$values)
```

You will see the message *"Distribution set from 'triangle' to 'uniform' because
'min' equals 'max'"* in a real run when this happens -- it is expected, not an
error.

## Why barrier credits add up

Log-reductions combine by **addition in log space**. A train that removes 2 log
at the plant and 3 log in the field removes 5 log overall, and 5 log means the
surviving concentration is `10^-5` of the inflow. Internally `ambre` builds a
wide `repeatID x eventID` table of each barrier's draw and sums the columns, so
the additivity holds *run by run*, preserving the correlation structure across
the Monte-Carlo sample. This is why you can compare a treatment strategy and a
barrier strategy on the same footing (see
`vignette("b-initial-vs-new-scenario", package = "ambre")`).

## Reproducibility

Because every draw is random, set a seed before a run you want to reproduce.
Same seed, same numbers:

```{r seed}
set.seed(1)
a <- create_random_distribution(type = "uniform", number_of_repeatings = 1,
                                number_of_events = 4, min = 0, max = 1, debug = FALSE)
set.seed(1)
b <- create_random_distribution(type = "uniform", number_of_repeatings = 1,
                                number_of_events = 4, min = 0, max = 1, debug = FALSE)
identical(a$events$values, b$events$values)
```

To read the uncertainty a run expresses, summarise a quantity across `repeatID`
with quantiles -- exactly what the final DALY boxplots do, and what
`vignette("d-interpreting-risk", package = "ambre")` walks through.