Reference

x <- symbol('x')
as.character(x)
#> [1] "x"
x
#> [caracas]: x
as_r(x)
#> expression(x)

2*x
#> [caracas]: 2⋅x
y <- symbol('y')
sqrt(3*x^y)
#> [caracas]:       ____
#>                 ╱  y 
#>            √3⋅╲╱  x

z <- cos(x)^2 + sin(x)^2
z
#> [caracas]:    2         2   
#>            sin (x) + cos (x)
simplify(z)
#> [caracas]: 1
tex(z)
#> [1] "\\sin^{2}{\\left(x \\right)} + \\cos^{2}{\\left(x \\right)}"

z <- cos(x)*cos(y) - sin(x)*sin(y)
z
#> [caracas]: -sin(x)⋅sin(y) + cos(x)⋅cos(y)
simplify(z)
#> [caracas]: cos(x + y)
z <- cos(x + y)
z
#> [caracas]: cos(x + y)
expand(z)
#> [caracas]: cos(x + y)
expand_trig(z)
#> [caracas]: -sin(x)⋅sin(y) + cos(x)⋅cos(y)

x <- symbol('x')
y <- symbol('y')
z <- log(x*y)
z
#> [caracas]: log(x⋅y)
expand_log(z)
#> [caracas]: log(x) + log(y)

A <- matrix(c("x", 0, 0, "2*x"), 2, 2)
A
#>      [,1] [,2] 
#> [1,] "x"  "0"  
#> [2,] "0"  "2*x"
B <- as_symbol(A)
B
#> [caracas]: ⎡x   0 ⎤
#>            ⎢      ⎥
#>            ⎣0  2⋅x⎦
2*B
#> [caracas]: ⎡2⋅x   0 ⎤
#>            ⎢        ⎥
#>            ⎣ 0   4⋅x⎦
B*B # Component-wise / Hadamard product
#> [caracas]: ⎡ 2      ⎤
#>            ⎢x    0  ⎥
#>            ⎢        ⎥
#>            ⎢       2⎥
#>            ⎣0   4⋅x ⎦
dim(B)
#> [1] 2 2
sqrt(B)
#> [caracas]: ⎡√x    0  ⎤
#>            ⎢         ⎥
#>            ⎣0   √2⋅√x⎦
log(B)
#> [caracas]: ⎡log(x)    zoo   ⎤
#>            ⎢                ⎥
#>            ⎣ zoo    log(2⋅x)⎦
sum(B)
#> [caracas]: 3⋅x
B %*% t(B)
#> [caracas]: ⎡ 2      ⎤
#>            ⎢x    0  ⎥
#>            ⎢        ⎥
#>            ⎢       2⎥
#>            ⎣0   4⋅x ⎦
diag(B)
#> [caracas]: [x  2⋅x]
cbind(B, B)
#> [caracas]: ⎡x   0   x   0 ⎤
#>            ⎢              ⎥
#>            ⎣0  2⋅x  0  2⋅x⎦
rbind(B, B)
#> [caracas]: ⎡x   0 ⎤
#>            ⎢      ⎥
#>            ⎢0  2⋅x⎥
#>            ⎢      ⎥
#>            ⎢x   0 ⎥
#>            ⎢      ⎥
#>            ⎣0  2⋅x⎦

determinant(B)
#> [caracas]:    2
#>            2⋅x

A <- matrix(c("a", 0, 0, 0, "a", "a", "a", 0, 0), 3, 3)
B <- as_symbol(A)
eigen_val(B)
#> [[1]]
#> [[1]]$eigval
#> [caracas]: a
#> 
#> [[1]]$eigmult
#> [1] 2
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> [caracas]: 0
#> 
#> [[2]]$eigmult
#> [1] 1
eigen_vec(B)
#> [[1]]
#> [[1]]$eigval
#> [caracas]: 0
#> 
#> [[1]]$eigmult
#> [1] 1
#> 
#> [[1]]$eigvec
#> [caracas]: [-1  0  1]ᵀ
#> 
#> 
#> [[2]]
#> [[2]]$eigval
#> [caracas]: a
#> 
#> [[2]]$eigmult
#> [1] 2
#> 
#> [[2]]$eigvec
#> [caracas]: [1  0  0]ᵀ
eigen(eval(as_r(B), list(a = 2)))
#> eigen() decomposition
#> $values
#> [1] 2 2 0
#> 
#> $vectors
#>      [,1]          [,2]       [,3]
#> [1,]    1 -1.000000e+00 -0.7071068
#> [2,]    0  2.220446e-16  0.0000000
#> [3,]    0  2.220446e-16  0.7071068

B
#> [caracas]: ⎡a  0  a⎤
#>            ⎢       ⎥
#>            ⎢0  a  0⎥
#>            ⎢       ⎥
#>            ⎣0  a  0⎦
diag(B)
#> [caracas]: [a  a  0]
diag(B) <- "b"
B
#> [caracas]: ⎡b  0  a⎤
#>            ⎢       ⎥
#>            ⎢0  b  0⎥
#>            ⎢       ⎥
#>            ⎣0  a  b⎦
diag(B)[-2] <- "a"
B
#> [caracas]: ⎡a  0  a⎤
#>            ⎢       ⎥
#>            ⎢0  b  0⎥
#>            ⎢       ⎥
#>            ⎣0  a  a⎦

p <- as_symbol(paste0("p", 1:3))
y <- as_symbol(paste0("y", 1:3))
a <- as_symbol("a")
l <- sum(y*log(p))
l
#> [caracas]: y₁⋅log(p₁) + y₂⋅log(p₂) + 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₃)
tex(L)
#> [1] "a \\left(p_{1} + p_{2} + p_{3} - 1\\right) - y_{1} \\log{\\left(p_{1} \\right)} - y_{2} \\log{\\left(p_{2} \\right)} - y_{3} \\log{\\left(p_{3} \\right)}"
g <- der(L, c(p, a))
g
#> [caracas]: ⎡    y₁      y₂      y₃                  ⎤
#>            ⎢a - ──  a - ──  a - ──  p₁ + p₂ + p₃ - 1⎥
#>            ⎣    p₁      p₂      p₃                  ⎦
sol <- solve_sys(g, 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₃
sol[[1L]]$p1
#> [caracas]:      y₁     
#>            ────────────
#>            y₁ + y₂ + y₃
tex(sol[[1L]]$p1)
#> [1] "\\frac{y_{1}}{y_{1} + y_{2} + y_{3}}"

x <- symbol('x')
eq <- 2*x^2 - x
eq
#> [caracas]:    2    
#>            2⋅x  - x
subs(eq, x, "y")
#> [caracas]:    2    
#>            2⋅y  - y

p <- as_symbol(paste0("p", 1:3))
y <- as_symbol(paste0("y", 1:3))
a <- as_symbol("a")
l <- sum(y*log(p))
L <- -l + a*(sum(p) - 1)
g <- der(L, c(a, p))
sols <- solve_sys(g, c(a, p))
sol <- sols[[1L]]
sol
#> $a
#> [caracas]: y₁ + y₂ + y₃
#> 
#> $p1
#> [caracas]:      y₁     
#>            ────────────
#>            y₁ + y₂ + y₃
#> 
#> $p2
#> [caracas]:      y₂     
#>            ────────────
#>            y₁ + y₂ + y₃
#> 
#> $p3
#> [caracas]:      y₃     
#>            ────────────
#>            y₁ + y₂ + y₃
H <- der2(L, c(p, a))
H
#> [caracas]: ⎡ y₁             ⎤
#>            ⎢───   0    0   1⎥
#>            ⎢  2             ⎥
#>            ⎢p₁              ⎥
#>            ⎢                ⎥
#>            ⎢      y₂        ⎥
#>            ⎢ 0   ───   0   1⎥
#>            ⎢       2        ⎥
#>            ⎢     p₂         ⎥
#>            ⎢                ⎥
#>            ⎢           y₃   ⎥
#>            ⎢ 0    0   ───  1⎥
#>            ⎢            2   ⎥
#>            ⎢          p₃    ⎥
#>            ⎢                ⎥
#>            ⎣ 1    1    1   0⎦
H_sol <- subs_lst(H, sol)
H_sol
#> [caracas]: ⎡              2                                     ⎤
#>            ⎢(y₁ + y₂ + y₃)                                      ⎥
#>            ⎢───────────────         0                0         1⎥
#>            ⎢       y₁                                           ⎥
#>            ⎢                                                    ⎥
#>            ⎢                               2                    ⎥
#>            ⎢                 (y₁ + y₂ + y₃)                     ⎥
#>            ⎢       0         ───────────────         0         1⎥
#>            ⎢                        y₂                          ⎥
#>            ⎢                                                    ⎥
#>            ⎢                                                2   ⎥
#>            ⎢                                  (y₁ + y₂ + y₃)    ⎥
#>            ⎢       0                0         ───────────────  1⎥
#>            ⎢                                         y₃         ⎥
#>            ⎢                                                    ⎥
#>            ⎣       1                1                1         0⎦

A <- matrix(c("a", 0, 0, 0, "a", "a", "a", 0, 0), 3, 3)
B <- as_symbol(A)
B[, 2]
#> [caracas]: [0  a  a]ᵀ
B[, -2]
#> [caracas]: ⎡a  a⎤
#>            ⎢    ⎥
#>            ⎢0  0⎥
#>            ⎢    ⎥
#>            ⎣0  0⎦
B[1, ]
#> [caracas]: [a  0  a]ᵀ
B[1, , drop = FALSE] # Note this is a 1x3 matrix
#> [caracas]: [a  0  a]
B[, 2] <- "x"
B
#> [caracas]: ⎡a  x  a⎤
#>            ⎢       ⎥
#>            ⎢0  x  0⎥
#>            ⎢       ⎥
#>            ⎣0  x  0⎦

sympy <- get_sympy()

sympy$diff("2*a*x", "x")
#> 2*a
sympy$solve("x**2 - 1", "x")
#> [[1]]
#> -1
#> 
#> [[2]]
#> 1

# Multinomial likelihood
p <- as_symbol(paste0("p", 1:3))
y <- as_symbol(paste0("y", 1:3))
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₃)
print(L, ascii = TRUE)
#> [caracas]: a*(p1 + p2 + p3 - 1) - y1*log(p1) - y2*log(p2) - y3*log(p3)
g <- der(L, c(p, a))
sol <- solve_sys(g, 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₃
print(sol, simplify = FALSE)
#> [[1]]
#> [[1]]$p1
#> [caracas]:      y₁     
#>            ────────────
#>            y₁ + y₂ + y₃
#> 
#> [[1]]$p2
#> [caracas]:      y₂     
#>            ────────────
#>            y₁ + y₂ + y₃
#> 
#> [[1]]$p3
#> [caracas]:      y₃     
#>            ────────────
#>            y₁ + y₂ + y₃
#> 
#> [[1]]$a
#> [caracas]: y₁ + y₂ + y₃

as.character(g)
#> [1] "[a - y1/p1, a - y2/p2, a - y3/p3, p1 + p2 + p3 - 1]"
as_character_matrix(g)
#>      [,1]        [,2]         [,3]         [,4]               
#> [1,] "a - y1/p1" " a - y2/p2" " a - y3/p3" " p1 + p2 + p3 - 1"