Computer algebra systems in R:

# Computer algebra systems in <code>R:</code>
## <code>caracas</code> and <code>Ryacas</code>
### Mikkel Meyer Andersen and Søren Højsgaard
### e-RUM 2020

---

---

Ask a question in Hopin session chat:

`@Mikkel Meyer Andersen #Q2S Question...`

---

# Background

---

#### `Ryacas`

* Based on 'Yet Another Computer Algebra system' ([Yacas](http://www.yacas.org/))

--
* `Ryacas` on CRAN around 2007
  + Authors: Rob Goedman, Gabor Grothendieck, **Søren Højsgaard**, Ayal Pinkus, Grzegorz Mazur
  + 2017: **Mikkel Meyer Andersen** became maintainer 
  + [JOSS paper *Andersen and Højsgaard* (2019)](https://joss.theoj.org/papers/10.21105/joss.01763)

--
* Links
 + Stable version: <https://CRAN.R-project.org/package=Ryacas>
 + Development version: <https://github.com/r-cas/ryacas/>
 + Online documentation: <http://r-cas.github.io/ryacas/>

#### `caracas`

* Based on [SymPy](https://www.sympy.org/) (large computer algebra library for Python), possible with `reticulate`

--
* '*cara*': face in Spanish (Castellano) / '*cas*': computer algebra system

--
* Initiated in 2019 by **Søren Højsgaard** and **Mikkel Meyer Andersen**
  + Supported by a grant from the R Consortium

--
* Links
 + Stable version: <https://CRAN.R-project.org/package=caracas>
 + Development version: <https://github.com/r-cas/caracas/>
 + Online documentation: <http://r-cas.github.io/caracas/>

---

### Pros/cons

#### `Ryacas`

* Tight two-way integration
* "Easily" extensible
* Not as feature rich and well-tested as SymPy

#### `caracas`

* `SymPy` is a large project with many developers
* Communication is back-and-forth between R and Python via `reticulate`

#### Similarities

* Both works by creating symbols (`Ryacas::ysym()`/`caracas::symbol()`) and 
  use S3 generics to use these as normal R objects.

*We focus on `caracas` in this talk and in the future developments.*

---

# Examples

---

### Getting started

```r
#devtools::install_github("r-cas/caracas")
install.packages("caracas")
```

```r
library(caracas)
packageVersion("caracas")
```

```
## [1] '1.0.1'
```

* Calculus: derivatives, integrals, sums etc.
* Linear algebra
* Solving equations

---

### Symbols

```r
*x <- symbol("x")
y <- symbol("y")
eq <- x^2 + 3*x + 4*y + y^4
```

```r
der(eq, x) # derivative
```

```
## [caracas]: 2⋅x + 3
```

```r
H <- der2(eq, c(x, y)) # Hessian
# print() / as.character() / tex() : $$H = `r tex(H)`$$ 
```

`$$H = \left[\begin{matrix}2 & 0\\0 & 12 y^{2}\end{matrix}\right]$$`

```r
eval(as_r(H), list(x = 1, y = 1))
```

```
##      [,1] [,2]
## [1,]    2    0
## [2,]    0   12
```

---

### Sums

```
sumf(expr, var, [from, to], doit = TRUE)
```

```r
k <- symbol("k")
res <- sumf(k, k, 0, "n", doit = FALSE)
res2 <- doit(res)
```

```r
# $$`r tex(res)` = `r tex(res2)`$$
```

`$$\sum_{k=0}^{n} k = \frac{n^{2}}{2} + \frac{n}{2}$$`

```r
simplify(res2)
```

```
## [caracas]: n⋅(n + 1)
##            ─────────
##                2
```

---

### Integrals

```
intf(expr, var, [from, to], doit = TRUE)
```

```r
res <- intf(1/x, x)
res
```

```
## [caracas]: log(x)
```

```r
res <- intf(1/x, x, doit = FALSE)
tex(res)
```

```
## [1] "\\int \\frac{1}{x}\\, dx"
```

```r
doit(res)
```

```
## [caracas]: log(x)
```

```r
as_r(res) # doit() is automatically called
```

```
## expression(log(x))
```

---

### Limits

```
limf(expr, [var, to], dir = "-", doit = TRUE)
```

```r
limf(sin(x)/x, x, 0)
```

```
## [caracas]: 1
```

```r
limf(1/x, x, 0)
```

```
## [caracas]: ∞
```

```r
l <- limf(1/x, x, 0, dir = "-", doit = FALSE)
tex(l)
```

```
## [1] "\\lim_{x \\to 0^-} \\frac{1}{x}"
```

```r
doit(l)
```

```
## [caracas]: -∞
```

---

### Solving equations

```
solve_lin(A, b)           # Ax = b; also inv(A) for inverse of A
solve_sys(lhs, rhs, vars) # lhs = rhs for vars; rhs omitted finds roots
```

```r
lhs <- cbind(3*x*y - y, x); rhs <- cbind(-5*x, y+4)
```

`$$\begin{align}
3 x y - y &= - 5 x \\
x &= y + 4
\end{align}$$`

```r
sol <- solve_sys(lhs, rhs, c(x, y))
sol
```

```
## Solution 1:
##   x =  2/3 
##   y =  -10/3 
## Solution 2:
##   x =  2 
##   y =  -2
```

```r
sol[[1]]$x # maybe with as_r(...)
```

```
## [caracas]: 2/3
```

---

### Linear algebra

```
as_symbol() # Converts R object to caracas symbol
```

```r
A <- as_symbol(matrix(c(2, 1, 4, "x"), 2, 2))
A
```

```
## [caracas]: ⎡2  4⎤
##            ⎢    ⎥
##            ⎣1  x⎦
```

```r
t(A)
```

```
## [caracas]: ⎡2  1⎤
##            ⎢    ⎥
##            ⎣4  x⎦
```

```r
A[1, ]
```

```
## [caracas]: [2  4]ᵀ
```

---

```r
determinant(A)
```

```
## [caracas]: 2⋅x - 4
```

```r
Ai <- inv(A) # or solve_lin(A)
Ai
```

```
## [caracas]: ⎡  2⋅x      -8   ⎤
##            ⎢───────  ───────⎥
##            ⎢4⋅x - 8  4⋅x - 8⎥
##            ⎢                ⎥
##            ⎢  -1        2   ⎥
##            ⎢───────  ───────⎥
##            ⎣2⋅x - 4  2⋅x - 4⎦
```

```r
simplify(Ai)
```

```
## [caracas]: ⎡    x       -2  ⎤
##            ⎢─────────  ─────⎥
##            ⎢2⋅(x - 2)  x - 2⎥
##            ⎢                ⎥
##            ⎢   -1        1  ⎥
##            ⎢ ───────   ─────⎥
##            ⎣ 2⋅x - 4   x - 2⎦
```

---

## Solving linear systems of equations

```r
b <- as_symbol(c("x", 3))
y <- A %*% b
y
```

```
## [caracas]: [2⋅x + 12  4⋅x]ᵀ
```

```r
solve_lin(A, y)
```

```
## [caracas]: ⎡2⋅x⋅(2⋅x + 12)     32⋅x     8⋅x     2⋅x + 12⎤
##            ⎢────────────── - ───────  ─────── - ────────⎥
##            ⎣   4⋅x - 8       4⋅x - 8  2⋅x - 4   2⋅x - 4 ⎦ᵀ
```

```r
solve_lin(A, y) %>% simplify() # == b
```

```
## [caracas]: [x  3]ᵀ
```

---

# Example: Multinomial likelihood

---

### Multinomial likelihood

The multinomial likelihood for three categories is

`$$\text{lik}(p \mid y) \propto p_1^{y_1} p_2^{y_2} p_3^{y_3}.$$`

Taking `$\log$` we get (together with constraint function `$g$`)

`\begin{align}
  l(p) &= y_1 \log(p_1) + y_2 \log(p_2) + y_3 \log(p_3) \\
  g(p) &= p_1 + p_2 + p_3 - 1 .
\end{align}`

Maximise `$l(p)$` for `$g(p) = 0$` using Lagrange multipliers to get the unconstrainted optimisation problem
`\begin{align}
  L = -l(p) + \lambda g(p)
\end{align}`

```r
p <- as_symbol(c("p1", "p2", "p3"))
y <- as_symbol(c("y1", "y2", "y3"))
a <- as_symbol("a")
l <- sum(y*log(p))
L <- -l + a*(sum(p) - 1)
L
```

```
## [caracas]: a⋅(p₁ + p₂ + p₃ - 1) - y₁⋅log(p₁) - y₂⋅log(p₂) - y₃⋅log(p₃)
```

---

### Multinomial likelihood

```r
gL <- der(L, c(p, a)) # gradient
sol <- solve_sys(gL, c(p, a))
sol
```

```
## Solution 1:
##   p1 =       y₁     
##        ────────────
##        y₁ + y₂ + y₃ 
##   p2 =       y₂     
##        ────────────
##        y₁ + y₂ + y₃ 
##   p3 =       y₃     
##        ────────────
##        y₁ + y₂ + y₃ 
##   a  =  y₁ + y₂ + y₃
```

---
class: inverse, center, middle

# Example: 1-way ANOVA

---

## 1-way ANOVA

Design matrix for three groups each with two observations:

```r
groups <- 1:3
obs_in_each_group <- 2
f <- factor(rep(groups, each = obs_in_each_group))
X <- as_symbol(model.matrix(~ f))
X # Design matrix
```

```
## [caracas]: ⎡1  0  0⎤
##            ⎢       ⎥
##            ⎢1  0  0⎥
##            ⎢       ⎥
##            ⎢1  1  0⎥
##            ⎢       ⎥
##            ⎢1  1  0⎥
##            ⎢       ⎥
##            ⎢1  0  1⎥
##            ⎢       ⎥
##            ⎣1  0  1⎦
```

---

## 1-way ANOVA

```r
y <- as_symbol(matrix(paste0("y", seq_along(f))))
*beta_hat <- inv(t(X) %*% X) %*% t(X) %*% y
print(beta_hat, rowvec = FALSE)
```

```
## [caracas]: ⎡      y₁   y₂      ⎤
##            ⎢      ── + ──      ⎥
##            ⎢      2    2       ⎥
##            ⎢                   ⎥
##            ⎢  y₁   y₂   y₃   y₄⎥
##            ⎢- ── - ── + ── + ──⎥
##            ⎢  2    2    2    2 ⎥
##            ⎢                   ⎥
##            ⎢  y₁   y₂   y₅   y₆⎥
##            ⎢- ── - ── + ── + ──⎥
##            ⎣  2    2    2    2 ⎦
```

---
class: inverse, center, middle

# Thank you!

Mikkel Meyer Andersen (<mikl@math.aau.dk>) 
 
Søren Højsgaard (<sorenh@math.aau.dk>)

## References

`caracas` and `Ryacas`: <https://github.com/r-cas> and CRAN