R Under development (unstable) (2024-09-23 r87189 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> 
> library(polyapost)
Loading required package: rcdd
If you want correct answers, use rational arithmetic.
See the Warnings sections in help pages for
    functions that do computational geometry.
> 
> set.seed(42)
> 
> d <- 10
> r <- floor((d - 1) / 4)
> 
> i <- 1:d
> a2 <- rbind(i, as.numeric(i <= (d + 1) / 2) - as.numeric(i >= (d + 1) / 2))
> b2 <- c((d + 1) / 2, 0)
> a1 <- rbind(- as.numeric(i <= (d + 1) / 2 - r) -
+     as.numeric(i >= (d + 1) / 2 + r))
> b1 <- - 1 / 2
> dimnames(a1) <- NULL
> dimnames(a2) <- NULL
> 
> hout <- hitrun(rep(1, d), nbatch = 100,
+     a1 = a1, b1 = b1, a2 = a2, b2 = b2, debug = TRUE)
> names(hout)
 [1] "batch"        "initial"      "final"        "current"      "proposal"    
 [6] "z"            "u1"           "u2"           "s1"           "s2"          
[11] "log.green"    "initial.seed" "final.seed"   "time"         "alpha"       
[16] "nbatch"       "blen"         "nspac"        "origin"       "basis"       
[21] "amat"         "bvec"         "debug"        "a1"           "b1"          
[26] "a2"           "b2"           "split.time"  
> 
> # check random numbers
> set.seed(42)
> identical(.Random.seed, hout$initial.seed)
[1] TRUE
> myz <- matrix(NA_real_, nrow(hout$current), ncol(hout$current))
> myu1 <- rep(NA_real_, nrow(hout$current))
> myu2 <- rep(NA_real_, nrow(hout$current))
> for (i in seq(along = myu1)) {
+     myz[i, ] <- rnorm(ncol(hout$current))
+     myu1[i] <- runif(1)
+ }
> identical(.Random.seed, hout$final.seed)
[1] TRUE
> identical(myz, hout$z)
[1] TRUE
> identical(myu1, hout$u1)
[1] TRUE
> identical(myu2, hout$u2)
[1] TRUE
> identical(.Random.seed, hout$final.seed)
[1] TRUE
> 
> # check green ratios (trivial for this problem)
> all(hout$log.green == 0)
[1] TRUE
> 
> # check construction of basis, origin, amat, bvec, rip
> hrep4 <- makeH(a1 = a1, b1 = b1)
> hrep4 <- addHin(- diag(d), rep(0, d), hrep4)
> hrep3 <- makeH(a2 = a2, b2 = b2)
> hrep3 <- addHeq(rep(1, d), 1, hrep3)
> hrep4 <- d2q(hrep4)
> hrep3 <- d2q(hrep3)
> dim(hrep3)
[1]  3 12
> dim(hrep4)
[1] 11 12
> vrep3 <- scdd(hrep3)$output
> is.line <- vrep3[ , 1] == "1" & vrep3[ , 2] == "0"
> is.point <- vrep3[ , 1] == "0" & vrep3[ , 2] == "1"
> all(is.point | is.line)
[1] TRUE
> sum(is.point) == 1
[1] TRUE
> sum(is.line) == d - nrow(hrep3)
[1] TRUE
> foo <- vrep3[ , - c(1, 2)]
> origin <- foo[is.point, ]
> basis <- foo[is.line, ]
> basis <- t(basis)
> basis <- qgram(basis)
> identical(q2d(origin), hout$origin)
[1] TRUE
> identical(q2d(basis), hout$basis)
[1] TRUE
> 
> amat <- qneg(hrep4[ , - c(1, 2)])
> bvec <- hrep4[ , 2]
> bvec <- qmq(bvec, qmatmult(amat, cbind(origin)))
> amat <- qmatmult(amat, basis)
> identical(q2d(amat), hout$amat)
[1] TRUE
> identical(q2d(bvec), hout$bvec)
[1] TRUE
> 
> 
> origin <- q2d(origin)
> basis <- q2d(basis)
> amat <- q2d(amat)
> bvec <- q2d(bvec)
> initial <- hout$initial
> all(qsign(qmq(bvec, as.vector(qmatmult(amat, cbind(initial))))) > 0)
[1] TRUE
> 
> # check bounds
> mys1 <- rep(Inf, nrow(hout$current))
> mys2 <- rep(-Inf, nrow(hout$current))
> for (i in seq(along = myu1)) {
+    z <- hout$z[i, ]
+    x <- hout$current[i, ]
+    ax <- as.numeric(amat %*% x)
+    az <- as.numeric(amat %*% z)
+    bnd <- (bvec - ax) / az
+    mys1[i] <- max(bnd[az < 0])
+    mys2[i] <- min(bnd[az > 0])
+ }
> all.equal(mys1, hout$s1)
[1] TRUE
> all.equal(mys2, hout$s2)
[1] TRUE
> 
> # check proposal
> myproposal <- matrix(NA, nrow(hout$current), ncol(hout$current))
> for (i in seq(along = myu1)) {
+    x <- hout$current[i, ]
+    z <- hout$z[i, ]
+    smin <- hout$s1[i]
+    smax <- hout$s2[i]
+    u <- hout$u1[i]
+    myproposal[i, ] <- x + z * (u * smin + (1 - u) * smax)
+ }
> all.equal(myproposal, hout$proposal)
[1] TRUE
> identical(hout$current[- 1, ], hout$proposal[- nrow(hout$proposal), ])
[1] TRUE
> identical(hout$current[1, ], hout$initial)
[1] TRUE
> identical(hout$proposal[nrow(hout$proposal), ], hout$final)
[1] TRUE
> 
> # check path is feasible (reduced coordinates)
> foo <- hout$current %*% t(amat)
> foo <- sweep(foo, 2, bvec)
> all(foo <= 0)
[1] TRUE
> 
> # check transformation
> foo <- hout$proposal %*% t(basis)
> foo <- sweep(foo, 2, origin, "+")
> all.equal(foo, hout$batch, tol = 1e-13)
[1] TRUE
> 
> # check path is feasible (original coordinates)
> foo <- hout$batch %*% t(a1)
> foo <- sweep(foo, 2, b1)
> bar <- hout$batch %*% t(a2)
> bar <- sweep(bar, 2, b2)
> all(foo <= sqrt(.Machine$double.eps))
[1] TRUE
> all(abs(bar) <= sqrt(.Machine$double.eps))
[1] TRUE
> 
> # now for non-uniform
> 
> alpha <- rep(1.3, d)
> ludfun <- function(x) {
+     if (any(x <= 0)) return(-Inf)
+     return(sum((alpha - 1) * log(x)))
+ }
> 
> hout <- hitrun(alpha, nbatch = 100,
+     a1 = a1, b1 = b1, a2 = a2, b2 = b2, debug = TRUE)
> 
> foo <- hout$current %*% t(basis)
> foo <- sweep(foo, 2, origin, "+")
> bar <- hout$proposal %*% t(basis)
> bar <- sweep(bar, 2, origin, "+")
> my.log.green <- apply(bar, 1, ludfun) - apply(foo, 1, ludfun)
> all.equal(my.log.green, hout$log.green, tol = 1e-13)
[1] TRUE
> identical(is.na(hout$u2), hout$log.green >= 0)
[1] TRUE
> my.accept.1 <- is.na(hout$u2) | hout$u2 < exp(hout$log.green)
> foo <- hout$proposal
> bar <- hout$current[- 1, ]
> bar <- rbind(bar, hout$final)
> baz <- foo - bar
> my.accept.2 <- apply(baz == 0, 1, all)
> identical(my.accept.1, my.accept.2)
[1] TRUE
> 
> # now check restart property
> 
> .Random.seed <- hout$initial.seed
> hout1 <- hitrun(alpha, nbatch = 50, a1 = a1, b1 = b1, a2 = a2, b2 = b2)
> hout2 <- hitrun(hout1)
> identical(hout$batch, rbind(hout1$batch, hout2$batch))
[1] TRUE
> identical(hout$final, hout2$final)
[1] TRUE
> identical(hout$final.seed, hout2$final.seed)
[1] TRUE
> 
> # now check batching and spacing
> 
> hout3 <- hitrun(hout, nbatch = 17, nspac = 3, blen = 11)
> hout4 <- hitrun(hout, nspac = 1, blen = 1,
+     nbatch = hout3$nbatch * hout3$blen * hout3$nspac)
> mybatch <- hout4$batch[seq(1, hout4$nbatch) %% hout3$nspac == 0, ]
> dim(mybatch)
[1] 187  10
> mybatch <- array(as.vector(mybatch),
+     c(hout3$blen, hout3$nbatch, ncol(mybatch)))
> dim(mybatch)
[1] 11 17 10
> mybatch <- apply(mybatch, c(2, 3), mean)
> dim(mybatch)
[1] 17 10
> all.equal(mybatch, hout3$batch)
[1] TRUE
> 
> # now check with outmat
> 
> i <- 1:d
> outmat <- rbind(i, i^2)
> dimnames(outmat) <- NULL
> 
> hout <- hitrun(alpha, nbatch = 101, blen = 17, outmat = outmat,
+     a1 = a1, b1 = b1, a2 = a2, b2 = b2, debug = TRUE)
> nrow(outmat) == ncol(hout$batch)
[1] TRUE
> mynext <- rbind(hout$current[- 1, ], hout$final)
> 
> foo <- mynext %*% t(basis)
> foo <- sweep(foo, 2, origin, "+")
> foo <- foo %*% t(outmat)
> foo <- array(as.vector(foo), c(hout$blen, hout$nbatch, ncol(foo)))
> foo <- apply(foo, c(2, 3), mean)
> all.equal(foo, hout$batch)
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   0.64    0.06    0.68