#' Compute a basis with 4ti2
#'
#' 4ti2 provides several executables that can be used to generate bases for a
#' configuration matrix A. See the references for details.
#'
#' @param A The configuration matrix
#' @param format How the basis (moves) should be returned. if "mat", the moves
#' are returned as the columns of a matrix.
#' @param dim The dimension to be passed to [vec2tab()] if format = "tab" is
#' used; a vector of the number of levels of each variable in order
#' @param all If TRUE, all moves (+ and -) are given. if FALSE, only the +
#' moves are given as returned by the executable.
#' @param dir Directory to place the files in, without an ending /
#' @param quiet If FALSE, messages the 4ti2 output
#' @param shell Messages the shell code used to do the computation
#' @param dbName The name of the model in the markov bases database,
#' http://markov-bases.de, see examples
#' @param ... Additional arguments to pass to the function, e.g. \code{p =
#' "arb"} specifies the flag \code{-parb}; not setting this issues a common
#' warning
#' @return a matrix containing the Markov basis as its columns (for easy
#' addition to tables)
#' @name lattice-bases
#' @references Drton, M., B. Sturmfels, and S. Sullivant (2009). \emph{Lectures
#' on Algebraic Statistics}, Basel: Birkhauser Verlag AG.
#' @examples
#'
#'
#' if (has_4ti2()) {
#'
#'
#' # basic input and output for the 3x3 independence example
#' (A <- rbind(
#' kprod(diag(3), ones(1,3)),
#' kprod(ones(1,3), diag(3))
#' ))
#' markov(A, p = "arb")
#'
#'
#'
#' # you can get the output formatted in different ways:
#' markov(A, p = "arb", all = TRUE)
#' markov(A, p = "arb", "vec")
#' markov(A, p = "arb", "tab", c(3, 3))
#' tableau(markov(A, p = "arb"), dim = c(3, 3)) # tableau notation
#'
#'
#'
#' # you can add options by listing them off
#' # to see the options available to you by function,
#' # go to http://www.4ti2.de
#' markov(A, p = "arb")
#'
#'
#'
#' # the basis functions are automatically cached for future use.
#' # (note that it doesn't persist across sessions.)
#' A <- rbind(
#' kprod( diag(4), ones(1,4), ones(1,4)),
#' kprod(ones(1,4), diag(4), ones(1,4)),
#' kprod(ones(1,4), ones(1,4), diag(4))
#' )
#' system.time(markov(A, p = "arb"))
#' system.time(markov(A, p = "arb"))
#'
#' # the un-cashed versions begin with an "f"
#' # (think: "forgetful" markov)
#' system.time(fmarkov(A, p = "arb"))
#' system.time(fmarkov(A, p = "arb"))
#'
#'
#'
#' # you can see the command line code by typing shell = TRUE
#' # and the standard output wiht quiet = FALSE
#' # we illustrate these with fmarkov because otherwise it's cached
#' (A <- rbind(
#' kprod(diag(2), ones(1,4)),
#' kprod(ones(1,4), diag(2))
#' ))
#' fmarkov(A, p = "arb", shell = TRUE)
#' fmarkov(A, p = "arb", quiet = FALSE)
#'
#'
#'
#' # compare the bases for the 3x3x3 no-three-way interaction model
#' A <- rbind(
#' kprod( diag(3), diag(3), ones(1,3)),
#' kprod( diag(3), ones(1,3), diag(3)),
#' kprod(ones(1,3), diag(3), diag(3))
#' )
#' str( zbasis(A, p = "arb")) # 8 elements = ncol(A) - qr(A)$rank
#' str( markov(A, p = "arb")) # 81 elements
#' str(groebner(A, p = "arb")) # 110 elements
#' str( graver(A)) # 795 elements
#'
#'
#' # the other bases are also cached
#' A <- rbind(
#' kprod( diag(3), ones(1,3), ones(1,2)),
#' kprod(ones(1,3), diag(3), ones(1,2)),
#' kprod(ones(1,3), ones(1,3), diag(2))
#' )
#' system.time( graver(A))
#' system.time( graver(A))
#' system.time(fgraver(A))
#' system.time(fgraver(A))
#'
#'
#'
#' # LAS ex 1.2.1, p.12 : 2x3 independence
#' (A <- rbind(
#' kprod(diag(2), ones(1,3)),
#' kprod(ones(1,2), diag(3))
#' ))
#'
#' markov(A, p = "arb", "tab", c(3, 3))
#' # Prop 1.2.2 says that there should be
#' 2*choose(2, 2)*choose(3,2) # = 6
#' # moves (up to +-1)
#' markov(A, p = "arb", "tab", c(3, 3), TRUE)
#'
#'
#'
#' # LAS example 1.2.12, p.17 (no 3-way interaction)
#' (A <- rbind(
#' kprod( diag(2), diag(2), ones(1,2)),
#' kprod( diag(2), ones(1,2), diag(2)),
#' kprod(ones(1,2), diag(2), diag(2))
#' ))
#' plot_matrix(A)
#' markov(A, p = "arb")
#' groebner(A, p = "arb")
#' graver(A)
#' tableau(markov(A, p = "arb"), dim = c(2,2,2))
#'
#'
#'
#'
#'
#' # using the markov bases database, must be connected to internet
#' # commented out for predictable and fast cran checks time
#' # A <- markov(dbName = "ind3-3")
#' # B <- markov(rbind(
#' # kprod(diag(3), ones(1,3)),
#' # kprod(ones(1,3), diag(3))
#' # ), p = "arb")
#' # all(A == B)
#'
#'
#'
#'
#'
#' # possible issues
#' # markov(diag(1, 10))
#' # zbasis(diag(1, 10), "vec")
#' # groebner(diag(1, 10), "vec", all = TRUE)
#' # graver(diag(1, 10), "vec", all = TRUE)
#' # graver(diag(1, 4), "tab", all = TRUE, dim = c(2,2))
#'
#'
#'
#' }
#'
#'
#' @param exec don't use this parameter
#' @param memoise don't use this parameter
basis <- function(exec, memoise = TRUE){
## stuff in basis
extension <- switch(exec,
markov = ".mar",
groebner = ".gro",
hilbert = ".hil",
graver = ".gra",
zbasis = ".lat",
zsolve = ".zfree"
)
commonName <- switch(exec,
markov = "markov",
groebner = "grobner",
hilbert = "hilbert",
graver = "graver",
zbasis = "lattice",
zsolve = "solve"
)
## memoise or not
mem_or_not <- if(memoise) memoise::memoise else (function(x) x)
## create the function to return
mem_or_not(function(A, format = c("mat", "vec", "tab"), dim = NULL,
all = FALSE, dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL,
...
){
if (!has_4ti2()) missing_4ti2_stop()
## check args
format <- match.arg(format)
if (format == "tab" && missing(dim)) {
stop('If format = "tab" is specified, dim must be also.', call. = FALSE)
}
## compute other args
opts <- as.list(match.call(expand.dots = FALSE))[["..."]]
if (is.null(opts)) {
opts <- ""
} else {
opts <- str_c("-", names(opts), "", unlist(opts))
opts <- str_c(opts, collapse = " ")
}
## make dir to put 4ti2 files in (within the tempdir) timestamped
dir.create(scratch_dir <- file.path(dir, time_stamp()))
## make 4ti2 file
if(!missing(A)) write.latte(A, file.path(scratch_dir, "PROJECT.mat"))
## switch to temporary directory
user_working_directory <- getwd()
setwd(scratch_dir); on.exit(setwd(user_working_directory), add = TRUE)
## create/retrieve markov basis
if (is.null(dbName)) {
## run 4ti2 if needed
if (is_mac() || is_unix()) {
system2(
file.path(get_4ti2_path(), exec),
paste(opts, file.path(scratch_dir, "PROJECT")),
stdout = glue("{exec}_out"),
stderr = glue("{exec}_err")
)
# generate shell code
shell_code <- glue(
"{file.path(get_4ti2_path(), exec)} {paste(opts, file.path(scratch_dir, 'PROJECT'))} > {exec}_out 2> {exec}_err"
)
if(shell) message(shell_code)
} else if (is_win()) {
matFile <- file.path(scratch_dir, "PROJECT")
matFile <- chartr("\\", "/", matFile)
matFile <- str_c("/cygdrive/c", str_sub(matFile, 3))
system2(
"cmd.exe",
glue("/c env.exe {file.path(get_4ti2_path(), exec)} {opts} {matFile}"),
stdout = glue("{exec}_out"),
stderr = glue("{exec}_err")
)
# generate shell code
shell_code <- glue(
"cmd.exe /c env.exe {file.path(get_4ti2_path(), exec)} {opts} {matFile} > {exec}_out 2> {exec}_err"
)
if(shell) message(shell_code)
}
if(!quiet) cat(readLines(glue("{exec}_out")), sep = "\n")
std_err <- readLines(glue("{exec}_err"))
if(any(std_err != "")) warning(str_c(std_err, collapse = "\n"), call. = FALSE)
} else { # if the model name is specified
download.file(
paste0("http://markov-bases.de/data/", dbName, "/", dbName, extension),
destfile = "PROJECT.mar" # already in tempdir
)
}
## fix case of no graver basis
if(exec == "graver"){
if(paste0("PROJECT", extension) %notin% list.files(scratch_dir)){
warning(sprintf("%s basis empty, returning 0's.", capitalize(commonName)), call. = FALSE)
return(fix_graver(A, format, dim))
}
}
## figure out what files to keep them, and make 4ti2 object
basis <- t(read.latte(paste0("PROJECT", extension)))
## fix case of no basis
basisDim <- dim(basis)
noBasisFlag <- FALSE
if(any(basisDim == 0)){
noBasisFlag <- TRUE
warning(sprintf("%s basis empty, returning 0's.", capitalize(commonName)), call. = FALSE)
basisDim[basisDim == 0] <- 1L
basis <- rep(0L, prod(basisDim))
dim(basis) <- basisDim
}
## format
if(all && !noBasisFlag) basis <- cbind(basis, -basis)
# out
if(format == "mat"){
return(basis)
} else {
lbasis <- as.list(rep(NA, ncol(basis)))
for(k in 1:ncol(basis)) lbasis[[k]] <- basis[,k]
if(format == "vec") return(lbasis)
if(format == "tab") return(lapply(lbasis, vec2tab, dim = dim))
}
})
}
# #' @export
# #' @rdname lattice-bases
# zsolve <- basis("zsolve", memoise = TRUE)
#' @export
#' @rdname lattice-bases
zbasis <- basis("zbasis", memoise = TRUE)
#' @export
#' @rdname lattice-bases
markov <- basis("markov", memoise = TRUE)
#' @export
#' @rdname lattice-bases
groebner <- basis("groebner", memoise = TRUE)
#' @export
#' @rdname lattice-bases
hilbert <- basis("hilbert", memoise = TRUE)
#' @export
#' @rdname lattice-bases
graver <- basis("graver", memoise = TRUE)
# #' @export
# #' @rdname lattice-bases
# fzsolve <- basis("zsolve", memoise = FALSE)
#' @export
#' @rdname lattice-bases
fzbasis <- basis("zbasis", memoise = FALSE)
#' @export
#' @rdname lattice-bases
fmarkov <- basis("markov", memoise = FALSE)
#' @export
#' @rdname lattice-bases
fgroebner <- basis("groebner", memoise = FALSE)
#' @export
#' @rdname lattice-bases
fhilbert <- basis("hilbert", memoise = FALSE)
#' @export
#' @rdname lattice-bases
fgraver <- basis("graver", memoise = FALSE)
fix_graver <- function(A, format, dim){
if(format == "mat"){
return(matrix(0L, nrow = ncol(A)))
} else {
lbasis <- list(rep(0L, ncol(A)))
if(format == "vec") return(lbasis)
if(format == "tab") return(lapply(lbasis, vec2tab, dim = dim))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.