## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7)
has_lpsolve <- requireNamespace("lpSolve", quietly = TRUE)

## -----------------------------------------------------------------------------
library(combreg)

A <- rbind(
  c(1, 1, 0),
  c(0, 1, 1)
)
con <- crr_constraints(A, b = c(1, 1))
con

## -----------------------------------------------------------------------------
is_tum(A)                                     # TRUE: this A is TUM
is_tum(rbind(c(1, 1, 0), c(1, 0, 1), c(0, 1, 1)))  # FALSE (odd cycle)

## -----------------------------------------------------------------------------
candidates <- rbind(
  c(1, 0, 1),   # A y = (1, 1) <= (1, 1): feasible
  c(1, 1, 0)    # A y = (2, 1): violates row 1
)
is_feasible(con, candidates)

## -----------------------------------------------------------------------------
d <- 3
A_simplex <- matrix(1, nrow = 1, ncol = d)
con_simplex <- crr_constraints(A_simplex, b = 1)
con_simplex
is_feasible(con_simplex, rbind(diag(d)[1, ], rep(1, d)))  # e_1 ok, all-ones not

## ----eval = has_lpsolve-------------------------------------------------------
sim <- simulate_crr(n = 60, p = 2, constraints = con, seed = 1)
head(sim$Y)

## ----eval = has_lpsolve-------------------------------------------------------
fit <- crr(sim$Y, sim$X, con,
           kernel = "exponential",
           n_iter = 400, warmup = 200, seed = 1)
fit

## ----eval = has_lpsolve-------------------------------------------------------
head(summary(fit))
plot(fit, pars = c("beta[1,1]", "beta[2,1]"))

## ----eval = has_lpsolve-------------------------------------------------------
fit_un <- crr(sim$Y, sim$X, method = "unconstrained",
              n_iter = 400, warmup = 200, seed = 1)

rmse <- function(est, truth) sqrt(mean((est - truth)^2))
c(constrained   = rmse(coef(fit), sim$beta),
  unconstrained = rmse(coef(fit_un), sim$beta))

## ----eval = has_lpsolve-------------------------------------------------------
fit$zeta_block_tuned

## ----eval = has_lpsolve-------------------------------------------------------
fit_fixed <- crr(sim$Y, sim$X, con,
                 n_iter = 300, warmup = 150, seed = 1,
                 control = crr_control(zeta_block = 100))
fit_fixed$zeta_block_tuned   # block sizes are capped at d

## ----eval = has_lpsolve-------------------------------------------------------
fit_custom <- crr(sim$Y, sim$X, con,
                  kernel = "half_gaussian",         # half-Gaussian dual kernel
                  prior  = crr_prior(sd = 10),      # wider N(0, 100) prior
                  n_iter = 400, warmup = 200,       # total / discarded sweeps
                  thin   = 2,                       # keep every 2nd draw
                  chains = 2,                       # independent chains (serial)
                  seed   = 1,                        # reproducible RNG stream
                  control = crr_control(
                    n_iter_hit_and_run = 20,        # inner hit-and-run steps
                    zeta_block = "adaptive",        # auto-tune the block size
                    n_threads = 1))                 # OpenMP threads (result-invariant)
fit_custom

## ----eval = has_lpsolve-------------------------------------------------------
df <- data.frame(x1 = sim$X[, 1], x2 = sim$X[, 2])
fit_formula <- crr(sim$Y, ~ 0 + x1 + x2, con, data = df,
                   n_iter = 300, warmup = 150, seed = 1)
coef(fit_formula)