R/lin-alg.R

Defines functions matrix_ diag_ vec as_diag `diag<-.caracas_symbol` `diag<-.default` `diag<-` diag.caracas_symbol diag.default diag mat_pow reciprocal_matrix t.caracas_symbol dim.caracas_symbol as_character_matrix number_cols number_rows scalar_to_matrix symbol_is_matrix stopifnot_matrix

Documented in as_character_matrix as_diag diag diag_ diag.caracas_symbol dim.caracas_symbol mat_pow matrix_ reciprocal_matrix t.caracas_symbol vec

stopifnot_matrix <- function(x){
  # MatrixSymbol
  try({
    res <- x$pyobj$is_MatrixExpr
    
    if (is.logical(res) && isTRUE(res)) {
      return(invisible(NULL))
    }
    
    if (reticulate::py_to_r(res)) {
      return(invisible(NULL))
    }
  }, silent = TRUE)
  
  if (!symbol_is_matrix(x)) {
      stop("'x' must be a matrix")
  }
}

symbol_is_matrix <- function(x) {
  # MatrixSymbol
  try({
    res <- x$pyobj$is_MatrixExpr
    
    if (is.logical(res) && isTRUE(res)) {
      return(TRUE)
    }
    
    if (reticulate::py_to_r(res)) {
      return(TRUE)
    }
  }, silent = TRUE)
  
  xstr <- as.character(x)
  
  if (grepl("^Matrix\\(\\[", xstr)) {
    return(TRUE)
  }
  
  # FIXME: From der() and der2()
  # if (grepl("^\\[\\[", xstr)) {
  #   return(TRUE)
  # }
  
  return(FALSE)
}

# Still needed?
#symbol_is_list_of_lists_matrix <- function(x) {
#   if (grepl("^\\[\\[", as.character(x))) {
#     return(TRUE)
#   } 
#   
#   return(FALSE)
# }

# symbol_is_vector <- function(x) {
#   xstr <- as.character(x)
#   
#   if (grepl("^\\[", xstr)) {
#     return(TRUE)
#   }
#   
#   return(FALSE)
# }

scalar_to_matrix <- function(scalar, dims) {
  matrix(rep(scalar, prod(dims)), nrow = dims[1L], ncol = dims[2L])
}


# ensure_symbol_is_matrix <- function(x) {
#   if (!inherits(x, "caracas_symbol")) {
#     stop(paste0("'x' ", TXT_NOT_CARACAS_SYMBOL))
#   }
#   
#   if (!symbol_is_matrix(x)) {
#     stop(paste0("'x' is not a matrix"))
#   }
# }

# On purpose: Does not ensure sympy nor check that x is a matrix
number_rows <- function(x) {
  rows <- x$rows
  
  if (inherits(rows, "python.builtin.int")) {
    rows <- reticulate::py_to_r(rows)
  }
  
  if (inherits(rows, "sympy.core.numbers.Integer")) {
    # FIXME: For MatrixSymbols; avoid as.character()?
    rows <- as.integer(as.character(rows))
  }
  
  return(rows)
}

# On purpose: Does not ensure sympy nor check that x is a matrix
number_cols <- function(x) {
  cols <- x$cols
  
  if (inherits(cols, "python.builtin.int")) {
    cols <- reticulate::py_to_r(cols)
  }
  
  if (inherits(cols, "sympy.core.numbers.Integer")) {
    # FIXME: For MatrixSymbols; avoid as.character()?
    cols <- as.integer(as.character(cols))
  }
  
  return(cols)
}

#' Get matrix as character matrix
#' 
#' @param x caracas symbol
#' 
#' @concept linalg
#' @examples
#' if (has_sympy()) {
#'   s  <- as_sym("[[r1, r2, r3], [u1, u2, u3]]")
#'   s2 <- apply(as_character_matrix(s), 2, function(x) (paste("1/(", x, ")")))
#'   as_sym(s2)
#' }
#' 
#' @export
as_character_matrix <- function(x) {
  y <- as_expr_worker(x, as_character = TRUE)
  return(eval(parse(text = y)))
}


#' Dimensions of a caracas symbol
#' 
#' @param x caracas symbol
#' 
#' @concept linalg
#' 
#' @export
dim.caracas_symbol <- function(x) {
  ensure_sympy()
  stopifnot_symbol(x)
  
  if (!symbol_is_matrix(x)) { 
    return(NULL) 
  }

  rows <- number_rows(x$pyobj)
  cols <- number_cols(x$pyobj)

  return(c(rows, cols))
}



#' Transpose of matrix
#'
#' @param x If `caracas_symbol` treat as such, else
#' call [base::t()].
#'
#' @concept linalg
#' @export
t.caracas_symbol <- function(x) {
  ensure_sympy()
  stopifnot_symbol(x)
  stopifnot_matrix(x)
  
  xT <- x$pyobj$T
  
  return(construct_symbol_from_pyobj(xT))
}


#' Elementwise reciprocal matrix
#'
#' @param x Object `x`
#' @param numerator The numerator in the result.
#'
#' @concept linalg
#' 
#' @examples
#' if (has_sympy()) {
#'   s <- as_sym("[[r1, r2, r3], [u1, u2, u3]]")
#'   reciprocal_matrix(s, numerator = 7)
#' }
#' 
#' @export
reciprocal_matrix <- function(x, numerator = 1){
  ensure_sympy()

  if (!symbol_is_matrix(x)) {
    stop("'x' must be sympy matrix\n")
  }
  
  numerator <- as.character(numerator)
  
  if (grepl("^-?[0-9]+$", numerator)) {
    # S(): Sympify
    # Means that this will be e.g. (S(1))/(1) = 1 instead of 1/1 = 1.0 (numeric)
    numerator <- paste0("S(", numerator, ")")
  }

  rx <- apply(as_character_matrix(x), 2, function(xx) {
    paste0("(", numerator, ")", "/(", xx, ")")
  })
  
  return(as_sym(rx))
}


#' Matrix power
#'
#' @param x A `caracas_symbol`, a matrix.
#' @param pow Power to raise matrix `x` to
#' 
#' @examples
#' if (has_sympy() && sympy_version() >= "1.6") {
#'   M <- matrix_(c("1", "a", "a", 1), 2, 2)
#'   M
#'   mat_pow(M, 1/2)
#' }
#' 
#' @concept linalg
#' @export
mat_pow <- function(x, pow = "1") {
  ensure_sympy()
  stopifnot_symbol(x)
  stopifnot_matrix(x)
  
  #pow_str <- deparse1(substitute(pow)) # to avoid 1/2 gets converted to 0.5
  pow_str <- paste0(deparse(substitute(pow))) # to avoid 1/2 gets converted to 0.5

  if (!("pow" %in% names(x$pyobj))) {
    stop("pow() not supported by SymPy version")
  }
  
  x_pow <- x$pyobj$pow(pow_str)
  
  return(construct_symbol_from_pyobj(x_pow))
}


#' Matrix diagonal
#'
#' @param x Object `x`
#' @param \dots Passed on
#' 
#' @concept linalg
#' 
#' @export
diag <- function(x, ...) {
  UseMethod("diag")
}

#' @export
diag.default <- function(x, ...) {
  return(base::diag(x, ...))
}


#' Matrix diagonal
#' 
#' @param x Object `x`
#' @param \dots Not used
#' 
#' @concept linalg
#' 
#' @export
diag.caracas_symbol <- function(x, ...) {
  ensure_sympy()
  stopifnot_symbol(x)
  stopifnot_matrix(x)  
  
  y <- eval_to_symbol(x$pyobj$diagonal())
  
  return(y)
}


#' Replace matrix diagonal
#'
#' @param x Object `x`
#' @param value Replacement value
#' 
#' @name diag-set
#' 
#' @concept linalg
#' 
#' @export
`diag<-` <- function(x, value) {
  UseMethod("diag<-")
}

#' @export
`diag<-.default` <- function(x, value) {
  return(base::`diag<-`(x, value))
}

#' Replace diagonal
#' 
#' @param x A `caracas_symbol`.
#' @param value Replacement value
#' 
#' @examples 
#' if (has_sympy()) {
#'   A <- matrix(c("a", 0, 0, 0, "a", "a", "a", 0, 0), 3, 3)
#'   B <- as_sym(A)
#'   B
#'   diag(B)
#'   diag(B) <- "b"
#'   B
#'   diag(B)
#' }
#' 
#' @concept vectors
#' 
#' @export
`diag<-.caracas_symbol` <- function(x, value) {
  ensure_sympy()
  stopifnot_symbol(x)
  stopifnot_matrix(x)
  
  xmat <- convert_to_r_mat(x)
  
  if (inherits(value, "caracas_symbol")) {
    value <- c(as_character_matrix(value))
  }
  
  diag(xmat) <- value
  
  y <- as_sym(xmat)
  return(y)
}

#' Construct diagonal matrix from vector
#' 
#' @param x Matrix with 1 row or 1 column that is the 
#' diagonal in a new diagonal matrix
#' 
#' @examples 
#' if (has_sympy()) {
#'   d <- as_sym(c("a", "b", "c"))
#'   D <- as_diag(d)
#'   D
#' }
#' 
#' @concept linalg
#' 
#' @export
as_diag <- function(x) {
  ensure_sympy()
  stopifnot_symbol(x)
  stopifnot_matrix(x)
  
  if (!(nrow(x) == 1L || ncol(x) == 1L)) {
    stop("Either supply a row vector or a column vector with at least one element")
  }
  
  if (nrow(x) > 1L) {
    x <- t(x)
  }
  
  n <- ncol(x)
  
  D <- as_sym(diag(n))
  
  for (i in seq_len(n)) {
    D[i, i] <- as.character(x[1L, i])
  }
  
  return(D)
}


#' Stacks matrix to vector
#' 
#' @param x Matrix
#' 
#' @examples 
#' if (has_sympy()) {
#'   A <- as_sym(matrix(1:9, 3))
#'   vec(A)
#' }
#' 
#' @concept linalg
#' 
#' @export
vec <- function(x) {
  ensure_sympy()
  stopifnot_symbol(x)
  stopifnot_matrix(x)
  
  stopifnot(length(dim(x)) == 2)
  x2 <- do.call(rbind, lapply(seq_len(ncol(x)), function(j) x[, j]))
  
  return(x2)
}



#' Symbolic diagonal matrix
#'
#' @param x Character vector with diagonal
#' @param n Number of times `x` should be repeated
#' @param declare_symbols Passed on to `as_sym()` when constructing symbolic matrix
#' @param \dots Passed on to `rep(x, n, ...)`
#'
#' @concept linalg
#' 
#' @examples
#' if (has_sympy()) {
#'   diag_(c("a", "b", "c"))
#'   diag_("a", 2)
#' }
#' 
#' @export
diag_ <- function(x, n = 1L, declare_symbols = TRUE, ...){
  ensure_sympy()
  
  if (!is.character(x)) {
    stop("'x' must be a character")
  }
  
  x <- rep(x, n, ...)
  n <- length(x)
  
  A <- matrix("0", nrow = n, ncol = n)
  diag(A) <- x
  y <- as_sym(A, declare_symbols = declare_symbols)
  
  return(y)
}

#' Symbolic matrix
#'
#' @param \dots Passed on to [matrix()]
#' @param declare_symbols Passed on to `as_sym()` when constructing symbolic matrix
#'
#' @concept linalg
#' 
#' @examples
#' if (has_sympy()) {
#'   matrix_(1:9, nrow = 3)
#'   matrix_("a", 2, 2)
#' }
#' 
#' @export
matrix_ <- function(..., declare_symbols = TRUE){
  ensure_sympy()
  
  A <- matrix(...)
  y <- as_sym(A, declare_symbols = declare_symbols)
  
  return(y)
}

Try the caracas package in your browser

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

caracas documentation built on Feb. 11, 2022, 9:07 a.m.