Skip to contents

Workflow: working with results and scenarios

Source: vignettes/articles/workflow.Rmd
workflow.Rmd

Once a model is built (see Model bricks and TOPIA I) and solved (see Solver backends and TOPIA II), a handful of utility functions do the rest of the day-to-day work: pull objects and results out, edit a piece and re-solve, organize scenarios on disk, and compare runs. This article tours that layer.

We use the packaged single-region TOPIA kit as a running example, and keep every store — scenarios, models, repositories and the registry — under a temporary folder. All four matter: the stores and the registry are written beside whatever path is redirected, so redirecting only the scenario path leaves models/ and a registry beside the article.

wf <- file.path(tempdir(), "wf")
set_scenarios_path(wf)                                  # where scenarios go
set_models_path(file.path(wf, "models"))                # ... and models
set_repositories_path(file.path(wf, "repositories"))    # ... and repositories
set_registry_file(file.path(wf, "energyRt_registry.csv"))

um  <- topia$modules$electricity$R1
mod <- newModel("TOPIA", data = um$repo,
                calendar = topia$modules$calendars$s4_h24,
                region   = um$regions,
                horizon  = topia$modules$horizons$base, discount = 0.05)
mod <- setGeoscale(mod, um$geoscale)
# interpolate + solve on the bundled GLPK; echo = FALSE keeps the log quiet
scen <- interpolate_model(mod, name = "BASE") |>
  solve_scenario(solver = solver_options$glpk, echo = FALSE)

Selecting objects: getObject()

getObject() returns the building-block objects held in a repository, model or scenario, filtered by class, name, description, region or any object slot. It is the object-level counterpart of getData() (which returns their data). By default it returns a named list keyed by object name; drop = TRUE unwraps a single match into the object itself.

repo <- um$repo
names(getObject(repo, class = "technology"))        # all technologies
#> [1] "ECOA" "EGAS" "ESOL" "EWIN" "EBIO"
names(getObject(repo, class = "supply", region = "W1"))
#> [1] "SUP_BIO" "RES_SOL" "RES_WIN"
getObject(repo, name = "ECOA", drop = TRUE)@invcost # the ECOA object itself
#>   vintage cluster region year invcost wacc payback eac retcost
#> 1    <NA>    <NA>   <NA>   NA    2000   NA      NA  NA      NA

Region matching reads @region slots and the region/src/dst columns of any data.frame slot, so it works for every class; region-agnostic objects (e.g. commodities) match every region. getObject() accepts a scenario too, so the same query works before or after solving.

Extracting data: getData()

getData() pulls parameter and variable data out of a scenario as tidy data.frames. With merge = TRUE it returns one long frame (a named list otherwise); every frame carries a scenario and a name column.

getData(scen, "vObjective", merge = TRUE)$value      # total system cost, MEUR
#> [1] 18351.9

gen <- getData(scen, "vTechOut", comm = "ELC", merge = TRUE)
head(gen[, c("scenario", "tech", "region", "year", "timeslice", "value")], 4)
#> # A tibble: 4 × 6
#>   scenario tech  region  year timeslice value
#>   <chr>    <chr> <chr>  <int> <chr>     <dbl>
#> 1 BASE     ECOA  W1      2025 WIN_h00   0.799
#> 2 BASE     ECOA  W1      2025 WIN_h01   0.736
#> 3 BASE     ECOA  W1      2025 WIN_h02   0.698
#> 4 BASE     ECOA  W1      2025 WIN_h03   0.673

The ... accept set filters, exact (comm = "ELC") or regex (comm_ = "^EL"). timeframe = c("lowest","highest","all") controls how timeslice-indexed values are returned. Rename dimensions or recode values on the way out with newNames = / newValues = (or afterwards with renameSets() / revalueSets()), and discover which parameters carry a given set with findData():

names(findData(scen, setsNames_ = "tech"))           # parameters indexed by 'tech'
#>  [1] "pTechWeatherAf"     "pTechWeatherAfs"    "pTechWeatherAfc"   
#>  [4] "pTechCap2act"       "pTechEac"           "pTechWacc"         
#>  [7] "pTechPayback"       "pTechEmisComm"      "pTechOlife"        
#> [10] "pTechFixom"         "pTechInvcost"       "pTechStock"        
#> [13] "pTechStockNew"      "pTechStockSurv"     "pTechVarom"        
#> [16] "pTechAf"            "pTechRampUp"        "pTechRampDown"     
#> [19] "pTechAfs"           "pTechGinp2use"      "pTechCinp2ginp"    
#> [22] "pTechUse2cact"      "pTechAct2AInp"      "pTechCap2AInp"     
#> [25] "pTechAct2AOut"      "pTechCap2AOut"      "pTechNCap2AInp"    
#> [28] "pTechNCap2AOut"     "pTechPho2AInp"      "pTechPho2AOut"     
#> [31] "pTechRet2AInp"      "pTechRet2AOut"      "pTechCinp2AInp"    
#> [34] "pTechCout2AInp"     "pTechCinp2AOut"     "pTechCout2AOut"    
#> [37] "pTechCact2cout"     "pTechCinp2use"      "pTechCvarom"       
#> [40] "pTechAvarom"        "pTechShare"         "pTechAfc"          
#> [43] "pTechCap"           "pTechNewCap"        "pTechRet"          
#> [46] "pTechRetCost"       "vTechAInp"          "vTechAOut"         
#> [49] "vTechAct"           "vTechCap"           "vTechEac"          
#> [52] "vTechFixom"         "vTechInp"           "vTechInv"          
#> [55] "vTechNewCap"        "vTechOut"           "vTechRetCost"      
#> [58] "vTechRetiredNewCap" "vTechRetiredStock"  "vTechPhaseOut"     
#> [61] "vTechStockCap"      "vTechStockPhaseOut" "vTechVarom"        
#> [64] "vTechEmsFuel"

find_in_model(mod, "ECOA") text-searches every object slot for a value — handy for locating where a name is used.

Editing an object: update()

update() edits the slots of a single model object (a technology, commodity, demand, …). It does not operate on a whole model or scenario — you update the object, put it back into the repository with add(..., overwrite = TRUE), and re-interpolate/solve.

ECOA <- getObject(repo, name = "ECOA", drop = TRUE)
ECOA <- update(ECOA, invcost = data.frame(invcost = 2500))  # pricier coal capex
repo_hi <- add(repo, ECOA, overwrite = TRUE)                # swap it back in

mod_hi  <- newModel("TOPIA_HI", data = repo_hi,
                    calendar = topia$modules$calendars$s4_h24,
                    region = um$regions, horizon = topia$modules$horizons$base,
                    discount = 0.05)

(For an already-interpolated scenario, update_parameter(scen, param, data) writes rows straight into its parameter store.)

Scenario folders and structure

(The Scenario management article explains this layer in depth — stores, manifests, the registry, runs, and variants; this section is the quick tour.)

The scenarios directory is an option — read it with get_scenarios_path(), set it with set_scenarios_path() (default "scenarios/"). Each scenario gets its own folder, named automatically as name-model-calendar-horizon (parts joined by dashes, which valid names cannot contain; duplicate parts are dropped, so a scenario and model sharing a name appear once):

get_scenarios_path()
#> [1] "C:\\Users\\admin\\AppData\\Local\\Temp\\RtmpKs1y8G/wf"
basename(scen@path)
#> [1] "BASE-TOPIA-s4_h24-base"

Saving a scenario writes an Arrow-backed folder that mirrors the scenario’s structure — a thin scen.RData shell, a scenario.yml manifest, each large data slot as a dataset, and one folder per solver run:

<scenarios-path>/<name-model-calendar-horizon>/
├── scenario.yml        # manifest: model ref/hash, calendar, horizon,
│                       #   `default:` — the default run (user-changeable)
├── scen.RData          # thinned S4 scenario shell
├── class · format · layout   # "scenario" · "parquet" · "3"
├── logfile.csv         # timestamped save log
├── model/data/<repo>/<object>/<slot>/   # model data (embedded scenarios only)
├── modInp/parameters/<param>/           # interpolated input parameters
└── runs/               # one folder per solve
    ├── glpk/           #   base-problem run, named by its solve label
    │   ├── run.yml     #     provenance: solver, cmdline, times, objective
    │   ├── input/      #     the exchange tables handed to the solver
    │   ├── output/     #     the solver's raw dump (read_solution() reads it)
    │   └── modOut/variables/<var>/      # the solved variables
    └── julia_highs/    #   another backend's run of the same problem

Beyond the per-scenario traces above, an optional operation log records the sequence of a working session: set_log_file("energyRt_log.csv") makes interpolate_model(), solve_scenario() and solve_recursive() append one line each (operation, object, status, objective, duration), and read_log() returns the sequence as a tibble. It is off by default.

Solving the same scenario with several backends (or option sets — pass run = "<label>" to solve_scenario()) keeps every run side by side; scenario_runs(scen) lists them with status and objective, and read_solution(scen, run = "<label>") switches the active solution. A saved scenario references its model in the model store when the model is a sealed version (seal_model()), and carries a copy of it otherwise (see ?save_scenario, mode, and the Scenario management article).

A scenario knows whether its data is in RAM or on disk via isInMemory(); when on disk, the slots are empty and read lazily from the folder on demand. Folders from older energyRt versions are still readable and migrate in place with upgrade_scenario_layout(path).

Saving and loading

save_scenario() spills the scenario to its folder and returns the (now on-disk) object; load_scenario(path, env = NULL) reads the shell back. Its data stays on disk until requested — getData() reads it lazily, or obj2mem() pulls the whole scenario back into memory.

saved <- save_scenario(scen, verbose = FALSE)
ld    <- load_scenario(saved@path, env = NULL, verbose = FALSE)

A saved scenario is also reachable by its registered name — load_scenario("BASE", env = NULL) resolves it through the project registry, and the cached accessor getScenario("BASE") is the everyday short form (see the Scenario management article):

identical(getScenario("BASE")@name, ld@name)
#> [1] TRUE
basename(saved@path)                                 # the scenario folder
#> [1] "BASE-TOPIA-s4_h24-base"
isInMemory(saved)                                    # FALSE -- data is on disk
#> [1] TRUE
getData(ld, "vObjective", merge = TRUE)$value        # lazy read, no full load
#> [1] 18351.9
ld <- obj2mem(ld)                                    # rehydrate fully into RAM
isInMemory(ld)
#> [1] TRUE

Every save_scenario() / save_model() / save_repository() also records its object in the project registry — a single CSV (default energyRt_registry.csv at the project root, see get_registry_file()) indexing saved models, repositories, scenarios, and runs. Repositories shared by several models can be stored once with save_repository() and referenced by each model (save_model(embed_repos = )), just as scenarios can reference a stored model. load_registry() reads it, find_in_registry() filters it (by type, name, or model hash), and refresh_registry() rebuilds it by rescanning the scenarios/ and models/ stores — the on-disk manifests are the source of truth, so the registry can always be reconstructed.

Comparing scenarios

Layer a policy lever onto the base model to get a second scenario, then pass a named list of scenarios to getData() — the scenario column makes the comparison a one-liner. (Scenarios already saved compare by name alone: getData(c("BASE", "CO2CAP"), ...) loads them from the registry on the fly.)

scen_cap <- interpolate_model(mod, "CO2CAP", um$CO2_CAP) |>   # add the CO2 cap lever
  solve_scenario(solver = solver_options$glpk, echo = FALSE)
emis <- getData(list(BASE = scen, CO2CAP = scen_cap), "vEmsFuelTot",
                comm = "CO2", merge = TRUE)

emis |>
  group_by(scenario, year) |> summarise(ktCO2 = sum(value), .groups = "drop") |>
  ggplot(aes(factor(year), ktCO2, fill = scenario)) +
  geom_col(position = "dodge") +
  labs(x = "year", title = "CO2 emissions by scenario") + theme_bw()

sapply(list(BASE = scen, CO2CAP = scen_cap),
       function(s) round(getData(s, "vObjective", merge = TRUE)$value[1]))
#>   BASE CO2CAP 
#>  18352  32669

A scenario that adds a lever to an interpolated one can be built on it: base = scen reuses BASE’s interpolated tables and computes only what the cap adds, and the result equals interpolating from scratch:

scen_cap2 <- interpolate_model(mod, "CO2CAP", um$CO2_CAP, base = scen)
sum(compare_inputs(scen_cap2, scen_cap)$summary$status != "identical")
#> [1] 0

The same holds for new technologies or commodities and for a sample of the base’s calendar or regions; anything else is interpolated in full, with a message saying why. Saved on a sealed base, such a scenario references the base’s tables instead of copying them (see Scenario management, “Scenarios built on a base”).

To reason about model size rather than results, model_size() estimates the parameter/variable/constraint counts of an interpolated scenario, and size() reports its in-memory footprint:

model_size(scen)
#> model_size: BASE
#>   parameters : 200 value, 359 maps, 15 sets
#>   param rows : 5,972
#>   estimate   : ~20,012 variables, ~20,704 constraints (from gating maps)
#>   top parameters by rows:
#>     pTimesliceAgg      1,188
#>     pTechCinp2use      1,152
#>     pWeather           1,004
#>     pTimesliceWeight   404
#>     pDemand            384
#>     pExportRowPrice    384
#>     pExportRow         384
#>     pStorageInpEff     384
#>     pStorageOutEff     384
#>     pTimesliceShare    101
#>     pTechFixom         20
#>     pTechInvcost       20
#>     pTechEac           17
#>     pTechStock         14
#>     pWacc              12
size(scen)
#> [1] "10.8 Mb"

For a systematic check that a build is correct and efficient, compare_interp_settings(mod) and compare_solve_settings(mod, solvers = ...) interpolate (and solve) the model under different storage settings and tabulate build size, time and — crucially — confirm the objective is invariant across fold/sparse/prune. These run several builds, so they are best run interactively:

compare_interp_settings(mod)                          # size/time by setting
compare_solve_settings(mod, solvers = list(solver_options$glpk))

Sampling: solving a subset of the model

Interpolation can solve a sample of a model without editing the model itself — in time (a sampled calendar) and in space (a sampled geoscale). Both work the same way: build the sample object, pass it to interpolate_model(), and the scenario is declared, validated and solved on the sample.

Sampled calendars

A sampled calendar is a row subset of the full calendar’s timetable plus an explicit year_fraction — the fraction of the year the sample covers. Timeslice weights annualise the sample inside the solver (the top slice carries weight 1/year_fraction), so the objective is an estimate of the full model’s — and exactly equals it when the sampled slices are representative:

tt  <- topia$modules$calendars$s4_h24@timetable
tt2 <- tt[tt$SEASON %in% c("WIN", "SUM"), ]      # two of the four seasons
cal2 <- newCalendar(timetable = tt2, name = "s2h24",
                    year_fraction = sum(tt2$share))

scen_t <- interpolate_model(mod, "S2", cal2) |>
  solve_scenario(solver = solver_options$glpk, echo = FALSE)
sapply(list(full = scen, sampled = scen_t),
       function(s) round(getData(s, "vObjective", merge = TRUE)$value[1]))
#>    full sampled 
#>   18352   20427

Two traps: pass year_fraction explicitly (it is not recomputed from the subset), and subset the existing timetable rather than rebuilding one with make_timetable() (which would renormalise the shares to a full year). Declarations are validated against the model’s calendar, so objects naming out-of-sample timeslices interpolate cleanly — their out-of-sample rows are simply filtered.

A sample of an already interpolated full-calendar scenario can be built on it, interpolate_model(mod, "S2", cal2, base = scen): the full scenario’s tables are taken filtered to the sampled timeslices instead of being recomputed.

Spatial sampling

The spatial mirror: pass a filtered geoscale — a geoscales::filter_geoscale() subset of the model’s regions — and the scenario becomes a sub-territory model: only the sampled regions are declared, out-of-sample parameter rows are filtered, and trade routes crossing the sample boundary are dropped (with a message). Unlike calendar sampling there is no reweighting: regional quantities are extensive, so the objective is the sub-territory’s own, and disjoint samples add up to the full model (absent cross-boundary trade).

TOPIA is single-region, so a compact three-region demo:

library(geoscales)

r3 <- c("R1", "R2", "R3")
gs <- geoscale_from_leaftable(
  data.frame(zone = c("Z12", "Z12", "Z3"), region = r3, km2 = 1),
  geoframes = c("zone", "region"), name = "demo")

m3 <- newModel(
  "M3", region = r3, horizon = newHorizon(2020), discount = 0.05,
  repo = newRepository(
    "r",
    newCommodity("COA", timeframe = "ANNUAL"),
    newCommodity("ELC", timeframe = "ANNUAL"),
    newSupply("SUP_COA", commodity = "COA",
              supply = data.frame(region = r3, cost = 5)),
    newTechnology("ECOA", input = list(comm = "COA"),
                  output = list(comm = "ELC"),
                  af = data.frame(af.up = 0.9),
                  invcost = data.frame(invcost = 1000),
                  olife = list(olife = 30), cap2act = 1),
    newDemand("DEM_ELC", commodity = "ELC",
              demand = data.frame(region = r3, demand = c(10, 20, 30))),
    newTrade("TR23", commodity = "ELC", cap2act = 1,
             routes = data.frame(src = "R2", dst = "R3"),
             trade  = data.frame(src = "R2", dst = "R3", ava.up = 4),
             invcost = data.frame(invcost = 10))))

Filtering the geoscale to R1 + R2 and interpolating drops the R2 -> R3 route on the way:

gs12 <- filter_geoscale(gs, "region", c("R1", "R2"),
                        drop_empty_geoframes = TRUE)
scen12 <- interpolate_model(m3, "S12", gs12)
#> spatial sample: sub-territory solve on 2 of 3 regions
#> spatial sample: kept 2 region(s); dropped 1 object(s) and 1 trade route(s) outside the sample (verbose = TRUE lists them)
as.character(scen12@settings@region)
#> [1] "R1" "R2"

A dropped route can optionally leave a priced boundary stub at the kept endpoint — an export (kept source) or import (kept destination) object bounded by cap.up (or, failing that, the route’s ava.up) — via boundary_prices; routes without a price row are just dropped:

bp <- data.frame(src = "R2", dst = "R3", price = 50, cap.up = 4)
scen12p <- interpolate_model(m3, "S12P", gs12, boundary_prices = bp)
#> spatial sample: sub-territory solve on 2 of 3 regions
#> spatial sample: kept 2 region(s); dropped 1 object(s) and 1 trade route(s) outside the sample (verbose = TRUE lists them)
grep("^EXP", as.character(scen12p@modInp@sets$expp), value = TRUE)
#> [1] "EXP_TR23_R22R3"

The same surgery is available up front, as a model-to-model verb:

m12 <- subset_model_regions(m3, c("R1", "R2"), boundary_prices = bp)
as.character(m12@config@region)
#> [1] "R1" "R2"

The filtered geoscale records its coverage (fraction of each weight’s parent total — see geoscales::geoscale_coverage()), and its name is mangled ("demo[region:R1+R2]") so different samples of one parent never collide in caches or joins. A pruned geoscale (a coarser level of the model’s own hierarchy) requests a full-territory solve at the parent level instead. That is a different operation – it aggregates the model’s data rather than subsetting it – so interpolate_model() refuses a pruned geoscale and points at aggregate_model_regions(), which builds the coarser model directly.

A trade window instead of a flat price

A stub priced with a single number is a perfectly elastic counterparty: the region can buy or sell any amount the bound allows without moving the price. For supply-curve work that is usually too generous. boundary_window() builds a stepped curve instead — a rising cost of imports, a falling revenue on exports — and sizes its quantity from the region’s own demand:

bw <- boundary_window(m3, regions = c("R1", "R2"),
                      share = 0.1, nsteps = 3, price = 40)
head(bw, 2)
#>   src dst year timeslice price cap.up nsteps price_lo price_hi
#> 1  R2  R3   NA      <NA>    40      2      3       20       40

share = 0.1 caps the boundary at 10% of the region’s demand in each timeslice, which is where the import/export bound lives; summed over the year that is 10% of annual demand. The window deliberately refuses to guess a price — it can size the quantity from the model, but what the neighbours would charge is a judgement, so price is required.

The result is an ordinary boundary_prices table with three extra columns (nsteps, price_lo, price_hi), so it goes where a flat one goes; omit them and you get exactly the flat stub as before.

m12w <- subset_model_regions(m3, c("R1", "R2"), boundary_prices = bw)
stub <- m12w@data[["boundary_stubs"]]@data[[1]]
# `@cluster` declares WHAT the steps are; the split itself is enforced in the
# quantity bound, so read it there rather than from a descriptive copy
as.data.frame(stub@cluster)[, c("cluster", "desc", "order")]
#>   cluster                                         desc order
#> 1      S1 step 1 of 3: 0-0.3333 of quantity at 36.6667     1
#> 2      S2 step 2 of 3: 0.3333-0.6667 of quantity at 30     2
#> 3      S3 step 3 of 3: 0.6667-1 of quantity at 23.3333     3

Region by region

solve_by_region() runs the whole sweep — one region at a time by default, or one group at a time — storing each as a variant of a single scenario:

br <- solve_by_region(m3, name = "sweep")           # autarky
br <- solve_by_region(m3, name = "sweep_w",         # with a trade window
                      trade = "auto", price = 40, share = 0.1)
br$runs
getData(br, "vTechCap", merge = TRUE)               # stacked, region column

With trade = "none" the regions are disjoint and nothing crosses, so the per-region objectives add up to the full model’s — the same property the sample objectives have above. Any trade window breaks it: each region then faces an exterior the others cannot see. br$additive records which case you are in, and print() says so rather than leaving it to be discovered.

on_error = "continue" keeps the sweep going past a region that fails to solve, which matters at 30-odd regions: a single infeasible one should not cost you the rest.

Sampling the calendar repeatedly

A single sampled calendar gives one cheap estimate. solve_by_sample() runs a whole ensemble of them — consecutive blocks tiling the year, disjoint random draws, or bootstrap draws — and keeps every run:

ens <- solve_by_sample(m3, name = "ens", sample_size = 2,
                       method = "sequential")     # two-season blocks
ens$runs
sample_summary(ens, "vTechCap", probs = c(0.05, 0.5, 0.95))

Samples are drawn as whole members of a calendar level — seasons, weeks, days. A hand-edited sample need not be a complete cross: its calendar is irregular and keeps the full calendar’s timeslice names. calendar_samples() builds the specification on its own, so it can be inspected or edited first.

Two rules follow from how a sampled calendar works, and both are enforced rather than left to be discovered. A sample must keep at least two members at every level: a single season (or week, or day) collapses its level, and the only way to build such a calendar renames the leaf timeslices so the model’s own data would stop matching them. And every sample is annualised — pTimesliceWeight = 1/year_fraction — so each one estimates the whole year. The ensemble is therefore a set of replicates whose spread is the sampling uncertainty; summing them is meaningless, unlike the region sweep above, whose disjoint pieces do add up.

Combining a sequence’s results

Each driver stores one variant per unit of work and nothing merged is written to disk. Combination happens in memory, on read, and the rule differs by driver because the units mean different things.

getData() on a sequence stacks and does not aggregate. It binds every unit’s rows and adds a provenance column — region for solve_by_region(), sample for solve_by_sample() — so the caller decides what to do with them. Two methods are exceptions: getData() on a recursive sequence stitches, taking each year from the step that decided it, and on a guided sequence it delegates to the final full-horizon solve.

rc <- solve_recursive(m3, name = "rec")  # milestone by milestone

getData(br, "vTechCap", merge = TRUE)    # stacked, with a `region` column
getData(ens, "vTechCap", merge = TRUE)   # stacked, with a `sample` column
sample_summary(ens, "vTechCap", by = "tech")
recursive_objective(rc)
driver what a unit is how results combine
solve_by_region() a disjoint piece of the territory quantities add; objectives sum only under trade = "none", which the result reports as additive
solve_by_sample() a replicate estimate of the whole year a distribution — mean and quantiles via sample_summary(). Summing is wrong, because every sample is already annualised
solve_recursive() a milestone window stitched, each year contributed once. Per-step objectives are not comparable, so recursive_objective() recomputes a discounted total
solve_guided() a stage of a ladder the answer is $final; guided_gap() prices the guidance against perfect foresight

Two things the combination does not do. Regional objectives sum to the whole system only in autarky: with transmission the per-region sum exceeds the full model, and the difference is the network’s value, so additive = TRUE means “may be summed”, not “equals the full model”. And under trade = "auto" each region faces an exterior the others cannot see — one region’s export and another’s import are independent decisions at independently set prices, and nothing reconciles them.

Seeing what a sequence covered

A sequence records its design per unit, and in a different shape for each driver. sequence_coverage() puts all four on one frame — one row per unit and element covered, with dim naming the scale the units were drawn on — and autoplot() draws that frame on the scale itself.

sequence_coverage(ens)             # unit, dim, key, role, count
autoplot(ens)                      # the calendar, tiled by sample
autoplot(ens, "wall", year = 2020) # the same, as a wall calendar
autoplot(br)                      # the geoscale, filled by unit
autoplot(rc)                      # the horizon, decided against look-ahead

Three properties show in the drawing that a table of runs leaves implicit. Elements no unit covered keep a grey fill, so a gap in an ensemble is a hole in the picture rather than an absent row. An element two units claimed is filled as (shared) instead of being given to the first of them, and the subtitle counts both that and the within-unit multiplicity a bootstrap draw carries. And a recursive step’s look-ahead milestones are filled apart from the years it decided, which is the overlap between steps seen for what it is.

The wall calendar needs a day-resolution timeframe (YDAY or MDAY) and rolls each unit up to whole days; the tile grid, over the calendar’s two finest timeframes, works for any calendar. A region sample is mapped on the finest geoframe that contains every solved region — a model’s regions need not sit at the finest frame of its geoscale — and falls back to a plain list of regions where the geoscale carries no geometry.

The spread of a result across an ensemble is a histogram rather than a map:

autoplot(ens, "histogram")                                  # the objective
autoplot(ens, "histogram", var = "vTechCap", by = "tech")

Bars are weighted by each sample’s year_fraction, since sequential sampling can leave one block covering more of the year than its siblings, and the dashed rule is the weighted mean sample_summary() reports — one weighting rule, so the plot and the table agree about the same ensemble. Where ncap_floor > 0 the free pass is overlaid on the floored one. A variable must reduce to one value per sample within each group; a dimension left inside a bar is named rather than summed away, since totalling capacity over years would answer a different question. There is no histogram for solve_by_region(), whose units are shards of one system: their values are addends, and their spread would report only that regions differ in size.

Watching a sequence run

A sequence solve is the longest-running operation in the package and reports as one. set_progress_bar() registers a progressr handler, after which each driver advances a bar one unit at a time; verbose = TRUE prints a line per unit, as an alternative or in addition. Only the outermost sequence claims the bar, so the interpolation inside a driver’s loop does not open a second one over the first.

set_progress_bar()
ens2 <- solve_by_sample(m3, name = "ens_bar", sample_size = 2,
                        verbose = FALSE)

See also