## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)
set.seed(20260507)

## ----load---------------------------------------------------------------------
library(leafwax)

example_record_path <- function(filename) {
  installed_path <- system.file(
    "extdata", "example_records", filename,
    package = "leafwax"
  )
  if (nzchar(installed_path)) {
    return(installed_path)
  }

  source_paths <- file.path(
    c("inst", file.path("..", "inst")),
    "extdata", "example_records", filename
  )
  source_path <- source_paths[file.exists(source_paths)][1]
  if (!is.na(source_path)) {
    return(source_path)
  }

  stop("Could not locate example record: ", filename, call. = FALSE)
}

sugan_path <- example_record_path("LS13WASU_C29_d2H.csv")
sugan <- read.csv(sugan_path)
sugan$d2h_wax <- sugan$d2H_wax
sugan$age     <- sugan$age_yrBP

sonk_path <- example_record_path("LS14LASO_C29_d2H.csv")
sonk <- read.csv(sonk_path)
sonk$d2h_wax <- sonk$d2H_wax
sonk$age     <- sonk$age_yrBP

c(sugan_n        = nrow(sugan),
  sugan_range_pm = round(diff(range(sugan$d2h_wax)), 1),
  sonk_n         = nrow(sonk),
  sonk_range_pm  = round(diff(range(sonk$d2h_wax)), 1))

## ----plot-wax, fig.height = 6, echo = TRUE------------------------------------
op <- par(mfrow = c(2, 1), mar = c(4.5, 5.4, 2.2, 1), mgp = c(3.4, 0.8, 0))

wax_ylim <- extendrange(c(sugan$d2h_wax, sonk$d2h_wax), f = 0.05)

plot_wax <- function(record, title, boundary, ylim) {
  ord <- order(record$age)
  plot(
    record$age[ord],
    record$d2h_wax[ord],
    type = "o", pch = 16, col = "black",
    xlim = rev(range(record$age)),
    ylim = ylim,
    xlab = "Age (yr BP)",
    ylab = expression(delta^2 * H[wax] ~ "(‰)"),
    main = title
  )
  abline(v = boundary, lty = 2, col = "red")
}

plot_wax(sugan, "Lake Sugan (LS13WASU, C29 n-alkane)",
         boundary = 800, ylim = wax_ylim)
plot_wax(sonk, "Sonk11D (LS14LASO, C29 n-alkane)",
         boundary = 5000, ylim = wax_ylim)

par(op)

## ----slope--------------------------------------------------------------------
sugan_lon <- 93.95;   sugan_lat <- 38.8667
sonk_lon  <- 75.1961; sonk_lat  <- 41.7939

slope_sugan <- suppressWarnings(local_effective_slope(
  longitude = sugan_lon, latitude = sugan_lat,
  model_name = "baseline_sp", n_draws = 100,
  verbose = FALSE
))

slope_sonk <- suppressWarnings(local_effective_slope(
  longitude = sonk_lon, latitude = sonk_lat,
  model_name = "baseline_sp", n_draws = 100,
  verbose = FALSE
))

rbind(
  sugan = quantile(slope_sugan, c(0.025, 0.5, 0.975)),
  sonk  = quantile(slope_sonk,  c(0.025, 0.5, 0.975))
)

## ----invert-------------------------------------------------------------------
recon_sugan <- suppressWarnings(invert_d2H(
  d2H_wax    = sugan$d2h_wax,
  d2H_wax_sd = rep(3, nrow(sugan)),
  longitude  = rep(sugan_lon, nrow(sugan)),
  latitude   = rep(sugan_lat, nrow(sugan)),
  model_name = "baseline_sp",
  n_posterior_draws = 100,
  slope        = slope_sugan,
  record_id    = "LS13WASU",
  return_full  = TRUE,
  verbose      = FALSE
))

recon_sonk <- suppressWarnings(invert_d2H(
  d2H_wax    = sonk$d2h_wax,
  d2H_wax_sd = rep(3, nrow(sonk)),
  longitude  = rep(sonk_lon, nrow(sonk)),
  latitude   = rep(sonk_lat, nrow(sonk)),
  model_name = "baseline_sp",
  n_posterior_draws = 100,
  slope        = slope_sonk,
  record_id    = "LS14LASO",
  return_full  = TRUE,
  verbose      = FALSE
))

## ----detect-------------------------------------------------------------------
rho_sugan <- estimate_temporal_autocorrelation(
  sugan$d2h_wax, sugan$age, method = "ar1"
)

dc_sugan <- detect_change(
  reconstruction    = recon_sugan,
  age               = sugan$age,
  baseline_interval = c(-100, 800),
  test_intervals    = list(older = c(800, 1700)),
  sigma_residual    = 16,
  sigma_analytical  = 3,
  rho_t             = rho_sugan,
  beta_eff          = stats::median(slope_sugan),
  confidence        = 0.95,
  magnitudes        = c(10, 30, 50)
)

rho_sonk <- estimate_temporal_autocorrelation(
  sonk$d2h_wax, sonk$age, method = "ar1"
)

dc_sonk <- detect_change(
  reconstruction    = recon_sonk,
  age               = sonk$age,
  baseline_interval = c(4000, 5000),
  test_intervals    = list(early_holocene = c(5000, 6000)),
  sigma_residual    = 16,
  sigma_analytical  = 3,
  rho_t             = rho_sonk,
  beta_eff          = stats::median(slope_sonk),
  confidence        = 0.95,
  magnitudes        = c(10, 30, 50)
)

list(
  sugan = list(rho_t = round(rho_sugan, 3),
               threshold_permil = round(dc_sugan$threshold, 1),
               intervals        = dc_sugan$intervals),
  sonk = list(rho_t = round(rho_sonk, 3),
              threshold_permil = round(dc_sonk$threshold, 1),
              intervals        = dc_sonk$intervals)
)

## ----assess-------------------------------------------------------------------
build_claim <- function(beta_eff, rho_t, baseline, test, magnitude_precip) {
  list(
    level             = 4,
    interval_baseline = baseline,
    interval_test     = test,
    sigma_analytical  = 3,
    rho_t             = rho_t,
    confidence        = 0.95,
    beta_eff          = beta_eff,
    magnitude_precip  = magnitude_precip,
    # Level 2 integrity gates: both must be ruled out by independent
    # record-specific evidence regardless of which Level 2 path is used.
    sediment_source_ruled_out = list(
      value    = TRUE,
      evidence = "grain size and mineralogy stable across the interval"
    ),
    depositional_artifact_ruled_out = list(
      value    = TRUE,
      evidence = "continuous varved sequence with no erosional unconformity"
    ),
    # Level 2 path (a): corroborating evidence against vegetation
    # reorganization.
    corroborating_proxies = list(
      regional_proxy = "regional records show coeval shift"
    ),
    vegetation_stationary = list(
      value    = TRUE,
      evidence = "n-alkane chain length distributions stable across the boundary"
    ),
    seasonal_source_stationary = list(
      value    = TRUE,
      evidence = "regional d18O record shows no seasonality shift"
    ),
    evapotranspirative_stationary = list(
      value    = TRUE,
      evidence = "leaf-water proxy stable; no aridity transition"
    )
  )
}

sugan_record <- data.frame(
  d2h_wax     = sugan$d2h_wax,
  age         = sugan$age,
  d2h_wax_err = rep(3, nrow(sugan))
)
sonk_record <- data.frame(
  d2h_wax     = sonk$d2h_wax,
  age         = sonk$age,
  d2h_wax_err = rep(3, nrow(sonk))
)

verdict_sugan <- suppressWarnings(assess_claim(
  record         = sugan_record,
  claim          = build_claim(stats::median(slope_sugan),
                                rho_sugan,
                                c(-100, 800), c(800, 1700),
                                magnitude_precip = 30),
  reconstruction = recon_sugan
))

verdict_sonk <- suppressWarnings(assess_claim(
  record         = sonk_record,
  claim          = build_claim(stats::median(slope_sonk),
                                rho_sonk,
                                c(4000, 5000), c(5000, 6000),
                                magnitude_precip = 30),
  reconstruction = recon_sonk
))

c(sugan_highest_level  = verdict_sugan$highest_level,
  sugan_supported_at_4 = verdict_sugan$asserted_supported,
  sonk_highest_level   = verdict_sonk$highest_level,
  sonk_supported_at_4  = verdict_sonk$asserted_supported)

## ----envelope-----------------------------------------------------------------
# Hypothetical OIPC value at Sonk11D. In a real analysis, extract from
# the OIPC raster at (sonk_lon, sonk_lat) using terra::extract().
sonk_oipc_ref <- -75   # per mil, illustrative only

# Hypothetical regional pollen scenario: a 30 percentage-point
# woody-to-grass transition across the interval, with a moderate C4
# increase. Names must be exactly tree, shrub, grass, C4 (case-sensitive).
env_sonk <- suppressWarnings(compute_vegetation_envelope(
  oipc_ref   = sonk_oipc_ref,
  from       = c(tree = 0.4, shrub = 0.3, grass = 0.2, C4 = 0.05),
  to         = c(tree = 0.1, shrub = 0.2, grass = 0.5, C4 = 0.20),
  model_name = "full_interact_sp",
  n_draws    = 100,
  verbose    = FALSE
))

c(envelope_median   = env_sonk$envelope_median,
  envelope_p975_abs = env_sonk$envelope_p975_abs)

## ----magnitude-path-claim-----------------------------------------------------
sonk_claim_magnitude <- list(
  level             = 2,
  interval_baseline = c(4000, 5000),
  interval_test     = c(5000, 6000),
  sigma_analytical  = 3,
  rho_t             = rho_sonk,
  confidence        = 0.95,
  sediment_source_ruled_out = list(
    value    = TRUE,
    evidence = "grain size + mineralogy stable across the interval"
  ),
  depositional_artifact_ruled_out = list(
    value    = TRUE,
    evidence = "continuous varved sequence; no unconformity at the boundary"
  ),
  oipc_ref = sonk_oipc_ref,
  level2_vegetation_path = list(
    vegetation_scenario = list(
      from = c(tree = 0.4, shrub = 0.3, grass = 0.2, C4 = 0.05),
      to   = c(tree = 0.1, shrub = 0.2, grass = 0.5, C4 = 0.20),
      evidence = "hypothetical 30-pp woody-to-grass scenario for the vignette"
    )
  )
)

verdict_sonk_magnitude <- suppressWarnings(assess_claim(
  record = sonk_record,
  claim  = sonk_claim_magnitude
))

c(highest_level = verdict_sonk_magnitude$highest_level,
  l2_passed     = verdict_sonk_magnitude$levels$passed[2])

## ----plot, fig.height = 6, echo = TRUE----------------------------------------
op <- par(mfrow = c(2, 1), mar = c(4.5, 5.4, 2.2, 1), mgp = c(3.4, 0.8, 0))

precip_ylim <- extendrange(
  c(recon_sugan$summary$d2h_precip_lower,
    recon_sugan$summary$d2h_precip_upper,
    recon_sonk$summary$d2h_precip_lower,
    recon_sonk$summary$d2h_precip_upper),
  f = 0.05
)

plot_recon <- function(rec, ages, title, boundary, ylim) {
  ord <- order(ages)
  plot(
    ages[ord],
    rec$summary$d2h_precip_median[ord],
    type = "n",
    xlim = rev(range(ages)),
    ylim = ylim,
    xlab = "Age (yr BP)",
    ylab = expression(delta^2 * H[precipitation] ~ "(‰)"),
    main = title
  )
  polygon(
    c(ages[ord], rev(ages[ord])),
    c(rec$summary$d2h_precip_lower[ord],
      rev(rec$summary$d2h_precip_upper[ord])),
    border = NA, col = adjustcolor("steelblue", alpha.f = 0.3)
  )
  lines(ages[ord], rec$summary$d2h_precip_median[ord],
        type = "o", pch = 16, col = "black")
  abline(v = boundary, lty = 2, col = "red")
}

plot_recon(recon_sugan, sugan$age,
           "Lake Sugan (LS13WASU): small interval contrast",
           boundary = 800, ylim = precip_ylim)
plot_recon(recon_sonk, sonk$age,
           "Sonk11D (LS14LASO): large 4-6 ka contrast",
           boundary = 5000, ylim = precip_ylim)

par(op)