lattice-bases: Compute a basis with 4ti2

Description

4ti2 provides several executables that can be used to generate bases for a configuration matrix A. See the references for details.

Usage

basis(exec, memoise = TRUE)
zbasis(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
markov(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
groebner(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
hilbert(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
graver(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
fzbasis(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
fmarkov(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
fgroebner(A, format = c("mat", "vec", "tab"), dim = NULL,
  all = FALSE, dir = tempdir(), quiet = TRUE, shell = FALSE,
  dbName = NULL, ...)
fhilbert(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)
fgraver(A, format = c("mat", "vec", "tab"), dim = NULL, all = FALSE,
  dir = tempdir(), quiet = TRUE, shell = FALSE, dbName = NULL, ...)

Value

a matrix containing the Markov basis as its columns (for easy addition to tables)

Arguments

exec: don't use this parameter
memoise: don't use this parameter
A: The configuration matrix
format: How the basis (moves) should be returned. if "mat", the moves are returned as the columns of a matrix.
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
all: If TRUE, all moves (+ and -) are given. if FALSE, only the + moves are given as returned by the executable.
dir: Directory to place the files in, without an ending /
quiet: If FALSE, messages the 4ti2 output
shell: Messages the shell code used to do the computation
dbName: The name of the model in the markov bases database, http://markov-bases.de, see examples
...: Additional arguments to pass to the function, e.g. p = "arb" specifies the flag -parb; not setting this issues a common warning

References

Drton, M., B. Sturmfels, and S. Sullivant (2009). Lectures on Algebraic Statistics, Basel: Birkhauser Verlag AG.

Examples

Run this code



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))



}