## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  collapse = TRUE, comment = "#>"
)
library(gdpar)
set.seed(1)

## ----suite--------------------------------------------------------------------
suite <- gdpar_geometry_suite()
data.frame(
  target  = names(suite),
  pathology = vapply(suite, `[[`, character(1), "pathology"),
  remedy    = vapply(suite, `[[`, character(1), "geometry_remedy")
)

## ----thresholds---------------------------------------------------------------
str(gdpar_geometry_thresholds())

## ----diag-heavy, eval=FALSE---------------------------------------------------
# # Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
# diag <- gdpar_geometry_diagnostic(suite$G4_quasi_deterministic, n_grid = 3)
# diag$pathology      # "quasi_deterministic"
# diag$culprit        # the localized direction(s)

## ----euclid-------------------------------------------------------------------
P   <- diag(c(1, 100))
tgt <- gdpar_geom_target(
  log_prob      = function(t) -0.5 * drop(t %*% P %*% t),
  grad_log_prob = function(t) -drop(P %*% t), dim = 2L)

fit0 <- gdpar_geom_hmc(tgt, gdpar_geom_metric_euclidean(dim = 2L),
                       n_iter = 400L, n_warmup = 200L, epsilon = 0.12,
                       L = 20L, seed = 1L)
c(accept = round(fit0$accept_rate, 3),
  sd1 = round(sd(fit0$draws[, 1]), 3),   # truth 1
  sd2 = round(sd(fit0$draws[, 2]), 3))   # truth 0.1

## ----dense, eval=FALSE--------------------------------------------------------
# # A dense mass equal to the posterior precision whitens a straight canyon.
# metric_dense <- gdpar_geom_metric_euclidean(M = P)

## ----subriem------------------------------------------------------------------
canyon <- gdpar_geom_target(
  log_prob      = function(th) -0.5 * (th[1]^2 + 50 * th[2]^2),
  grad_log_prob = function(th) -c(th[1], 50 * th[2]), dim = 2L)

metric_sr <- gdpar_geom_metric_subriemannian(
  canyon, fisher = function(th) diag(c(1, 50)))

fit_sr <- gdpar_geom_hmc(canyon, metric = metric_sr, epsilon = 0.5, L = 12L,
                         n_iter = 800L, n_warmup = 300L, seed = 3L)
c(accept = round(fit_sr$accept_rate, 3),
  divergent = fit_sr$n_divergent,
  floor_sd = round(sd(fit_sr$draws[, 1]), 3),   # truth 1
  wall_sd  = round(sd(fit_sr$draws[, 2]), 3))   # truth 1/sqrt(50) ~ 0.141

## ----budget-------------------------------------------------------------------
str(gdpar_geom_orchestrate_budget())
str(gdpar_geom_orchestrate_criteria())

## ----orch-heavy, eval=FALSE---------------------------------------------------
# # Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
# res <- gdpar_geom_orchestrate(suite$G0_isotropic, n_grid = 1)
# res$status     # "resolved" at euclidean_diagonal for the isotropic control

## ----bridge-target------------------------------------------------------------
geom_target <- gdpar_geom_target(
  log_prob      = function(t) -0.5 * sum(t^2),
  grad_log_prob = function(t) -t, dim = 3L)
geom_target$dim
geom_target$grad_log_prob(c(1, 2, 3))

## ----bridge-heavy, eval=FALSE-------------------------------------------------
# # Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
# fit    <- gdpar(gdpar_bf(y ~ a(x), sigma ~ a(z)), data = d,
#                 family = gdpar_family("gaussian"), skip_id_check = TRUE)
# bridge <- gdpar_geom_bridge(fit)
# res    <- gdpar_geom_orchestrate(bridge$target, bridge$geom_target,
#                                  reference = bridge$reference)
# 
# # Or, in one call, building the model through the shared .gdpar_K_build() seam:
# res2 <- gdpar_geom_fit(gdpar_bf(y ~ a(x), sigma ~ a(z)), data = d,
#                        family = gdpar_family("gaussian"), skip_id_check = TRUE)
# res2$status

## ----laplace-good-------------------------------------------------------------
A   <- matrix(c(2, 0.8, 0.8, 1), 2, 2)
mu  <- c(1, -0.5)
tgt <- gdpar_geom_target(
  log_prob      = function(t) -0.5 * as.numeric(t(t - mu) %*% A %*% (t - mu)),
  grad_log_prob = function(t) -as.numeric(A %*% (t - mu)),
  hessian       = function(t) -A, dim = 2L)
lap <- gdpar_geom_laplace(tgt, draws = 500L, seed = 1L)
c(label = lap$fit_quality_label, max_mode_err = max(abs(lap$mode - mu)))

## ----laplace-fallback-heavy, eval=FALSE---------------------------------------
# # Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
# res <- gdpar_geom_orchestrate(bridge$target, bridge$geom_target,
#                               reference = bridge$reference,
#                               laplace_fallback = TRUE, laplace_draws = 2000L)
# res$status                 # "certified_limit_laplace"
# res$laplace$fit_quality_label   # on the real 9.2.O canyon: "very_poor"