---
title: "Getting started with apsimeval"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with apsimeval}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)
```

## The interface

Everything below can also be done from a browser interface:

```{r}
library(apsimeval)
run_app()                                   # opens on the working directory
run_app(out_dir = "~/apsim/sundarbans/outputs")
run_app(out_dir = "outputs", project = "sundarbans.apsimproj")
run_app(out_dir = "outputs", poll_ms = 1000, launch.browser = FALSE)
```

`run_app()` needs the optional packages: `install.packages(c("shiny", "DT", "bslib"))`.
Leave it open while APSIM runs — with **Live sync** on it re-reads only the
files that changed and redraws, keeping your current settings.

## Reading output

```{r}
sim <- read_out_dir("outputs")              # every .out in the folder
apsim_vars("outputs/N120.out")              # variables and units, header only
sim <- read_out_dir("outputs", vars = c("wagt", "wrr"))
```

## Treatments from simulation names

```{r}
detect_factors(sim$sim)                     # guesses the separator and fields
preview_factors(sim$sim, parts = c("site", "crop", "water", "N"))
sim <- set_factors(sim, parts = c("site", "crop", "water", "N"),
                   numeric_from = "N")
```

Irregular names take a regex, hopeless ones a lookup table:

```{r}
sim <- set_factors(sim, regex = "(?<site>[A-Za-z]+)[-_]N(?<N>\\d+)")
map <- factor_template(sim$sim, factors = c("water", "N"))   # fill in and pass back
sim <- set_factors(sim, map = map)
```

## Calibration and validation

```{r}
obs <- map_obs(read_obs_raw("observed.xlsx"),
               date_col = "Sampling_Date", value_cols = c("Grain_yield", "Biomass"),
               treatment_col = "Treatment",
               # Spreadsheet headings are not APSIM variable names. Without
               # this, nothing pairs and nothing plots.
               var_map = c(Grain_yield = "wrr", Biomass = "wagt"))

check_variables(obs, sim)                   # catches variable name mismatches
check_levels(obs, sim)                      # catches treatment label mismatches

long <- as_long(sim, vars = "wrr")
pr   <- pair_obs(long, obs, tol_days = 3)

cut <- suggest_cutoff(pr)                   # first season calibrates
pr  <- split_phase(pr, "date", at = cut)

# Statistics are never pooled across variables: RMSE would mix units and R2
# would be inflated by the gap between the variables rather than model skill.
gof("obs", "pred", data = pr, by = "variable")
gof("obs", "pred", data = pr, by = c("variable", "phase"))
gof("obs", "pred", data = pr, by = c("variable", "phase", "N"))
```

The series stays continuous across the cutoff:

```{r}
plot_ts(long, obs = obs, cutoff = cut)
# Colouring is a presentation choice and does not affect the statistics; the
# corner box follows `stats_by`, which defaults to the faceting column.
plot_one2one(pr, group = "N", stats_by = character(0))   # pooled box
plot_one2one(pr, group = "N", facet = "phase", sd_band = TRUE)
```

## Scenarios and sensitivity

```{r}
ag <- aggregate_sim(sim, by = c("sim", "year", "water", "N"),
                    vars = c(wrr = "max"), crop_only = TRUE)
plot_exceedance(ag, "wrr", group = "N", facet = "water", ref = 7000)
plot_bar(ag, "wrr", "N", letters = TRUE)

sens_anova(ag, "wrr", c("site", "water", "N"), order = 2)
sens_oat(ag, "wrr", c("site", "water", "N"))
```

## Figures

```{r}
set_labels(list(wrr = "Grain yield"))
set_palette(c(N0 = "#D55E00", N120 = "#0072B2"))
set_colour_mode("grey")                     # "colour", "grey", or "mono"

fig <- build_figure(list(
  panel_spec("one2one", group = "N", facet = "phase"),
  panel_spec("exceedance", var = "wrr", group = "N", ref = 7000),
  panel_spec("bar", var = "wrr", group = "N", letters = TRUE)
), sim = sim, paired = pr, ncol = 3)

save_pub(fig, "figure_2", width = "double", height = 80,
         dpi = 600, format = c("tiff", "pdf"))
```

## Projects

```{r}
p <- new_project("sundarbans_rice", "outputs")
p <- set_project_factors(p, sim)
p$labels <- list(wrr = "Grain yield"); p$colour_mode <- "grey"
save_project(p, "sundarbans_rice.apsimproj")

p <- load_project("sundarbans_rice.apsimproj")
sim <- apply_factors(read_out_dir("outputs"), p)   # after a re-run
```

## Replicated observed data

Field data usually arrives with one row per plot:

| Sampling_Date | Site | Treatment | Replication | Biomass |
|---|---|---|---|---|
| 2016-08-01 | Gosaba | N120 | 1 | 1000 |
| 2016-08-01 | Gosaba | N120 | 2 | 1100 |

Name the replication column and the replicates are averaged per date and
treatment, keeping the spread for error bars:

```{r}
raw <- read_obs_raw("observed.xlsx")
attr(raw, "guess")                # suggested date / treatment / rep / value columns

obs <- map_obs(raw, date_col = "Sampling_Date", value_cols = "Biomass",
               treatment_col = c("Site", "Treatment"),   # joined to match sim names
               rep_col = "Replication")
```

`obs` gains `obs_sd` and `obs_n`. A date with one replicate gets `obs_sd = NA`
rather than zero, so no misleading zero-length bar is drawn.

```{r}
plot_ts(long, obs = obs, cutoff = cut)          # error bars drawn automatically
plot_ts(long, obs = obs, cutoff = cut, obs_sd = NULL)   # suppress them
plot_one2one(pr, obs_sd = "obs_sd")             # horizontal bars on the scatter
```

If your sheet already holds means and a separate SD column (`Biomass_SD`),
leave `rep_col` empty: the SD column is detected and paired automatically.

## Timeline graphs

`plot_series()` handles several variables and several treatments at once.

Different units go on a secondary axis:

```{r}
long <- as_long(sim, vars = c("wagt", "pond"), keep = c("water", "N"))
plot_series(long, vars = c("wagt", "pond"), secondary = "pond")
```

Or give each variable its own row, which avoids the dual-axis problem entirely:

```{r}
plot_series(long, vars = c("wagt", "wrr", "pond"), layout = "stacked")
```

All treatments side by side — the usual way to judge a treatment effect,
since the shapes are compared rather than just the endpoints:

```{r}
plot_series(long, vars = c("wagt", "pond"), secondary = "pond",
            facet_by = "sim", ncol = 2)             # every simulation
plot_series(long, vars = "wagt", facet_by = "water") # or a factor column
plot_series(long, vars = "wagt", facet_by = "water",
            colour_by = "treatment")                # colour instead of panels
```

Combine layouts for a variable-by-treatment grid, and add the cutoff, observed
points and replicate error bars as usual:

```{r}
plot_series(long, vars = c("wagt", "wrr"), layout = "stacked", facet_by = "sim")
plot_series(long, vars = "wagt", facet_by = "sim", obs = obs, cutoff = cut)
```

Note `keep =` in `as_long()`: treatment columns have to survive the pivot or
there is nothing left to facet by.

A caution on the secondary axis. Two series with different units can be made
to look correlated by the scaling alone, so `layout = "stacked"` is the safer
choice when the variables are not meant to be compared point by point. The
right axis is always labelled in original units, and legend entries drawn
against it are marked.

### Observed data on a subset of the variables

Observations rarely cover every plotted variable — biomass and yield are
measured, AWD ponding usually is not. Pass whatever you have; variables
without observations are drawn as simulated-only lines and the omission is
reported:

```{r}
plot_series(long, vars = c("wagt", "wrr", "pond"), secondary = "pond", obs = obs)
#> No observed data for: pond (simulated only).
```

Each variable's points carry their own marker shape and their own replicate
error bars, computed exactly as elsewhere in the package (`rep_col` in
`map_obs()`, or an existing `obs_sd` column). Points belonging to a secondary
axis variable are rescaled with the line, and their error bars are rescaled as
lengths rather than positions.

When panelling by treatment, observed records need to know which panel they
belong to. If the observed data has no column matching `facet_by`, it is
derived from the simulation label:

```{r}
plot_series(long, vars = "wagt", facet_by = "water", obs = obs)  # obs has only `sim`
```

Records matching no simulation are dropped and counted in a message rather
than being silently discarded.