# tests/testthat/test-solve.R In caracas: Computer Algebra

```# inv
# solve_lin
# solve_sys

test_that("inv/solve_lin", {
skip_if_no_sympy()

A <- as_sym(matrix(c(1, "x", 0, 2, -1, 3, 4, 2, 5), nrow = 3, ncol = 3))
b <- as_sym(c(4, 1, 7))
x <- A %*% b

# -------------------------
# inv(A)
# -------------------------
Ai <- inv(A)
maybe_b <- simplify(Ai %*% x)
expect_equal(as.character(b), as.character(maybe_b))

# -------------------------
# solve_lin(A)
# -------------------------
Ai2 <- solve_lin(A)
expect_equal(as.character(Ai), as.character(Ai2))

# -------------------------
# solve_lin(A, x)
# -------------------------
maybe_b <- simplify(solve_lin(A, x))
expect_equal(as.character(b), as.character(maybe_b))
})

test_that("solve_sys(lhs, vars)", {
skip_if_no_sympy()

###################################################################
# Single variable
###################################################################
#------------------------------------------------------------------
sol1 <- solve_sys(as_sym("x**2 + 1"), "x")

expect_equal(length(sol1), 2L)
expect_equal(unname(unlist(lapply(sol1, lapply, as.character))),
c("-1i", "1i"))
expect_equal(names(sol1[[1L]]), "x")
expect_equal(as.character(sol1[[1L]]\$x), "-1i")
expect_equal(names(sol1[[2L]]), "x")
expect_equal(as.character(sol1[[2L]]\$x), "1i")
#------------------------------------------------------------------

#------------------------------------------------------------------
x <- symbol("x")
sol2 <- solve_sys(as_sym("x**2 + 1"), x)

expect_equal(as.character(sol1), as.character(sol2))

expect_equal(length(sol2), 2L)
expect_equal(unname(unlist(lapply(sol2, lapply, as.character))),
c("-1i", "1i"))
expect_equal(names(sol2[[1L]]), "x")
expect_equal(as.character(sol2[[1L]]\$x), "-1i")
expect_equal(names(sol2[[2L]]), "x")
expect_equal(as.character(sol2[[2L]]\$x), "1i")
#------------------------------------------------------------------

###################################################################
# Multiple variables
###################################################################
# Must be as a row vector (1 x m matrix), not a list....
lhs <- t(as_sym(matrix(c("x**2 + 1", "y+3"))))
sol1 <- solve_sys(lhs, c("x", "y"))
sol1_ord <- order(unlist(lapply(sol1, function(l) Arg(as_expr(l\$x)))))
sol1 <- sol1[sol1_ord]

expect_equal(length(sol1), 2L)
expect_equal(sort(names(sol1[[1L]])), c("x", "y"))
expect_equal(as.character(sol1[[1L]]\$x), "-1i")
expect_equal(as.character(sol1[[1L]]\$y), "-3")
expect_equal(sort(names(sol1[[2L]])), c("x", "y"))
expect_equal(as.character(sol1[[2L]]\$x), "1i")
expect_equal(as.character(sol1[[2L]]\$y), "-3")

y <- symbol("y")
sol2 <- solve_sys(lhs, list(x, y))
sol2_ord <- order(unlist(lapply(sol2, function(l) Arg(as_expr(l\$x)))))
sol2 <- sol2[sol2_ord]

expect_equal(length(sol2), 2L)
expect_equal(sort(names(sol2[[1L]])), c("x", "y"))
expect_equal(as.character(sol2[[1L]]\$x), "-1i")
expect_equal(as.character(sol2[[1L]]\$y), "-3")
expect_equal(sort(names(sol2[[2L]])), c("x", "y"))
expect_equal(as.character(sol2[[2L]]\$x), "1i")
expect_equal(as.character(sol2[[2L]]\$y), "-3")
})

test_that("solve_sys(lhs, rhs, vars)", {
skip_if_no_sympy()

sol1 <- solve_sys(as_sym("x**2"), as_sym("-1"), "x")
#sol1

expect_equal(length(sol1), 2L)
expect_equal(unname(unlist(lapply(sol1, lapply, as.character))),
c("-1i", "1i"))
expect_equal(names(sol1[[1L]]), "x")
expect_equal(as.character(sol1[[1L]]\$x), "-1i")
expect_equal(names(sol1[[2L]]), "x")
expect_equal(as.character(sol1[[2L]]\$x), "1i")
})

test_that("solve system of non-linear equations", {
skip_if_no_sympy()

# Multinomial likelihood
p <- as_sym(paste0("p", 1:3))
y <- as_sym(paste0("y", 1:3))
a <- as_sym("a")
l <- sum(y*log(p))
L <- -l + a*(sum(p) - 1)
g <- der(L, list(a, p))
expect_match(as.character(g),
"[p1 + p2 + p3 - 1, a - y1/p1, a - y2/p2, a - y3/p3]",
fixed = TRUE)

sols <- solve_sys(g, list(a, p))
expect_equal(length(sols), 1L)

sol <- sols[[1L]]
expect_equal(sort(names(sol)), sort(c("p1", "p2", "p3", "a")), fixed = TRUE)
expect_match(as.character(sol\$p1), "y1/(y1 + y2 + y3)", fixed = TRUE)
expect_match(as.character(sol\$p2), "y2/(y1 + y2 + y3)", fixed = TRUE)
expect_match(as.character(sol\$p3), "y3/(y1 + y2 + y3)", fixed = TRUE)
expect_match(as.character(sol\$a), "y1 + y2 + y3", fixed = TRUE)

H <- der2(L, list(p, a))
H_sol <- subs_lst(H, sol)
expect_match(as.character(H_sol),
"[[(y1 + y2 + y3)^2/y1, 0, 0, 1], [0, (y1 + y2 + y3)^2/y2, 0, 1], [0, 0, (y1 + y2 + y3)^2/y3, 1], [1, 1, 1, 0]]",
fixed = TRUE)
})

test_that("input variations", {
skip_if_no_sympy()

# Analysis of variance - one-way ANOVA
ngrp <- 3   # number of groups
spg  <- 2   # number of subjects per group
g <- seq_len(ngrp)
f <- factor(rep(g, each = spg))
y <- as_sym(matrix(paste0("y", seq_along(f))))
X <- as_sym(model.matrix(~ f))
XtX <- t(X) %*% X
XtXinv <- inv(XtX)
Xty <- t(X) %*% y
beta_hat <- XtXinv %*% Xty

beta <- as_sym(paste0("beta", 1:3))
res <- y - X %*% beta
logL <- - RSS / 2

Score <- der(logL, beta) %>% matrify()

beta_lst <- listify(beta)
sol1 <- solve_sys(Score, beta_lst)
sol2 <- solve_sys(Score, beta)

expect_equal(length(sol1), 1L)
expect_equal(length(sol2), 1L)

sol1 <- sol1[[1L]]
sol2 <- sol2[[1L]]

expect_equal(paste0(sapply(sol1, as.character), collapse = ";"),
"y1/2 + y2/2;-y1/2 - y2/2 + y3/2 + y4/2;-y1/2 - y2/2 + y5/2 + y6/2")

expect_equal(paste0(sapply(sol1, as.character), collapse = ";"),
paste0(sapply(sol2, as.character), collapse = ";"))

sol1_mat <- t(do.call(cbind, sol1))
expect_equal(as.character(sol1_mat),
"Matrix([[y1/2 + y2/2], [-y1/2 - y2/2 + y3/2 + y4/2], [-y1/2 - y2/2 + y5/2 + y6/2]])")
expect_equal(as.character(sol1_mat),
as.character(beta_hat))

})
```

## 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.