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

## ----setup, eval=FALSE--------------------------------------------------------
# library(AdapDiscom)
# library(MASS)

## ----eval=FALSE---------------------------------------------------------------
# # Number of variables
# p <- 10
# 
# # AR(1) structure
# Sigma_ar1 <- generate.cov(p, example = 1)
# print(Sigma_ar1[1:3, 1:3])
# 
# # Block diagonal structure
# Sigma_block <- generate.cov(p, example = 2)
# print(Sigma_block[1:3, 1:3])
# 
# # Kronecker product structure
# Sigma_kron <- generate.cov(p, example = 3)
# print(dim(Sigma_kron))

## ----eval=FALSE---------------------------------------------------------------
# set.seed(123)
# # Set parameters
# p <- 300
# n <- 440
# n.tuning <- 200
# n.test <- 400
# p1 <- p%/%3
# p2 <- p%/%3
# p3 <- p - p1 - p2
# pp <- c(p1, p2, p3)  # Block sizes
# sigma <- 1

## ----eval=FALSE---------------------------------------------------------------
# # Generate different covariance matrices for each block
# cov.mat1 <- generate.cov(p1, 1)  # AR(1) structure for block 1
# cov.mat2 <- generate.cov(p2, 2)  # Block diagonal for block 2
# cov.mat3 <- generate.cov(p3, 3)  # Kronecker product for block 3

## ----eval=FALSE---------------------------------------------------------------
# # Generate true beta coefficients
# beta1 <- rep(c(rep(0.5, 5), rep(0, 95)), p/100)
# beta.true <- beta1
# 
# # Generate complete data for all samples
# pre.x1 <- mvrnorm(n + n.tuning + n.test, rep(0, p%/%3), cov.mat1)
# pre.x2 <- mvrnorm(n + n.tuning + n.test, rep(0, p%/%3), cov.mat2)
# pre.x3 <- mvrnorm(n + n.tuning + n.test, rep(0, p%/%3), cov.mat3)
# pre.x <- cbind(pre.x1, pre.x2, pre.x3)
# pre.ep <- rnorm(n + n.tuning + n.test, 0, sigma)
# 
# # Training data
# n.com <- n/4
# n1 <- n2 <- n3 <- n4 <- n.com
# X.train <- pre.x[1:n, ]
# ep <- pre.ep[1:n]
# y.train <- X.train %*% beta1 + ep
# colnames(X.train) <- paste0("X", 1:p)
# 
# # Introduce missing data pattern
# X.train[(n1 + 1):(n1 + n2), (p1 + p2 + 1):(p1 + p2 + p3)] <- NA
# X.train[(n1 + n2 + 1):(n1 + n2 + n3), (p1 + 1):(p1 + p2)] <- NA
# X.train[(n1 + n2 + n3 + 1):(n1 + n2 + n3 + n4), (1:p1)] <- NA
# 
# # Tuning data
# X.tuning <- pre.x[(n + 1):(n + n.tuning), ]
# ep.tuning <- pre.ep[(n + 1):(n + n.tuning)]
# y.tuning <- X.tuning %*% beta1 + ep.tuning
# colnames(X.tuning) <- paste0("X", 1:p)
# 
# # Test data
# X.test <- pre.x[(n + n.tuning + 1):(n + n.tuning + n.test), ]
# ep.test <- pre.ep[(n + n.tuning + 1):(n + n.tuning + n.test)]
# y.test <- X.test %*% beta1 + ep.test
# colnames(X.test) <- paste0("X", 1:p)

## ----eval=FALSE---------------------------------------------------------------
# # Run AdapDiscom with default parameters
# result <- adapdiscom(
#   beta = beta.true,      # True coefficients (optional, for evaluation)
#   x = X.train,          # Training data
#   y = y.train,          # Training response
#   x.tuning = X.tuning,  # Tuning data
#   y.tuning = y.tuning,  # Tuning response
#   x.test = X.test,      # Test data
#   y.test = y.test,      # Test response
#   nlambda = 20,         # Number of lambda values
#   nalpha = 10,          # Number of alpha values
#   pp = pp,              # Block sizes
#   robust = 0,           # Classical estimation
#   standardize = TRUE,   # Standardize data
#   itcp = TRUE          # Include intercept
# )
# 
# # View results
# print(paste("R-squared:", round(result$R2, 4)))

## ----eval=FALSE---------------------------------------------------------------
# # Optimal parameters
# cat("Optimal lambda:", result$lambda, "\n")
# cat("Optimal alpha:", paste(round(result$alpha, 4), collapse = ", "), "\n")
# 
# # Performance metrics
# cat("Training error:", round(result$train.error, 4), "\n")
# cat("Test error:", round(result$test.error, 4), "\n")
# cat("R-squared:", round(result$R2, 4), "\n")
# 
# # Model selection
# cat("Number of selected variables:", result$select, "\n")
# 
# # If true beta provided, evaluation metrics
# if (!is.null(beta.true)) {
#   cat("Estimation error:", round(result$est.error, 4), "\n")
#   cat("False Positive Rate:", round(result$fpr, 4), "\n")
#   cat("False Negative Rate:", round(result$fnr, 4), "\n")
# }
# 
# # Estimated coefficients
# cat("First 6 estimated coefficients:\n")
# print(round(head(result$a1), 4))

## ----eval=FALSE---------------------------------------------------------------
# # Compare different methods
# methods <- c("adapdiscom", "discom", "fast_adapdiscom", "fast_discom")
# results <- list()
# 
# # AdapDiscom
# results$adapdiscom <- result # Already computed above
# 
# # DISCOM
# results$discom <- discom(
#   beta = beta.true, x = X.train, y = y.train,
#   x.tuning = X.tuning, y.tuning = y.tuning,
#   x.test = X.test, y.test = y.test,
#   nlambda = 20, nalpha = 10, pp = pp
# )
# 
# # Fast AdapDiscom
# results$fast_adapdiscom <- fast_adapdiscom(
#   beta = beta.true, x = X.train, y = y.train,
#   x.tuning = X.tuning, y.tuning = y.tuning,
#   x.test = X.test, y.test = y.test,
#   nlambda = 20, nalpha = 10, pp = pp
# )
# 
# # Fast DISCOM
# results$fast_discom <- fast_discom(
#   beta = beta.true, x = X.train, y = y.train,
#   x.tuning = X.tuning, y.tuning = y.tuning,
#   x.test = X.test, y.test = y.test,
#   nlambda = 20, pp = pp
# )
# 
# # Compare performance
# comparison <- data.frame(
#   Method = names(results),
#   Test_Error = round(sapply(results, function(x) x$test.error), 4),
#   R_Squared = round(sapply(results, function(x) x$R2), 4),
#   Selected_Vars = sapply(results, function(x) x$select),
#   Est_Error = round(sapply(results, function(x) x$est.error), 4),
#   FPR = round(sapply(results, function(x) x$fpr), 4),
#   FNR = round(sapply(results, function(x) x$fnr), 4),
#   Time = round(sapply(results, function(x) x$time), 2)
# )
# 
# print(comparison)

## ----eval=FALSE---------------------------------------------------------------
# # Run with robust estimation
# result_robust <- adapdiscom(
#   beta = beta.true,
#   x = X.train, y = y.train,
#   x.tuning = X.tuning, y.tuning = y.tuning,
#   x.test = X.test, y.test = y.test,
#   nlambda = 20, nalpha = 10, pp = pp,
#   robust = 1,        # Enable robust estimation
#   k.value = 1.5      # Huber tuning parameter
# )

## ----eval=FALSE---------------------------------------------------------------
# # Advanced usage with custom parameters
# result_advanced <- adapdiscom(
#   beta = beta.true,
#   x = X.train, y = y.train,
#   x.tuning = X.tuning, y.tuning = y.tuning,
#   x.test = X.test, y.test = y.test,
#   nlambda = 30,              # More lambda values
#   nalpha = 15,               # More alpha values
#   pp = pp,
#   robust = 0,
#   standardize = TRUE,
#   itcp = TRUE,
#   lambda.min.ratio = 0.001,  # Custom lambda range
#   k.value = 1.5
# )