---
title: "rim - R's interface to Maxima"
output:
  html_document: default
  pdf_document: default
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

Development version of the CRAN package. `rim` provides an interface to the powerful and fairly complete maxima computer algebra system

This repository is a fork from the original `RMaxima` created by Kseniia Shumelchyk and Hans W. Borcher [shumelchyk/rmaxima](https://github.com/shumelchyk/rmaxima), which is currently not maintained.

# Requirements
- Maxima (>= 5.42.1), this package has been tested with version between 5.42.1 and 5.46.0, but may also work with older versios 

# Installation

To install the current release version from CRAN:

```{r, eval = FALSE}
install.packages("rim")
```

If you want to install the latest version the easiest way is to install the R package `remotes` first and install the package from this github repo:

```{r, eval = FALSE}
install.packages("remotes")
remotes::install_github(repo = "rcst/rim")
```

# Usage

If you want to learn how to use Maxima you can check out a number of ressources from it's [project site](https://maxima.sourceforge.io/documentation.html).

## RMarkdown

This section demonstrates the use of `rim` via it's knitr engine. Alternatively, you can run each line of the Maxima code chunks via `maxima.get()`.

This page has been generated using `rim`'s knitr engine itself. 

Note that you can set the output format, e.g. `maxima.options(engine.format = "latex")` to get LaTeX-style expression rendering. When using pandoc to render the resulting file, MathJax is loaded or embedded automatically, when the ouput is a HTML document and so the Maxima output expressions are rendered nicely. Alternatively, expression outputs can be printed in MathML directly using `maxima.options(engine.format = "mathml")`. 

When attaching the package and Maxima is installed already, then it will be set up and started in the background automatically, when a command is send. In addition, if you have installed the R-package `knitr` (which get's installed by default with this package as well), it will register maxima as a `knitr` engine. If this worked, maxima code can be used and rendered in RMarkdown documents.

```{r}
# library(rim)
devtools::load_all()
maxima.options(engine.format = "latex", 
	       engine.label = TRUE,
	       inline.format = "latex", 
	       inline.label = FALSE)
```

For example, to generate this page you can download it's source `index.Rmd` from [here](https://raw.githubusercontent.com/rcst/rim/master/docs/index.Rmd) and then use `rmarkdown::render` to generate it.

```{r, eval = FALSE}
download.file(url = "https://raw.githubusercontent.com/rcst/rim/master/docs/index.Rmd", 
	      destfile = "index.Rmd")
rmarkdown::render("index.Rmd")
```

Now we can enter Maxima expression inside code chunks. Note that we need to end each line by `;` or `$` (suppressing output). For example we can type the following code chunk into our RMarkdown document file:

````

```{maxima}`r ''`
L: sqrt(1 - 1/R^2);
assume(R > 0)$
'integrate(x, x, 0, L) = integrate(x, x, 0, L);
```

````

In the above code chunk we define a variable `L` (that depends on another variable `R`). We tell Maximas database to assume `R` being larger than zero and suppress any output resulting from this command. The last line prints an unevaluated integral on the left-hand side of the equal sign and on the right-hand side evaluates that same definit integral. The tick quotation mark tells Maxima not to evaluate the suceeding expression. This, when rendered into a document will be printed as 

```{maxima first-example}
L: sqrt(1 - 1/R^2);
assume(R > 0)$
'integrate(x, x, 0, L) = integrate(x, x, 0, L);
```

The resulting output (if not suppressed) is printed after each line of code, and thus apprears as several. Input and output reference labels that are assign by Maxima are printed into he code chunk. Those can be used to refer to previous commands in other code chunks. For example

```{maxima rhs}
sqrt(rhs(%o3));
```

takes the right-hand side of the result of the equation from input label `%i3`, which is assigned to output label `%o3` and computes the square-root of it.

Of course exercising on such simple computations is of little benefit. The real benefit comes from more complicated expression and the effort that we would need to put if we wanted to typeset the result, such as this

```{maxima integration, output.var = "foo"}
integrate(1 / (1 + x^4), x);
```

`pandoc` automatically renders the LaTeX output format from Maxima by including MathJax JavaScript script. In general, `pandoc` takes care of how mathematical equations delimited by `$$` are rendered.


However, we can also change Maxima's output format to MathML, which works if the output document is a HTML document. 

```{r changing-options}
maxima.options(engine.format = "mathml")
maxima.options()
```

```{maxima simple-mathml}
sqrt(3/4);
```

```{maxima differentiate}
f(x) := e^(x^2)$
diff(f(x), x);
```

```{maxima repeat}
%;
```

Notice, that we can use the symbol `%` to refer to the last expression and that this works across code chunks.

You can also replay other previous expressions by referring to their output labels, as already demonstrated above. Whether or not reference labels are printed via the `knitr` engine has no influence on their existence in Maxima internally. 

`rim` also supports printing Maxima output as inline. This works by putting an `r` inline code chunk and using `maxima.inline()`. Example: `` `r '\x60r maxima.inline("log(%o1);")\x60'` `` results in `r maxima.inline("log(%o1);")`. 

```{r changing-options-back}
maxima.options(engine.format = "latex")
```

Comments can be added to Maxima code as text between `/*` and `*/` (example taken from the [Maxima manual](https://maxima.sourceforge.io/docs/manual/maxima_73.html)).

```{maxima comments}
/* aa is a variable of interest */  aa : 1234;
bb : aa^2; /* Value of bb depends on aa */ 
/* User-defined infix operator */  infix ("Q");
/* Parses same as a b c, not abc */  a/* foo */Q/* bar */c;
/* Comments /* can be nested /* to any depth */ */ */  1 + xyz;
```


### Chunk Options

The following chunk options are currently available:

- `echo`: Should the code chunk be printed (default is `TRUE`). If `eval=FALSE`, individual code lines are printed without without input labels, since these are assigned by Maxima upon evaluation
- `eval`: Should the code chunk be evaluated by Maxima (default is `TRUE`)
- `include`:  whether code chunk and results should be printed (default is `TRUE`). If `FALSE`, code is still evaluated
- `output.var`: A character string for a variable name that caputures the parsed output expression(s) (i.e., translated into R-code), see example below.

**NOTE:** The `results` chunk option is explicitly ignored from chunk header and implicitly controlled via `maxima.options(engine.format = ...)`. 


#### Capturing parsed Maxima output

One of the above Maxima code chunks had the option `output.var = "foo"` set.

````

```{maxima output-var, output.var = "foo"}`r ''`
integrate(1 / (1 + x^4), x);
```

````

This options tells the Maxima engine to captured the Maxima result(s) as parsed R-code (unevaluated expressions) into a (named) list that in this case is assigned to the variable name `foo`. 
Note, that this variable does **not** get assigned into the global environment i.e., `GlobalEnv` of the running R session. 
Instead, depending on whether the code chunk was run interactively (e.g., via RStudio) or as part of a rendering/ knitting process, it will be assigned into a separate environment or into the environment returned from `knitr::knit_gloabal()`. 
In any case, the environment used is the one returned by `rim_global()`.
To make sure rendering and executing works, it's recommended to use the environment returned by `rim_global()`.

The parsed expressions can subsequently be evaluated in R e.g., to evaluate an expression in the context of a data.frame

```{r using-output-var}
df <- data.frame(x = 1:10, y = letters[1:10])
df$z <- eval(rim_global()$foo[["o5"]], envir = df)
df
rim_global()$foo
```

In this way, Maxima can assist us with symolic computations in our statistical modelling in R.

# Demos

## Plots

```{maxima plot2d, fig.cap = "plot2d()", fig.align="center"}
r: (exp(cos(t))-2*cos(4*t)-sin(t/12)^5)$
plot2d([parametric, r*sin(t), r*cos(t), [t,-8*%pi,8*%pi]]);
```

```{maxima plot3d, fig.cap = "plot3d()", fig.align="center"}
plot3d(log (x^2*y^2), [x, -2, 2], [y, -2, 2],[grid, 29, 29],
       [palette, [gradient, red, orange, yellow, green]],
       color_bar, [xtics, 1], [ytics, 1], [ztics, 4],
       [color_bar_tics, 4]);
```

```{maxima draw, fig.cap = "draw()", fig.align="center"}
example1:
  gr3d (title          = "Controlling color range",
        enhanced3d     = true,
        color          = green,
        cbrange        = [-3,10],
        explicit(x^2+y^2, x,-2,2,y,-2,2)) $

example2:
  gr3d (title          = "Playing with tics in colorbox",
        enhanced3d     = true,
        color          = green,
        cbtics         = {["High",10],["Medium",05],["Low",0]},
        cbrange = [0, 10],
        explicit(x^2+y^2, x,-2,2,y,-2,2))$

example3:
  gr3d (title      = "Logarithmic scale to colors",
        enhanced3d = true,
        color      = green,
        logcb      = true,
        logz       = true,
        palette    = [-15,24,-9],
        explicit(exp(x^2-y^2), x,-2,2,y,-2,2))$

draw(
  dimensions = [500,1500],
  example1, example2, example3)$
```

```{maxima draw2d, fig.cap = "draw2d()", fig.align="center"}
draw2d(
  dimensions = [1000, 1000],
  proportional_axes = xy,
  fill_color        = sea_green,
  color             = aquamarine,
  line_width        = 6,
  ellipse(7,6,2,3,0,360))$
``` 

```{maxima draw3d, fig.cap = "draw3d()", fig.align="center"}
draw3d(
   dimensions = [1000, 1000],
   surface_hide      = true,
   axis_3d           = false,
   proportional_axes = xyz,
 
   color             = blue,
   cylindrical(z,z,-2,2,a,0,2*%pi), 
 
   color            = brown,
   cylindrical(3,z,-2,2,az,0,%pi),
 
   color            = green,
   cylindrical(sqrt(25-z^2),z,-5,5,a,0,%pi))$
```

## Distributions

This is basically a reproduced demo for `Sympy` from [Bryan Zhang's Blog](https://brianzhang01.github.io/2018/04/distributions-with-sympy/).

First we define some helper functions:

```{maxima dist-functions}
area(dist) := integrate(dist, x, minf, inf)$
mean(dist) := area(dist*x)$
EX2(dist) := area(dist*x^2)$
variance(dist) := EX2(dist) - mean(dist)^2$
mgf(dist) := area(dist*%e^(x*t))$
```

### Normal Distribution

```{maxima normal}
normal(x) := 
      (2*%pi*sigma^2)^(-1/2) * 
      exp(-(x-mu)^2/(2*sigma^2));

assume(sigma > 0)$

area(normal(x));
mean(normal(x));
variance(normal(x));
mgf(normal(x));
```

### Pareto Distribution

```{maxima pareto}
pareto(x) := (alpha * (x[m])^alpha) / x^(alpha+1);
assume(x[m] > 0)$ 
assume(alpha > 0)$
```


### Laplace Distribution

```{maxima laplace}
laplace(x) := (2*b)^-1 * exp(-abs(x - mu)/b);

load("abs_integrate")$

assume(b > 0)$

area(laplace(x));
mean(laplace(x));
variance(laplace(x));
```

### Exponential Distribution

```{maxima exponential}
expo(x) := unit_step(x) * lambda * exp(-lambda * x);

assume(lambda > 0)$

area(expo(x));
mean(expo(x));
variance(expo(x));
```

## Matrices

```{maxima matrices}
m: matrix([0, 1, a], [1, 0, 1], [1, 1, 0]);
transpose(m);
determinant(m);
f: invert(m), detout;
m . f;
expand(%);
factor(%);
```

## If-then-else

```{maxima x-y}
x: 1234;
y: 2345;
```

```{maxima if-then-else}
if x > y
  then x
  else y;
```
