---
title: 'The `mvp` package: fast multivariate polynomials in R'
author: "Robin K. S. Hankin"
output: html_vignette
bibliography: mvp.bib
link-citations: true
vignette: >
  %\VignetteIndexEntry{mvp}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r set-options, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 3.5, message = FALSE, warning = FALSE)
options(width = 80, tibble.width = Inf)
```

```{r out.width='20%', out.extra='style="float:right; padding:10px"',echo=FALSE}
knitr::include_graphics(system.file("help/figures/mvp.png", package = "mvp"))
```

To cite the `mvp` package in publications please use @hankin2022_mvp.
Multivariate polynomials are interesting and useful objects, with a
wide range of applications.  The `mvp` package provides some
functionality for fast manipulation of multivariate polynomials, using
the Standard Template library of `C++` [@stroustrup1997], commonly
known as the `STL`.  It is comparable in speed to the `spray` package
for sparse arrays, while retaining the symbolic capabilities of the
`mpoly` package [@kahle2013].  I present some timing results
separately, in `inst/timings.Rmd`.  The `mvp` package uses the
excellent print and coercion methods of `mpoly`.  The `mvp` package
provides improved speed over `mpoly`, the ability to handle negative
powers, and a more sophisticated substitution mechanism.

# The `STL map` class

A `map` is a sorted associative container that contains key-value
pairs with unique keys.  It is interesting here because search and
insertion operations have logarithmic complexity.  Multivariate
polynomials are considered to be the sum of a finite number of
*terms*, each multiplied by a coefficient.  A *term* is something like
$x^2y^3z$.  We may consider this term to be the map

```
{"x" -> 2, "y" -> 3, "z" -> 1}
```

where the map takes symbols to their (integer) power; it is understood
that powers are nonzero.  An `mvp` object is a map from terms to
their coefficients; thus $7xy^2 -3x^2yz^5$ would be

```
{{"x" -> 1, "y" -> 2} -> 7, {"x" -> 2, 'y" -> 1, "z" ->5} -> -3}
```  

and we understand that coefficients are nonzero.  In `C++` the
declarations would be

```
typedef vector <signed int> mypowers;  
typedef vector <string> mynames;  

typedef map <string, signed int> term; 
typedef map <term, double> mvp; 
```

Thus a `term` maps a string to a (signed) integer, and a `mvp`
maps terms to doubles.  One reason why the `map` class is fast is
that the order in which the keys are stored is undefined: the compiler
may store them in the order which it regards as most propitious.  This
is not an issue for the maps considered here as addition and
multiplication are commutative and associative.

Note also that constant terms are handled with no difficulty
(constants are simply maps from the empty map to its value), as is the
zero polynomial (which is simply an empty map).

## The package in use

Consider a simple multivariate polynomial $3xy+z^3+xy^6z$ and its
representation in the following R session:


```{r simplemult}
library("mvp",quietly=TRUE)
library("magrittr")
(p <- as.mvp("3 x y + z^3 + x y^6 z"))
```

Coercion and printing are accomplished by the `mpoly` package
(there is no way I could improve upon Kahle's work).  Note carefully
that the printed representation of the mvp object is created by the
`mpoly` package and the print method can rearrange both the terms
of the polynomial ($3xy+z^3+xy^6z = z^3+3xy+xy^6z$, for example) and
the symbols within a term ($3xy=3yx$, for example) to display the
polynomial in a human-friendly form.

However, note carefully that such rearranging does not affect the
mathematical properties of the polynomial itself.  In the `mvp`
package, the order of the terms is not preserved (or even defined) in
the internal representation of the object; and neither is the order of
the symbols within a single term.  Although this might sound odd, if
we consider a marginally more involved situation, such as

```{r stoatgoat}
(M <- as.mvp("3 stoat goat^6 -4 + 7 stoatboat^3 bloat -9 float boat goat gloat^6"))
dput(M)
```

it is not clear that any human-discernible ordering is preferable to
any other, and we would be better off letting the compiler decide a
propitious ordering.  In any event, the `mpoly` package can
specify a print order:


```{r printlex}
print(M,order="lex", varorder=c("stoat","goat","boat","bloat","gloat","float","stoatboat"))
```

Note in passing that variable names may be any character string, not
just single letters.

## Arithmetic operations

The arithmetic operations `*`, `+`, `-` and `^` work
as expected:

```{r arithasexpected}
(S1 <- rmvp(5,2,2,4))
(S2 <- rmvp(5,2,2,4))
S1 + S2
S1 * S2
S1^2
```

## Substitution

The package has two substitution functionalities.  Firstly, we can
substitute one or more variables for a numeric value.  Define a mvp
object:

```{r s3subs}
(S3 <- as.mvp("x + 5 x^4 y + 8 y^2 x z^3"))
```

And then we may substitute $x=1$:

```{r subsxequal1}
subs(S3, x = 1)
```

Note the natural R idiom, and that the return value is another mvp
object.  We may substitute for the other variables:

```{r naturalidiom}
subs(S3, x = 1, y = 2, z = 3)
```

(in this case, the default behaviour is to return the resulting
polynomial coerced to a scalar).  We can suppress the coercion using
the `drop` (formerly `lose`) argument:

```{r suppresscoercion}
subs(S3, x = 1, y = 2, z = 3, drop=FALSE)
```

The idiom
also allows one to substitute a variable for an `mvp` object:

```{r subavar}
subs(as.mvp("a+b+c"), a="x^6")
```

Note carefully that `subs()` depends on the order of substitution:

```{r subsorder}
subs(as.mvp("a+b+c"), a="x^6",x="1+a")
subs(as.mvp("a+b+c"), x="1+a",a="x^6")
```
  
### Pipes

Substitution works well with pipes:

```{r workswellwithpipes}
as.mvp("a+b") %>% subs(a="a^2+b^2") %>% subs(b="x^6")
```

## Vectorised substitution

Function `subvec()`allows one to substitute variables for numeric
values using vectorised idiom:

```{r subvec_example}
p <- rmvp(6,2,2,letters[1:3])
p
subvec(p,a=1,b=2,c=1:5)   # supply a named list of vectors
```	


## Differentiation

Differentiation is implemented.  First we have the `deriv()`
method:

```{r diffisimp}
(S <- as.mvp("a + 5 a^5*b^2*c^8 -3 x^2 a^3 b c^3"))
deriv(S, letters[1:3])
deriv(S, rev(letters[1:3]))  # should be the same.
```

Also a slightly different form: `aderiv()`, here used to evaluate
$\frac{\partial^6S}{\partial a^3\partial b\partial c^2}$:

```{r aderivisalso}
aderiv(S, a = 3, b = 1, c = 2)
```


Again, pipes work quite nicely:


```{r pipenice}
S %<>% aderiv(a=1,b=2) %>% subs(c="x^4") %>% `+`(as.mvp("o^99"))
S
```

## Taylor series

The package includes functionality to deal with Taylor and Laurent series:


```{r  taylorX}
(X <- as.mvp("1+x+x^2 y")^3)
trunc(X,3)         # truncate, retain only terms with total power <= 3
trunc1(X,x=3)    # truncate, retain only terms with  power of x <= 3
onevarpow(X,x=3) # retain only terms with power of x == 3
```

```{r secondordertaylor}
## second order taylor expansion of f(x)=sin(x+y) for x=1.1, about x=1:
sinxpy <- horner("x+y",c(0,1,0,-1/6,0,+1/120,0,-1/5040))  # sin(x+y)
dx <- as.mvp("dx")
t2 <- sinxpy  + aderiv(sinxpy,x=1)*dx + aderiv(sinxpy,x=2)*dx^2/2
(t2 %<>% subs(x=1,dx=0.1))  # (Taylor expansion of sin(y+1.1), left in symbolic form)
(t2 %>% subs(y=0.3))  - sin(1.4)  # numeric; should be small
```

Function `series()` will decompose an `mvp` object into a power series in a single variable:

```{r mvpseries}
p <- as.mvp("a^2 x b + x^2 a b + b c x^2 + a b c + c^6 x")
p
series(p,'x')
```

This works nicely with `subs()` if we wish to take a power series
about `x-v`, where `v` is any `mvp` object.  For example:

```{r subsxv}
p %>% subs(x="xmv+a+b") %>% series("xmv") 
```

is a series in powers of `x-a-b`.  We may perform a consistency check
by a second substitution, returning us to the original expression:

```{r seriespower}
p == p %>% subs(x="xmv+a+b") %>% subs(xmv="x-a-b")
```

If function `series()` is given a variable name ending in `_m_foo`,
where `foo` is any variable name, then this is typeset as `(x-foo)`.
For example:

```{r anyvar}
as.mvp('x^3 + x*a') %>% subs(x="x_m_a + a") %>% series("x_m_a")
```

So above we see the expansion of $x^2+ax$ in powers of $x-a$.  If we
want to see the expansion of a mvp in terms of a more complicated
expression then it is better to use a nonce variable `v`:

```{r betternonce}
as.mvp('x^2 + x*a+b^3') %>% subs(x="x_m_v + a^2+b") %>% series("x_m_v")
```

where it is understood that $v=a+b^2$.  Function `taylor()` is a
convenience wrapper that does some of the above in one step:

```{r wrapperonestep}
p <- as.mvp("1+x-x*y+a")^2
taylor(p,'x','a')
```

But it's not as good as I expected it to be and frankly it's overkill.

## Extraction

Given a multivariate polynomial, one often needs to extract certain
terms.  Because the terms of an `mvp` object have an
implementation-dependent order, this can be difficult.  But we can use
function `onevarpow()`:


```{r implementdiff}
P <- as.mvp("1 + z + y^2 + x*z^2 +  x*y")^4
onevarpow(P,x=1,y=2)
```



## Negative powers


The `mvp` package handles negative powers, although the idiom is not perfect and I'm still working on it.
There is the `invert()` function:

```{r negpowers}
(p <- as.mvp("1+x+x^2 y"))
invert(p)
```

In the above, `p` is a regular multivariate polynomial which
includes negative powers.  It obeys the same arithmetic rules as other
mvp objects:

```{r samerules}
p + as.mvp("z^6")
```

# The `disordR` package

It is possible to examine the coefficients of an `mvp` object:

```{r introtodisord}
a <- as.mvp("5 + 8*x^2*y - 13*y*x^2 + 11*z - 3*x*yz")
a
coeffs(a)
```

Above, note that the result of `coeffs()` is a `disord` object,
defined in the `disordR` package [@hankin2022_disordR].  The order of
the elements unspecified as the `STL map` class holds the keys and
values in an implementation-specific order.  This device stops the
user from illegal operations on the coefficients.  For example,
suppose we had another `mvp` object, `b`:


```{r disordpluscoeffs,error=TRUE}
b <- a*2
b
coeffs(a) + coeffs(b)
```

above, we get an error because the coefficients of `a` and `b` are
possibly stored in a different order and therefore vector addition
makes no sense.  However, we can operate on coefficients of a single
`mvp` object at will:


```{r disordRsingleobjectstuff}
coeffs(a) > 0
coeffs(a) + coeffs(a)^4
```

Extraction also works but subject to standard `disordR` idiom
restrictions:

```{r disordextract}
coeffs(a)[coeffs(a) > 0]
```

But "mixing" objects is forbidden:


```{r disordnotallowed,error=TRUE}
coeffs(a)[coeffs(b) > 0]
```

Extraction methods work, again subject to `disordR` restrictions:


```{r disordreplace}
coeffs(a)[coeffs(a)<0] <- coeffs(a)[coeffs(a)<0] + 1000 # add 1000 to every negative coefficient
a
```

In cases like this where the replacement object is complicated, using
`magrittr` would simplify the idiom and reduce the opportunity for
error:

```{r usemagrittr}
coeffs(b)[coeffs(b)%%3==1] %<>% `+`(100)  # add 100 to every element equal to 1 modulo 3
b
```

One good use for this is to "zap" small elements:

```{r mvpzapsmall}
x <- as.mvp("1 - 0.11*x + 0.005*x*y")^2
x
```

Then we can zap as follows:

```{r zapmvp}
cx <- coeffs(x)
cx[abs(cx) < 0.01] <- 0
coeffs(x) <- cx
x
```

(I should write a method for `zapsmall()` that does this)


# Multivariate generating functions

We can see the generating function for a chess knight:

```{r knightgener}
knight(2)
```

How many ways are there for a 4D knight to return to its starting
square after four moves?  Answer:

```{r constknight}
constant(knight(4)^4)
```

# References