R Under development (unstable) (2023-04-23 r84305 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

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(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)
> 
>  ways <- 7
> 
>  dat <- matrix(NA, nrow = 2^ways, ncol = ways)
>  for (i in 1:ways)
+     dat[ , i] = rep(rep(0:1, each = 2^(i - 1)), times = 2^(ways - i))
> 
>  colnames(dat) <- paste("v", 1:ways, sep = "")
>  dat <- as.data.frame(dat)
>  for (i in 1:ncol(dat)) dat[[i]] <- as.factor(dat[[i]])
>  
>  mu <- 5
>  y <- rpois(nrow(dat), mu)
>  dat <- cbind(dat, y = y)
> 
>  M <- model.matrix(y ~ (v1 + v2 + v3 + v4 + v5 + v6 + v7)^3, data = dat)
> 
>  v <- M
>  linearity <- y > 0
> 
>    stopifnot(is.numeric(v))
>    stopifnot(all(is.finite(v)))
>    stopifnot(is.matrix(v))
>    if (! missing(linearity)) {
+        stopifnot(is.logical(linearity))
+        stopifnot(length(linearity) == nrow(v))
+    } else {
+        linearity <- rep(FALSE, nrow(v))
+    }
> 
>    v <- d2q(v)
>    lresult <- rep(TRUE, nrow(v))
> 
>        vresult <- v[lresult, , drop = FALSE]
>        w <- apply(vresult, 2, qsum)
>        wminus <- qmq(rep("0", length(w)), w)
> 
>        hrep <- rbind(wminus, vresult)
>        fred <- c(1, rep(0, nrow(vresult)))
>        hrep <- cbind(as.character(fred), hrep)
>        fred <- c(0, as.numeric(linearity[lresult]))
>        hrep <- cbind(as.character(fred), hrep)
>        dimnames(hrep) <- NULL
> 
>  out <- lpcdd(hrep, w, minimize = FALSE)
>  print(out)
$solution.type
[1] "Optimal"

$primal.solution
 [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
[20] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
[39] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
[58] "0" "0" "0" "0" "0" "0" "0"

$dual.solution
  [1] "0"     "0"     "-4/5"  "-4/5"  "0"     "-4/5"  "0"     "0"     "-24/5"
 [10] "-4/5"  "0"     "0"     "-24/5" "0"     "-24/5" "-24/5" "32/5"  "-4"   
 [19] "0"     "0"     "-8/5"  "0"     "-8/5"  "-8/5"  "0"     "0"     "-8/5" 
 [28] "-8/5"  "0"     "-8/5"  "0"     "0"     "-12/5" "-4"    "0"     "0"    
 [37] "-8/5"  "0"     "-8/5"  "-8/5"  "0"     "0"     "-8/5"  "-8/5"  "0"    
 [46] "-8/5"  "0"     "0"     "-12/5" "0"     "-8/5"  "-8/5"  "0"     "-8/5" 
 [55] "0"     "0"     "-12/5" "-8/5"  "0"     "0"     "-12/5" "0"     "-12/5"
 [64] "-12/5" "0"     "-8/5"  "0"     "0"     "-12/5" "0"     "-12/5" "-12/5"
 [73] "0"     "0"     "-12/5" "-12/5" "0"     "-12/5" "0"     "0"     "0"    
 [82] "0"     "-12/5" "-12/5" "0"     "-12/5" "0"     "0"     "0"     "-12/5"
 [91] "0"     "0"     "0"     "0"     "0"     "0"     "-32/5" "0"     "-12/5"
[100] "-12/5" "0"     "-12/5" "0"     "0"     "0"     "-12/5" "0"     "0"    
[109] "0"     "0"     "0"     "0"     "-32/5" "4/5"   "0"     "0"     "-16/5"
[118] "0"     "-16/5" "-16/5" "0"     "0"     "-16/5" "-16/5" "0"     "-16/5"
[127] "0"     "0"     "12/5" 

$optimal.value
[1] "0"

> 
>  ##### check gradient of Lagrangian function zero
>  b.augmented <- hrep[ , 2]
>  v.augmented <- hrep[ , - c(1, 2)]
>  blurfle <- out$dual.solution
>  all(qpq(w, qmatmult(rbind(blurfle), v.augmented)) == "0")
[1] TRUE
> 
>  ##### check primal feasibility (trivial here since solution is zero, but)
>  blurfle <- out$primal.solution
>  foo <- qpq(b.augmented, qmatmult(v.augmented, cbind(blurfle)))
>  all(qsign(foo) >= 0)
[1] TRUE
>  
>  ##### check dual feasibility
>  blurfle <- qsign(out$dual.solution)
>  linearity.augmented <- c(FALSE, linearity)
>  length(linearity.augmented) == nrow(hrep)
[1] TRUE
>  all(blurfle[! linearity.augmented] >= 0)
[1] TRUE
> 
>  ##### check complementary slackness
>  all(qsign(foo) * qsign(blurfle) == 0)
[1] TRUE
> 
>  ##### now redo with ordinary computer arithmetic
>  aout <- lpcdd(q2d(hrep), q2d(w), minimize = FALSE)
>  names(aout)
[1] "solution.type"   "primal.solution" "dual.solution"   "optimal.value"  
>  names(out)
[1] "solution.type"   "primal.solution" "dual.solution"   "optimal.value"  
>  for (i in 2:length(out))
+      print(all.equal(aout[[i]], q2d(out[[i]])))
[1] TRUE
[1] TRUE
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   0.54    0.06    0.59