## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  collapse = TRUE, comment = "#>"
)
have_cmdstan <- requireNamespace("cmdstanr", quietly = TRUE)

## ----build-uniform, eval=TRUE-------------------------------------------------
library(gdpar)

spec_uniform <- amm_spec(
  p    = 2L,
  dims = dimwise(a = ~ x1 + x2, b = ~ x1),
  W    = W_basis(type = "polynomial", degree = 2)
)
print(spec_uniform)

## ----build-override, eval=TRUE------------------------------------------------
spec_mixed <- amm_spec(
  p    = 3L,
  dims = dimwise(a = ~ x1 + x2, b = ~ x1) |>
           override(k = 2L, a = ~ x1) |>
           override(k = 3L, b = NULL),
  W    = W_basis(type = "polynomial", degree = 2)
)
print(spec_mixed)

## ----build-chainable, eval=TRUE-----------------------------------------------
spec_chain <- amm_build(p = 3L) |>
  amm_set_a_uniform(~ x1 + x2) |>
  amm_set_a(k = 2L, ~ x1) |>
  amm_set_b_uniform(~ x1) |>
  amm_set_b(k = 3L, NULL) |>
  amm_set_W(W_basis(type = "polynomial", degree = 2)) |>
  amm_set_x_vars(c("x1", "x2")) |>
  as_amm_spec()
print(spec_chain)

## ----serialize, eval=TRUE-----------------------------------------------------
tmp <- tempfile(fileext = ".gdpar")
amm_save_spec(spec_uniform, tmp)
cat(readLines(tmp), sep = "\n")

spec_roundtrip <- amm_load_spec(tmp)
identical(spec_uniform$level, spec_roundtrip$level)

## ----families, eval=TRUE------------------------------------------------------
fam_mv1 <- gdpar_family_multi("gaussian", p = 2L)
print(fam_mv1)

fam_mv2 <- gdpar_family_multi(gdpar_family("poisson"), p = 3L)
print(fam_mv2)

fam_mv3 <- gdpar_family_multi(
  list(
    gdpar_family("gaussian"),
    gdpar_family("gaussian")
  ),
  p = 2L
)
print(fam_mv3)

## ----promotion-illustrative, eval=FALSE---------------------------------------
# fit <- gdpar(
#   formula = y ~ x1 + x2,
#   family  = gdpar_family("gaussian"),  # auto-promoted to gdpar_family_multi
#   amm     = amm_spec(p = 2L,
#                      dims = dimwise(a = ~ x1 + x2)),
#   data    = df
# )
# # Informational message emitted (class gdpar_family_promotion_message):
# #   "Auto-promoted univariate family 'gaussian' to gdpar_family_multi
# #    via gdpar_family_multi(family, p = 2) because amm$p > 1.
# #    Pass an explicit gdpar_family_multi to silence."

## ----W-materialize, eval=TRUE-------------------------------------------------
wb_poly_p2 <- W_basis(type = "polynomial", degree = 2, p = 2L)
print(wb_poly_p2)
wb_poly_p2$block_indices

wb_bs <- W_basis(type = "bspline", degree = 3, df = 4, p = 2L)
print(wb_bs)
wb_bs$block_indices

## ----W-per-k, eval=TRUE-------------------------------------------------------
subs <- as_per_k(wb_poly_p2)
length(subs)
subs[[1L]]$dim
subs[[2L]]$dim

## ----outcome-illustrative, eval=FALSE-----------------------------------------
# df$y <- cbind(y1, y2)              # matrix column added to existing df
# # or
# df <- data.frame(x1 = ..., x2 = ...)
# df$y <- y_matrix                   # y_matrix is a numeric matrix n x p

## ----gdpar-illustrative, eval=FALSE-------------------------------------------
# fit <- gdpar(
#   formula                     = y ~ x1 + x2,
#   family                      = gdpar_family_multi("gaussian", p = 2L),
#   amm                         = spec_uniform,
#   data                        = df,                # df$y is n x 2 matrix column
#   parametrization             = "auto",            # runs per-coord preflight
#   parametrization_aggregation = "any_ncp",         # default; conservative
#   iter_warmup                 = 1000L,
#   iter_sampling               = 1000L,
#   chains                      = 4L,
#   seed                        = 42L
# )

## ----preflight-shape, eval=FALSE----------------------------------------------
# # After the fit:
# rep <- fit$parametrization$report   # gdpar_preflight_report
# print(rep)                           # default: level = "global"
# print(rep, level = "dim")            # per-coordinate detail
# print(rep, level = "both")           # global + per-coord
# 
# preflight_per_dim(rep)               # tidy data.frame, one row per (component, dim)
# preflight_global_decision(rep)       # one row per component
# as.data.frame(rep)                   # same as preflight_per_dim()
# summary(rep)                         # named list with aggregates

## ----predict-illustrative, eval=FALSE-----------------------------------------
# arr <- predict(fit)                                 # type = "theta_i", summary = "draws"
# dim(arr)                                            # (S, n, p)
# 
# arr_resp <- predict(fit, type = "response")         # per-coord inverse link applied
# arr_se   <- predict(fit, summary = "mean_se")       # list of p data.frames
# arr_q    <- predict(fit, summary = "quantiles")     # list of p data.frames
# 
# arr_new  <- predict(fit, newdata = df_new)          # reconstruct from posterior draws

## ----coef-illustrative, eval=FALSE--------------------------------------------
# cf <- coef(fit)
# print(cf)                                       # level = "global", default
# print(cf, level = "coord")                      # per-coord means only
# print(cf, level = "full")                       # per-coord with full quantiles
# 
# summary(cf)                                     # aggregated counts + theta_ref mean
# format(cf)                                      # one-line representation

## ----asdf-illustrative, eval=FALSE--------------------------------------------
# df_coef <- as.data.frame(coef(fit))
# head(df_coef)
# 
# # Using dplyr with explicit namespace to avoid library() in vignettes:
# df_coef |>
#   dplyr::filter(component == "a") |>
#   dplyr::group_by(k) |>
#   dplyr::summarise(mean_abs = mean(abs(mean)),
#                    max_q95  = max(q95))

## ----smoke-data, eval=have_cmdstan--------------------------------------------
library(gdpar)

set.seed(42L)
n <- 300L
df <- data.frame(
  x1 = rnorm(n),
  x2 = rnorm(n)
)
# True theta_ref = c(0.5, -0.5), beta_a 2x2, sigma_y = 0.3 per coord
true_theta_ref <- c(0.5, -0.5)
true_beta_a <- matrix(c(0.8, -0.6,
                        0.4,  0.7), nrow = 2L, byrow = TRUE)
y_mat <- matrix(NA_real_, nrow = n, ncol = 2L)
for (k in seq_len(2L)) {
  y_mat[, k] <- true_theta_ref[k] +
                true_beta_a[k, 1L] * df$x1 +
                true_beta_a[k, 2L] * df$x2 +
                rnorm(n, sd = 0.3)
}
df$y <- y_mat
str(df)

## ----smoke-spec, eval=have_cmdstan--------------------------------------------
spec <- amm_spec(
  p    = 2L,
  dims = dimwise(a = ~ x1 + x2)
)
print(spec)

## ----smoke-fit, eval=have_cmdstan---------------------------------------------
fit <- gdpar(
  formula                     = y ~ x1 + x2,
  family                      = gdpar_family("gaussian"),  # auto-promoted
  amm                         = spec,
  data                        = df,
  parametrization             = "auto",
  parametrization_aggregation = "any_ncp",
  iter_warmup                 = 300L,
  iter_sampling               = 300L,
  chains                      = 2L,
  refresh                     = 0L,
  verbose                     = FALSE,
  seed                        = 42L
)
print(fit)

## ----smoke-params, eval=have_cmdstan------------------------------------------
fit$parametrization$cp_a
fit$parametrization$cp_W
fit$parametrization$cp_a_per_k
fit$parametrization$cp_W_per_k

## ----smoke-preflight, eval=have_cmdstan---------------------------------------
rep <- fit$parametrization$report
print(rep, level = "both")

## ----smoke-coef, eval=have_cmdstan--------------------------------------------
cf <- coef(fit)
print(cf, level = "coord")

df_coef <- as.data.frame(cf)
head(df_coef, 8L)

## ----smoke-predict, eval=have_cmdstan-----------------------------------------
q_list <- predict(fit, type = "response", summary = "quantiles")
length(q_list)
head(q_list$dim_1, 4L)
head(q_list$dim_2, 4L)

## ----smoke-diag, eval=have_cmdstan--------------------------------------------
diagnostics(fit)

## ----reporter-csv, eval=FALSE-------------------------------------------------
# csv_path <- system.file(
#   "benchmarks", "results", "cp_ncp_hit_rate_multi.csv",
#   package = "gdpar"
# )
# results <- utils::read.csv(csv_path, stringsAsFactors = FALSE)
# str(results)

## ----reporter-agg, eval=FALSE-------------------------------------------------
# sub <- results[!is.na(results$regime_truth) & !is.na(results$hit), ]
# hits_comp <- aggregate(hit ~ component, data = sub,
#                         FUN = function(x) mean(as.logical(x)))
# hits_sc <- aggregate(hit ~ scenario, data = sub,
#                       FUN = function(x) mean(as.logical(x)))
# print(hits_comp)
# print(hits_sc)

## ----reporter-plot, eval=FALSE------------------------------------------------
# # Using ggplot2 with explicit namespace to avoid library() in vignettes;
# # the actual reporter at inst/benchmarks/scripts/report_hit_rate_multi.R
# # follows the same convention.
# results$hit_label <- ifelse(
#   is.na(results$hit), "borderline",
#   ifelse(as.logical(results$hit), "hit", "miss")
# )
# ggplot2::ggplot(
#   results,
#   ggplot2::aes(x = factor(k), fill = hit_label)
# ) +
#   ggplot2::geom_bar(width = 0.7) +
#   ggplot2::facet_grid(scenario ~ component, scales = "free_x",
#                        space = "free_x") +
#   ggplot2::scale_fill_manual(values = c(
#     "hit"        = "#2c7bb6",
#     "miss"       = "#d7191c",
#     "borderline" = "#fdae61"
#   )) +
#   ggplot2::labs(x = "k (coordinate)", y = "count", fill = "verdict") +
#   ggplot2::theme_minimal()

## ----psis-loo-subject, eval = FALSE-------------------------------------------
# fit <- gdpar(
#   formula  = y ~ x1 + x2,
#   family   = gdpar_family_multi("gaussian", p = 2L),
#   amm      = amm_spec(p = 2L, dims = dimwise(a = ~ x1 + x2, b = NULL)),
#   data     = train_df,
#   refresh  = 0L, seed = 42L
# )
# lo <- gdpar_loo(fit)
# print(lo)
# # elpd_loo, se_elpd_loo, Pareto-k summary as in loo::loo

## ----psis-loo-cell, eval = FALSE----------------------------------------------
# lo_cell <- gdpar_loo(fit, aggregation = "cell")
# sum(lo_cell$diagnostics$pareto_k > 0.7) # per-cell concentration