R Under development (unstable) (2025-04-16 r88149 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #### Testing  groupBound() , clusterGroupBound() , etc
> ####          ==========     =================
> 
> ## Examples in  ../man/groupBound.rd  and ../man/clusterGroupBound.rd
> ## should *not*
> ## be duplicated here.  OTOH, want strict testing (*and* be speedy !)
> stopifnot(require(hdi))
Loading required package: hdi
Loading required package: scalreg
Loading required package: lars
Loaded lars 1.3

> 
> ## From Matrix/.../tests-tools-1.R :
> showSys.time <- function(expr, ...) {
+   ## prepend 'Time' for R CMD Rdiff
+   st <- system.time(expr, ...)
+   writeLines(paste("Time", capture.output(print(st))))
+   invisible(st)
+ }
> doExtras <- function ()
+ {
+   interactive() || nzchar(Sys.getenv("R_hdi_check_extra")) ||
+     identical("true", unname(Sys.getenv("R_PKG_CHECKING_doExtras")))
+ }
> do.Extras <- doExtras()
> b64 <- .Machine$sizeof.pointer == 8
> cat(if(b64) "64" else "32", "bit platform; do.Extras = ", do.Extras,"\n")
64 bit platform; do.Extras =  FALSE 
> 
> ##' @title Generate "Reference" Datasets where variables are in *groups*
> ##' @param pvec integer vector of number of variables   p[k], in k-th group, k = 1..np.
> ##' @param s0.v integer vector of non-zero coefficient s0[k], in k-th group, k = 1..np.
> ##' @param n number of observations
> ##' @param ... further arguments such as \code{xtype}, \code{btype}, passed to
> ##'  \code{\link{rXb}}
> ##' @param iteration see \code{\link{rXb}}
> ##' @param do2S ditto
> ##' @return a \code{\link{list}} with components \code{x} and \code{beta} similar to those
> ##' of \code{rXb()} and \code{p.vec} and \code{s0.vec} which are the
> ##' corresponding  arguments.
> ##' @author Martin Maechler
> rXbBlock <- function(p.vec, s0.vec, n, ..., iteration = NA, do2S = FALSE) {
+   stopifnot(is.numeric(p.vec), p.vec == round(p.vec),
+             is.numeric(s0.vec), s0.vec == round(s0.vec),
+             (np <- length(p.vec)) >= 1, length(s0.vec) == np,
+             0 <= s0.vec, s0.vec <= p.vec)
+   grL <-
+     lapply(1:np, function(k)
+       rXb(n = n, p = p.vec[k], s0 = s0.vec[k],
+           ..., iteration=iteration, do2S=do2S))
+   ## now concatenate the X matrices and the beta vectors:
+   list(x    = do.call(cbind, lapply(grL, `[[`, "x")),
+        beta = do.call(c,     lapply(grL, `[[`, "beta")),
+        p.vec = p.vec, s0.vec = s0.vec)
+ }
> 
> ##' @title From (X, beta, sigma),  generate  y = X b + eps
> ##' @param rXb a list, e.g. resulting from rXb (== rXb) or rXbBlock()
> ##' @param sigma a non-negative number
> ##' @return a list with components (x, y)
> ##' @author Martin Maechler
> rXy <- function(rXb, sigma = 1) {
+   stopifnot(is.list(rXb), is.matrix(X <- rXb$x), is.numeric(b <- rXb$beta),
+             ncol(X) == length(b),
+             length(sigma) == 1, is.finite(sigma), 0 <= sigma)
+   list(x = X, y = drop(X %*% b) + rnorm(nrow(X), sd = sigma))
+ }
> 
> ##' from a p-vector, return corresponding list of 'groups', i.e., group indices
> pvec2grp <- function(pv) {
+   stopifnot(is.numeric(pv), is.finite(pv), pv >= 0)
+   cv <- cumsum(c(1,pv))
+   mapply(seq, from=cv[-length(cv)], to = cv[-1L]-1L)
+ }
> 
> if(FALSE) { ## not used currently
+   suppressWarnings(RNGversion("3.5.0"))
+   set.seed(47) # <-> iteration = NA to use R's seed
+   d1 <- rXb(n = 64, p =  16, s0 =  7, iteration=NA, do2S = FALSE)
+   d2 <- rXb(n = 64, p = 256, s0 = 12, iteration=NA, do2S = FALSE)
+   stopifnot(identical(dim(d2$x), c(64L, 256L)), sum(d2$beta != 0) == 12)
+ 
+   yhat.1 <- with(d1, drop(x %*% beta))
+   yhat.2 <- with(d2, drop(x %*% beta))
+ }
> suppressWarnings(RNGversion("3.5.0"))
> set.seed(20)
> d1B <- rXbBlock(p = c(7, 8, 7), s0 = c(2, 4, 1), n = 16)
> p2 <- sample(c(32, 16, 16, 16, 12, 12, 5, 5, 5))
> s2 <- c(rmultinom(1, size = p2, prob = ifelse(p2 > 15, 1/4, 1/8)))
> rbind(p2, s2)
   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
p2   12   16   16    5   16    5    5   32   12
s2    0    2    3    0    2    0    1    2    2
> d2B <- rXbBlock(p = p2, s0 = s2, n = 500)
> d3B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "toeplitz", x.par = 0.95)
> d4B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "equi.corr")
> d5B <- rXbBlock(p = p2, s0 = s2, n = 128, xtype = "exp.decay", x.par = c(0.25, 4))
> 
> if(interactive() && exists("db5", mode="list")) ## save when produced with "new hdi"
+   save(d1B, d2B, d3B, d4B, d5B, file = "rXb_blk.rda")
> if(!exists("rXb", mode="function")) ## "old hdi" -- tests should give the same (?!) -- TODO
+   print(load("rXb_blk.rda"))
> 
> if(dev.interactive(TRUE)) {
+   require(Matrix)
+   ## if you look carefully, you recognize the three blocks of sizes (7, 8, 7):
+   image(as(zapsmall(cov(d1B$x), digits=1), "sparseMatrix"))
+   ## Here, have 9 blocks :
+   image(as(zapsmall(cov(d2B$x), digits=1), "sparseMatrix"))
+   image(as(zapsmall(cov(d3B$x), digits=1), "sparseMatrix")) # easy to see
+   image(as(zapsmall(cov(d4B$x), digits=1), "sparseMatrix")) #  (ditto)
+   image(as(zapsmall(cov(d5B$x), digits=1), "sparseMatrix"))
+                                         # decaying quickly => basically tri-diagonal blocks
+ }
> 
> ##- set.seed(107)
> ##- d1xy <- rXy(d1B, sigma = 0.5)
> ##- showSys.time(
> ##-   gb1 <- groupBound(x = d1xy$x, y = d1xy$y, group = 1:ncol(d1xy$x))
> ##- ) # 2.00 [lynne 2015]
> ##- set.seed(17)# also for random splits
> ##- d2xy <- rXy(d2B, sigma = 1/16)
> ##- showSys.time(
> ##-   gb2  <- groupBound(x = d2xy$x, y = d2xy$y, group = 1:ncol(d2xy$x))
> ##- )
> ##- showSys.time(
> ##-   gb2g <- groupBound(x = d2xy$x, y = d2xy$y, group = pvec2grp(d2B$p.vec),
> ##-                      silent = TRUE)
> ##- )
> ##- 
> ##- stopifnot(
> ##-   as.vector(gb1) == 0
> ##-  ,
> ##-   all.equal(as.vector(gb2),
> ##-             if(b64) 6.529515 else 6.598612, tol = 2e-7)
> ##-  ,
> ##-   all.equal(as.vector(gb2g),
> ##-             if(b64) c(0, 0, 1.924578, 0, 1.436938, 0, 0, 0, 0.3264166)
> ##-             else c(0, 0.01877908, 2.0447942, 0, 1.753219, 0,0,0, 0.3779607),
> ##-             tol = 6e-7)
> ##- )
> 
> if(do.Extras) { ## save time in regular  'R CMD check' ---------
+ 
+   ##--------------- 3 --
+   suppressWarnings(RNGversion("3.5.0"))
+   set.seed(27)# also for random splits
+   d3xy <- rXy(d3B, sigma = 1/16)
+   showSys.time(
+     gb3  <- groupBound(x = d3xy$x, y = d3xy$y, group = 1:ncol(d3xy$x))
+   )
+   stopifnot(
+     all.equal(as.vector(gb3), if(b64) 7.3760622 else 7.7908584, tol = 2e-7)
+   )
+   suppressWarnings(RNGversion("3.5.0"))
+   set.seed(19)# for random splits
+   showSys.time(
+     gb3g <- groupBound(x = d3xy$x, y = d3xy$y, group = pvec2grp(d3B$p.vec),
+                        silent = TRUE)
+   )
+   stopifnot(
+     all.equal(as.vector(gb3g),
+               if(b64) c(0, 0.5237811, 1.354151, 0, 0, 0, 0.8847217, 0, 0)
+               else    c(0, 1.1008886, 1.9249053,0, 0, 0, 0.78661001, 0, 0),
+               tol = 5e-7) # 9.4e-8
+   )
+ 
+   ##--------------- 4 --
+   suppressWarnings(RNGversion("3.5.0"))
+   set.seed(37)# also for random splits
+   d4xy <- rXy(d4B, sigma = 1/16)
+   showSys.time(
+     gb4  <- groupBound(x = d4xy$x, y = d4xy$y, group = 1:ncol(d4xy$x),
+                        silent = TRUE)
+   )
+   showSys.time(
+     gb4g <- groupBound(x = d4xy$x, y = d4xy$y, group = pvec2grp(d4B$p.vec),
+                        silent = TRUE)
+   )
+   stopifnot(
+     all.equal(as.vector(gb4), if(b64) 5.230821 else 5.3752713,  tol = 2e-7)
+     ## all.equal(as.vector(gb4), 8.475868, tol = 2e-7)
+    ,
+     all.equal(as.vector(gb4g),
+               if(b64) c(0, 0, 3.328455, 0, 0, 0, 0, 0, 0)
+               else    c(0, 0.043785673, 3.6154918, 0, 0, 0, 0, 0, 0),
+               tol = 8e-7) # 1.4e-7
+   )
+ 
+ }# only if(do.Extras)
> 
> proc.time()
   user  system elapsed 
   1.50    0.23    1.75