# R/markov.r In algstat: Algebraic statistics in R

#### Documented in markov

#' Compute a Markov Basis with 4ti2
#'
#' A Markov basis of a matrix A is computed with the markov function of 4ti2, obtained with the LattE-integrale bundle.
#'
#' @param mat a matrix; for example the output of hmat
#' @param format how the moves should be returned (if "mat", moves are columns)
#' @param dim the dimension to be used in vec2tab if format = "tab" is used, oftentimes 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.
#' @param dir directory to place the files in, without an ending /
#' @param opts options for markov
#' @param quiet show 4ti2 output
#' @param dbName the name of the model in the markov bases database, http://markov-bases.de, see examples
#' @return a matrix containing the Markov basis as its columns (for easy addition to tables)
#' @export markov
#' @references Drton, M., B. Sturmfels, and S. Sullivant (2009). \emph{Lectures on Algebraic Statistics}, Basel: Birkhauser Verlag AG.
#' @examples
#' \dontrun{
#'
#'
#'
#'
#' # 2x2 independence example
#' # following convention, the first index indicates rows
#' varlvls <- c(2,2)
#' facets <- list(1,2)
#' ( A <- hmat(varlvls, facets) )
#' markov(A)
#' markov(A, "vec")
#' markov(A, "tab", varlvls)
#' markov(A, "tab", varlvls, TRUE)
#'
#'
#'
#'
#' # 3x3 independence example
#' # following convention, the first index indicates rows
#' varlvls <- c(3,3)
#' facets <- list(1,2)
#' ( A <- hmat(varlvls, facets) )
#' markov(A)
#' markov(A, "vec")
#' markov(A, "tab", varlvls)
#' markov(A, "tab", varlvls, TRUE)
#'
#'
#'
#'
#' # LAS example 1.2.1, p.12 (2x3 independence)
#' varlvls <- c(2,3)
#' facets <- list(1, 2)
#' ( A <- hmat(varlvls, facets) )
#' markov(A, "tab", varlvls)
#' # Prop 1.2.2 says that there should be
#' 2*choose(2, 2)*choose(3,2) # = 6
#' # moves.
#' markov(A, "tab", varlvls, TRUE)
#'
#'
#'
#'
#'
#' # LAS example 1.2.12, p.17  (no 3-way interaction)
#' varlvls <- c(2,2,2)
#' facets <- list(c(1,2), c(1,3), c(2,3))
#' ( A <- hmat(varlvls, facets) )
#' markov(A)
#  tableau(markov(A), dim = varlvls)
#'
#'
#'
#'
#'
#'
#' # LAS example 1.2.12, p.16  (no 3-way interaction)
#' varlvls <- c(2,2,2,2)
#' facets <- list(c(1,2), c(1,4), c(2,3))
#' ( A <- hmat(varlvls, facets) )
#' markov(A)
#' markov(A, "tab", varlvls) # hard to understand
#' tableau(markov(A), varlvls)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' # using the markov bases database, must be connected to internet
#' # A <- markov(dbName = "ind3-3")
#' B <- markov(hmat(c(3,3), list(1,2)))
#' # all(A == B)
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' markov(diag(1, 10))
#'
#' }
#'
markov <- function(mat, format = c("mat", "vec", "tab"), dim = NULL,
all = FALSE, dir = tempdir(), opts = "-parb", quiet = TRUE,
dbName
){

format <- match.arg(format)

## redirect with the special case of when the identity is given
if((nrow(mat) == ncol(mat)) && all(mat == diag(1, nrow(mat)))){
warning("the identity matrix was supplied to markov, returning 0's.")
if(format == "mat") return(matrix(0, nrow = nrow(mat), ncol = 1))
if(format == "vec") return(list(matrix(0, nrow = nrow(mat), ncol = 1)))
if(format == "tab") return(vec2tab(matrix(0, nrow = nrow(mat), ncol = 1), dim))
}

## make dir to put 4ti2 files in (within the tempdir) timestamped
timeStamp <- as.character(Sys.time())
timeStamp <- chartr("-", "_", timeStamp)
timeStamp <- chartr(" ", "_", timeStamp)
timeStamp <- chartr(":", "_", timeStamp)
dir2 <- file.path2(dir, timeStamp)
suppressWarnings(dir.create(dir2))

## define a function to write the code to a file
formatAndWriteFile <- function(mat, codeFile = "markovCode.mat"){

# line numbers up, e.g. 1 1 0 0
out <- paste(nrow(mat), ncol(mat))
out <- paste0(out, "\n")
out <- paste0(out,
paste(apply(unname(mat), 1, paste, collapse = " "),
collapse = "\n")
)

# write code file
writeLines(out, con = file.path2(dir2, codeFile))
invisible(out)
}

## make 4ti2 file
if(!missing(mat)) formatAndWriteFile(mat)

## switch to temporary directory
oldWd <- getwd()
setwd(dir2)

## create/retrieve markov basis
if(missing(dbName)){

## run 4ti2 if needed
if(.Platform\$OS.type == "unix"){

system2(
file.path2(getOption("markovPath"), "markov"),
paste(opts, file.path2(dir2, "markovCode.mat")),
stdout = "markovOut", stderr = FALSE
)

} else { # windows

matFile <- file.path2(dir2, "markovCode.mat")
matFile <- chartr("\\", "/", matFile)
matFile <- str_c("/cygdrive/c", str_sub(matFile, 3))

system2(
"cmd.exe",
paste(
"/c env.exe",
file.path(getOption("markovPath"), "markov"),
opts, matFile
), stdout = "markovOut", stderr = FALSE
)

}

} else { # if the model name is specified

paste0("http://markov-bases.de/data/", dbName, "/", dbName, ".mar"),
destfile = "markovCode.mat.mar" # already in tempdir
)

}

## figure out what files to keep them, and make 4ti2 object
basis <- basis[-1]
basis <- lapply(
str_split(str_trim(basis), " "),
function(x) as.integer(x[nchar(x) > 0])
)
if(all) basis <- c(basis, lapply(basis, function(x) -x))

## migrate back to original working directory
setwd(oldWd)

# out
if(format == "mat"){
basis <- matrix(unlist(basis), ncol = length(basis[[1]]), byrow = TRUE)
return(t(basis))
} else if(format == "vec") {
return(basis)
} else { # format == "tab"
return(lapply(basis, vec2tab, dim = dim))
}
}


## Try the algstat package in your browser

Any scripts or data that you put into this service are public.

algstat documentation built on May 29, 2017, 10:34 p.m.