## ----include = FALSE----------------------------------------------------------
# Nothing here is evaluated when the vignette is built: every step reads data
# the package does not ship, and the tuning step trains for hours. The figures
# are files in figures/, drawn once from real data (see the package's GitHub
# repository, tools/vignette_figures.R).
knitr::opts_chunk$set(eval = FALSE, collapse = TRUE, comment = "#>")

## ----include = FALSE, eval = TRUE---------------------------------------------
# The application's figures are drawn from its run (tools/vignette_results.R)
# and committed with its numbers; until then a note stands in for each, so
# the vignette builds without them.
result_figure <- function(file, figure) {
  if (file.exists(file)) return(knitr::include_graphics(file))
  knitr::asis_output(paste0("*", figure, " is drawn from the application run once it completes.*\n"))
}

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 1. The workflow keeps data preparation, validation, tuning and final prediction explicit."----
knitr::include_graphics("figures/workflow.png")

## -----------------------------------------------------------------------------
# library(soilcnn)

## -----------------------------------------------------------------------------
# store <- dsm_prepare(
#   points     = "profiles.gpkg",        # an sf object, a data frame, or a file sf reads
#   target     = "soc_stock",
#   raster_dir = "predictors/",          # aligned single-band rasters, one per predictor
#   windows    = c(3, 9, 15),            # in pixels: 0.75, 2.25 and 3.75 km at 250 m
#   out_dir    = "outputs",
#   transform  = "log1p",                # the space the network trains in
#   target_min = 0
# )

## -----------------------------------------------------------------------------
# data <- dsm_load(store)

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 2. The full application dataset: profile locations and the predictor footprint. Uneven coverage motivates spatial validation and applicability diagnostics. The map uses a spherical Lambert azimuthal equal-area projection centred at 75 W / 15 S, with equal horizontal and vertical scales."----
knitr::include_graphics("figures/observations.png")

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 3. Landscape context and native-cell windows from the high-resolution example. The marker identifies the central cell. Bars show source raster values; elevation and clay labels are rounded for readability. Window outlines are enlarged in the detail, using solid (3), dashed (9) and dotted (15) borders. This illustration is separate from the 250 m trial."----
knitr::include_graphics("figures/patches.png")

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 4. Conceptual multiscale architecture. Predictor scaling is fitted on the training rows of each fold."----
knitr::include_graphics("figures/architecture.png")

## -----------------------------------------------------------------------------
# cv <- spatial_cv(k = 5, block_size = "auto", buffer = "auto", test_frac = 0.15,
#                  calibration_frac = 0.15)

## -----------------------------------------------------------------------------
# random_cv(k = 10)                         # ignores geography, on purpose
# holdout_cv(validation_frac = 0.2)         # a single split
# region_cv(group = data$points$biome)      # leave one region out
# knndm_cv(k = 5, predpoints = prediction_sample(terra::rast("predictors/elev.tif")))

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 5. Validation geometry: fold 1 of each design, with one test set and one calibration set for all five. Hollow points are not used in this fold. All panels use the same equal-area projection and extent."----
knitr::include_graphics("figures/designs.png")

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 6. The buffer in the full spatial plan. Grey training points stay outside the shaded exclusion zone; hollow points were removed from training by a buffer."----
knitr::include_graphics("figures/buffer.png")

## -----------------------------------------------------------------------------
# plan <- resolve_resampling(cv, data)
# print(plan)

## -----------------------------------------------------------------------------
# fit <- dsm_train(
#   data,
#   model       = "cnn",
#   resampling  = plan,
#   tune_length = 30,      # a budget, as in caret, not a lattice
#   n_seeds     = 3,       # what turns "A beat B" into a claim with an error bar
#   output_dir  = "outputs/tuning"
# )
# 
# fit$by_config                                   # mean and sd, one row per config
# print_noise_floor(seed_noise_floor(fit$comparison))
# print_one_se(one_se(fit$by_config))             # the simplest config within one SE

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 7. A didactic, simulated example of the one-standard-error rule. F has the highest mean; D is the smallest configuration above the 0.490 threshold. Bars show plus/minus one SE. These values are not results of the SOC application."----
knitr::include_graphics("figures/selection.png")

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 8. The configurations the spatial design tried in the application, with the one-SE threshold and the configuration it selects."----
result_figure("figures/trial_selection.png", "Figure 8")

## -----------------------------------------------------------------------------
# rf_context <- dsm_train(data, model = "rf", resampling = plan,
#                         features = c("centre", "window_mean"), tune_length = 4,
#                         output_dir = "outputs/tuning")
# rf_point   <- dsm_train(data, model = "rf", resampling = plan,
#                         features = "centre", tune_length = 4, output_dir = "outputs/tuning")

## -----------------------------------------------------------------------------
# caret_available("^xgb")
# register_model(caret_spec("xgbTree"))
# xgb <- dsm_train(data, model = "xgbTree", resampling = plan, tune_length = 4,
#                  output_dir = "outputs/tuning")

## -----------------------------------------------------------------------------
# freeze_selection(fit$run_dir, config_id = one_se(fit$by_config)$config_id[1])
# score_test_grid(fit$run_dir, data, device = setup_torch_device())

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 9. Each design's final model on the same test profiles. Main panels share a zoom covering the central distribution; the insets show every profile on log scales, their dotted square marking that zoom. Metrics use all profiles. The MAE differences relative to spatial blocks have 95% intervals from whole 100 km blocks."----
result_figure("figures/trial_designs.png", "Figure 9")

## -----------------------------------------------------------------------------
# final <- dsm_final(fit, seeds = 10)

## -----------------------------------------------------------------------------
# final$per_config[[1]]$intervals$summary

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 10. The 90% intervals of the five designs on the same test profiles: coverage against the nominal 90% level, and the mean width, for each calibration at a constant width and at one that follows the level and the DI."----
result_figure("figures/trial_intervals.png", "Figure 10")

## -----------------------------------------------------------------------------
# map <- dsm_predict(final, data, extent = c(-56, -50.5, -15, -13.8),
#                    intervals = c("cv", "split"))
# map

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 11. The application's map: the ensemble median, the width of the 90% conformal interval, and the area of applicability. Continuous colour scales use the 2nd and 98th percentiles of the same regular sample used in the text, computed before reprojection; values beyond the limits are clipped for display."----
result_figure("figures/trial_maps.png", "Figure 11")

## -----------------------------------------------------------------------------
# themes <- read.csv2("importance_themes.csv")[, c("channel", "variable")]
# perm <- dsm_importance(final, data, groups = themes)
# shap <- dsm_importance(final, data, shap_importance(), groups = themes)
# compare_importance(permutation = perm, shap = shap)
# plot(shap)

## ----eval = TRUE, echo = FALSE, out.width = "100%", fig.cap = "Figure 12. What the application's model learned, by theme: the drop in CCC when a theme is permuted and its mean absolute SHAP value in separate panels with their original units, and signed SHAP maps on one shared scale for the leading themes."----
result_figure("figures/trial_importance.png", "Figure 12")