m2r: Interface to 'Macaulay2'

#' Smith normal form
#'
#' For an integer matrix M, this computes the matrices D, P, and Q such that
#' \emph{D = PMQ}, which can be seen as an analogue of the singular value
#' decomposition. All are integer matrices, and P and Q are unimodular (have
#' determinants +- 1).
#'
#' @param mat a matrix (integer entries)
#' @param code return only the M2 code? (default: \code{FALSE})
#' @name snf
#' @seealso [m2_matrix()]
#' @return a list of \code{m2_matrix} objects with names \code{D}, \code{P}, and
#'   \code{Q}
#' @examples
#'
#' \dontrun{ requires Macaulay2
#'
#' ##### basic usage
#' ########################################
#'
#' M <- matrix(c(
#'    2,  4,   4,
#'   -6,  6,  12,
#'   10, -4, -16
#' ), nrow = 3, byrow = TRUE)
#'
#' snf(M)
#'
#' (mats <- snf(M))
#' P <- mats$P; D <- mats$D; Q <- mats$Q
#'
#' P %*% M %*% Q                # = D
#' solve(P) %*% D %*% solve(Q)  # = M
#'
#' det(P)
#' det(Q)
#'
#'
#' M <- matrix(c(
#'      1,    2,    3,
#'      1,   34,   45,
#'   2213, 1123, 6543,
#'      0,    0,    0
#' ), nrow = 4, byrow = TRUE)
#' (mats <- snf(M))
#' P <- mats$P; D <- mats$D; Q <- mats$Q
#' P %*% M %*% Q                # = D
#'
#'
#'
#' ##### understanding lattices
#' ########################################
#'
#'
#'
#'
#' # cols of m generate the lattice L
#' M <- matrix(c(2,-1,1,3), nrow = 2)
#' row.names(M) <- c("x", "y")
#' M
#'
#' # plot lattice
#' df <- expand.grid(x = -20:20, y = -20:20)
#' pts <- t(apply(df, 1, function(v) M %*% v))
#' w <- c(-15, 15)
#' plot(pts, xlim = w, ylim = w)
#'
#' # decompose m
#' (mats <- snf(M))
#' P <- mats$P; D <- mats$D; Q <- mats$Q
#'
#' # PMQ = D, the columns of MQ = P^(-1) D  are a simpler basis of
#' # the lattice generated by (the cols of) M
#' (basis <- solve(P) %*% D)
#'
#' # plot lattice generated by new basis
#' pts2 <- t(apply(df, 1, function(v) basis %*% v))
#' points(pts2, pch = "*", col = "red")
#'
#'
#'
#' ##### other options
#' ########################################
#'
#' snf.(M)
#' snf(M, code = TRUE)
#'
#'
#'
#'
#' }
#'






#' @rdname snf
#' @export
snf <- function (mat, code = FALSE) {

  # run m2
  args <- as.list(match.call())[-1]
  eargs <- lapply(args, eval, envir = parent.frame())
  pointer <- do.call(snf., eargs)
  if(code) return(invisible(pointer))

  # parse output
  parsed_out <- m2_parse(pointer)

  # list and out
  list(
    D = parsed_out[[1]],
    P = parsed_out[[2]],
    Q = parsed_out[[3]]
  )
}





#' @rdname snf
#' @export
snf. <- function (mat, code = FALSE) {

  # arg checking
  # if (is.m2_matrix(mat)) mat <- mat$rmatrix
  if (is.m2_pointer(mat)) {
    param <- m2_name(mat)
  } else {
    if (!is.integer(mat)) stopifnot(all(mat == as.integer(mat)))
    param <- paste0("matrix", listify_mat(mat))
  }

  # create code and message
  m2_code <- sprintf("smithNormalForm %s", param)
  if(code) { message(m2_code); return(invisible(m2_code)) }

  # run m2 and return pointer
  m2.(m2_code)
}