# R/calculus.R In caracas: Computer Algebra

#### Documented in derder2drop_remainderhessianintjacobianlimprod_scoresum_taylor

```bound_to_str <- function(b) {
if (is.character(b)) {
return(b)
}

if (inherits(b, "caracas_symbol")) {
return(as.character(b))
}

if (is.infinite(b)) {
if (b > 0) {
return("oo")
} else {
return("-oo")
}
}

bnd <- as.character(b)
#bnd <- deparse(eval(substitute(substitute(b)), parent.frame()))
#bnd <- gsub("pi", "Pi", bnd, fixed = TRUE)

return(bnd)
}

calc_verify_func <- function(f) {
if (!inherits(f, "caracas_symbol")) {
stop(paste0("'f' ", TXT_NOT_CARACAS_SYMBOL))
}
}

#' Limit of a function
#'
#' @param f Function to take limit of
#' @param var Variable to take limit for (either string or `caracas_symbol`)
#' @param val Value for `var` to approach
#' @param dir Direction from where `var` should approach `val`: `'+'` or `'-'`
#' @param doit Evaluate the limit immediately (or later with [doit()])
#'
#' @examples
#' if (has_sympy()) {
#'   x <- symbol("x")
#'   lim(sin(x)/x, "x", 0)
#'   lim(1/x, "x", 0, dir = '+')
#'   lim(1/x, "x", 0, dir = '-')
#' }
#'
#' @concept calculus
#'
#' @export
lim <- function(f, var, val, dir = NULL, doit = TRUE) {
calc_verify_func(f)
ensure_sympy()
var <- as.character(var)
verify_variable_name(var)

val_str <- bound_to_str(val)

y <- if (!is.null(dir) && length(dir) == 1L && dir %in% c('+', '-')) {
if (doit) {
get_sympy()\$limit(f\$pyobj, var, val_str, dir)
} else {
get_sympy()\$Limit(f\$pyobj, var, val_str, dir)
}
} else {
if (doit) {
get_sympy()\$limit(f\$pyobj, var, val_str)
} else {
get_sympy()\$Limit(f\$pyobj, var, val_str)
}
}

z <- construct_symbol_from_pyobj(y)

return(z)
}

#' Sum of a function
#'
#' @param f Function to take sum of
#' @param var Variable to take sum for (either string or `caracas_symbol`)
#' @param lower Lower limit
#' @param upper Upper limit
#' @param doit Evaluate the sum immediately (or later with [doit()])
#'
#' @examples
#' if (has_sympy()) {
#'   x <- symbol("x")
#'   s <- sum_(1/x, "x", 1, 10)
#'   as_expr(s)
#'   sum(1/(1:10))
#'   n <- symbol("n")
#'   simplify(sum_(x, x, 1, n))
#' }
#'
#' @concept calculus
#'
#' @export
sum_ <- function(f, var, lower, upper, doit = TRUE) {
calc_verify_func(f)
ensure_sympy()
var <- as.character(var)
verify_variable_name(var)

lwr_str <- bound_to_str(lower)
upr_str <- bound_to_str(upper)

y <- if (doit) {
get_sympy()\$summation(f\$pyobj, c(var, lwr_str, upr_str))
} else {
get_sympy()\$Sum(f\$pyobj, c(var, lwr_str, upr_str))
}

z <- construct_symbol_from_pyobj(y)

return(z)
}

#' Product of a function
#'
#' @param f Function to take product of
#' @param var Variable to take product for (either string or `caracas_symbol`)
#' @param lower Lower limit
#' @param upper Upper limit
#' @param doit Evaluate the product immediately (or later with [doit()])
#'
#' @examples
#' if (has_sympy()) {
#'   x <- symbol("x")
#'   p <- prod_(1/x, "x", 1, 10)
#'   p
#'   as_expr(p)
#'   prod(1/(1:10))
#'   n <- symbol("n")
#'   prod_(x, x, 1, n)
#' }
#'
#' @concept calculus
#'
#' @export
prod_ <- function(f, var, lower, upper, doit = TRUE) {
calc_verify_func(f)
ensure_sympy()
var <- as.character(var)
verify_variable_name(var)

lwr_str <- bound_to_str(lower)
upr_str <- bound_to_str(upper)

y <- if (doit) {
get_sympy()\$product(f\$pyobj, c(var, lwr_str, upr_str))
} else {
get_sympy()\$Product(f\$pyobj, c(var, lwr_str, upr_str))
}

z <- construct_symbol_from_pyobj(y)

return(z)
}

#' Integrate a function
#'
#' If no limits are provided, the
#' indefinite integral is calculated.
#' Otherwise, if both limits are provided,
#' the definite integral is calculated.
#'
#' @param f Function to integrate
#' @param var Variable to integrate with respect to (either string or `caracas_symbol`)
#' @param lower Lower limit
#' @param upper Upper limit
#' @param doit Evaluate the integral immediately (or later with [doit()])
#'
#' @examples
#' if (has_sympy()) {
#'   x <- symbol("x")
#'
#'   int(1/x, x, 1, 10)
#'   int(1/x, x, 1, 10, doit = FALSE)
#'   int(1/x, x)
#'   int(1/x, x, doit = FALSE)
#'   int(exp(-x^2/2), x, -Inf, Inf)
#'   int(exp(-x^2/2), x, -Inf, Inf, doit = FALSE)
#' }
#'
#' @concept calculus
#'
#' @export
int <- function(f, var, lower, upper, doit = TRUE) {
calc_verify_func(f)
ensure_sympy()
var <- as.character(var)
verify_variable_name(var)

with_limits <- FALSE

if (missing(lower) && missing(upper)) {
} else if (!missing(lower) && !missing(upper)) {
with_limits <- TRUE
} else {
stop("Either both lower and upper must be given or neither must be given")
}

if (with_limits) {
lwr_str <- bound_to_str(lower)
upr_str <- bound_to_str(upper)

y <- if (doit) {
get_sympy()\$integrate(f\$pyobj, c(var, lwr_str, upr_str))
} else {
get_sympy()\$Integral(f\$pyobj, c(var, lwr_str, upr_str))
}

z <- construct_symbol_from_pyobj(y)

return(z)
}

# No limits:

# Potentially replacing existing symbol:
var_symb <- reticulate::py_eval(paste0("symbols('", var, "')"), convert = FALSE)

y <- if (doit) {
get_sympy()\$integrate(f\$pyobj, var_symb)
} else {
get_sympy()\$Integral(f\$pyobj, var_symb)
}

z <- construct_symbol_from_pyobj(y)

return(z)
}

# expr: caracas symbol
# vars: array of symbols (potentially just 1)
der_worker <- function(expr, vars) {
# vars <- vars_to_array(vars)
ensure_sympy()
stopifnot(inherits(vars, "caracas_symbol"))

multi_var <- grepl("^\\[", as.character(vars))

if (!multi_var) {
z <- get_sympy()\$diff(expr\$pyobj, vars\$pyobj)
v <- construct_symbol_from_pyobj(z)
return(v)
}

z <- get_sympy()\$derive_by_array(expr\$pyobj, vars\$pyobj)
v <- construct_symbol_from_pyobj(z)
return(v)
}

# Convert vars to array; vars can be many different things...
vars_to_array <- function(vars) {
# vars should already be in expr, so
# they are known as symbols on the Python side

vars_chr <- if (is.character(vars)) {
# character vector
vars
} else if (inherits(vars, "caracas_symbol")) {
# A single caracas_symbol
as.vector(as_character_matrix(vars))
} else if (inherits(vars, "list")) {
# A list of caracas_symbols
for (i in seq_along(vars)) {
v <- vars[[i]]

if (!inherits(v, "caracas_symbol")) {
stop("Element in vars had wrong type (expected caracas_symbol)")
}
}

unlist(lapply(vars, function(l) as.vector(as_character_matrix(l))))
} else {
stop("Unexpected vars type")
}

if (length(vars_chr) == 1L) {
v <- eval_to_symbol(vars_chr)
return(v)
}

# > 1 variable:
v <- eval_to_symbol(paste0("[", paste0(vars_chr, collapse = ", "), "]"))
return(v)
}

#' Symbolic differentiation of an expression
#'
#' @param expr A `caracas_symbol`
#' @param vars variables to take derivate with respect to
#' @param simplify Simplify result
#'
#' @examples
#' if (has_sympy()) {
#'   x <- symbol("x")
#'   y <- symbol("y")
#'   f <- 3*x^2 + x*y^2
#'   der(f, x)
#'   g <- der(f, list(x, y))
#'   g
#'   dim(g)
#'   G <- matrify(g)
#'   G
#'   dim(G)
#'
#'   h <- der(g, list(x, y))
#'   h
#'   dim(h)
#'   as.character(h)
#'   H <- matrify(h)
#'   H
#'   dim(H)
#'
#'   g %>%
#'     der(list(x, y), simplify = FALSE) %>%
#'     der(list(x, y), simplify = FALSE) %>%
#'     der(list(x, y), simplify = FALSE)
#' }
#'
#' @concept calculus
#'
#' @export
der <- function(expr, vars, simplify = TRUE) {
ensure_sympy()

new_vars <- vars_to_array(vars)
d <- der_worker(expr, new_vars)

if (simplify) {
if (grepl("^\\[\\[\\[", as.character(d))) {
d <- unbracket(d)
d <- matrify(d)
}
}

return(d)
}

#' Symbolic differentiation of second order of an expression
#'
#' @param expr A `caracas_symbol`
#' @param vars variables to take derivate with respect to
#' @param simplify Simplify result
#'
#' @examples
#' if (has_sympy()) {
#'   x <- symbol("x")
#'   y <- symbol("y")
#'   f <- 3*x^2 + x*y^2
#'   der2(f, x)
#'   h <- der2(f, list(x, y))
#'   h
#'   dim(h)
#'   H <- matrify(h)
#'   H
#'   dim(H)
#' }
#'
#' @concept calculus
#'
#' @export
der2 <- function(expr, vars, simplify = TRUE) {
ensure_sympy()

d1 <- der(expr, vars, simplify = simplify)
d2 <- der(d1, vars, simplify = simplify)

return(d2)
}

#' Remove remainder term
#'
#' @param x Expression to remove remainder term from
#'
#' @examples
#' if (has_sympy()) {
#'   def_sym(x)
#'   f <- cos(x)
#'   ft_with_O <- taylor(f, x0 = 0, n = 4+1)
#'   ft_with_O
#'   ft_with_O %>% drop_remainder() %>% as_expr()
#' }
#'
#' @seealso [taylor()]
#'
#' @concept calculus
#'
#' @export
drop_remainder <- function(x) {
ensure_sympy()

ft <- sympy_func(x = x, fun = "removeO")

return(ft)
}

#' Taylor expansion
#'
#' @param f Function to be expanded
#' @param x0 Point to expand around
#' @param n Order of remainder term
#'
#' @examples
#' if (has_sympy()) {
#'   def_sym(x)
#'   f <- cos(x)
#'   ft_with_O <- taylor(f, x0 = 0, n = 4+1)
#'   ft_with_O
#'   ft_with_O %>% drop_remainder() %>% as_expr()
#' }
#'
#' @seealso [drop_remainder()]
#'
#' @concept calculus
#'
#' @export
taylor <- function(f, x0 = 0, n = 6) {
ensure_sympy()

ft <- f %>% sympy_func("series", x0 = x0, n = n)

return(ft)
}

#' Score and Hessian matrix
#'
#' Compute column vector of first derivatives and matrix of second
#' derivatives of univariate function.
#'
#' @name score_hessian
#' @param expr 'caracas expression'.
#' @param vars variables to take derivative with respect to.
#' @param simplify Try to simplify result using `simplify()`; may be time consuming.
#' @seealso [jacobian()], [der()]
#' @examples
#'
#' if (has_sympy()) {
#'   def_sym(b0, b1, x, x0)
#'   f <- b0 / (1 + exp(b1*(x-x0)))
#'   S <- score(f, c(b0, b1))
#'   S
#'   H <- hessian(f, c(b0, b1))
#'   H
#' }
#'
#' @concept calculus
#' @export
#' @rdname score_hessian
score <- function(expr, vars, simplify=TRUE){
ensure_sympy()
stopifnot_symbol(expr)
if (!symbol_is_atomic(expr)){
stop("'expr' must be atomic")
}
out <- der(expr, vars, simplify = simplify)
out <- to_matrix(out)
return(out)
}

#' @export
#' @rdname score_hessian
hessian <- function(expr, vars, simplify=TRUE){
ensure_sympy()
stopifnot_symbol(expr)
if (!symbol_is_atomic(expr)){
stop("'expr' must be atomic")
}
out <- der2(expr, vars, simplify = simplify)
return(to_matrix(out))
}

#' Compute Jacobian
#'
#' @param expr 'caracas expression'.
#' @param vars variables to take derivative with respect to
#'
#' @seealso [score()], [hessian()] [der()]
#' @examples
#' if (has_sympy()) {
#'   x <- paste0("x", seq_len(3))
#'   def_sym_vec(x)
#'   y1 <- x1 + x2
#'   y2 <- x1^2 + x3
#'   y <- c(y1, y2)
#'   jacobian(y, x)
#'   u <- 2 + 4*x1^2
#'   jacobian(u, x1)
#' }
#' @concept calculus
#' @export
jacobian <- function(expr, vars) {
ensure_sympy()
stopifnot_symbol(expr)
out <- der(expr, vars)

if (is.null(dim(out))) {
out <- matrify(out)
}

out <- t(out)
return(out)
}
```

## Try the caracas package in your browser

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

caracas documentation built on Oct. 17, 2023, 5:08 p.m.