README.Rmd

---
output:
  md_document:
    variant: markdown_github
---

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

```{r, echo = FALSE}

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#",
  fig.path = "tools/README-",
  dpi = 250
)
```

# mpoly

<!-- badges: start -->
[![CRAN status](https://www.r-pkg.org/badges/version/mpoly)](https://cran.r-project.org/package=mpoly)
[![Travis build status](https://travis-ci.org/dkahle/mpoly.svg?branch=master)]( https://app.travis-ci.com/dkahle/mpoly)
[![AppVeyor build status](https://ci.appveyor.com/api/projects/status/github/dkahle/mpoly?branch=master&svg=true)](https://ci.appveyor.com/project/dkahle/mpoly)
[![Coverage status](https://codecov.io/gh/dkahle/mpoly/branch/master/graph/badge.svg)](https://app.codecov.io/github/dkahle/mpoly?branch=master)
<!-- badges: end -->


## Specifying polynomials

__mpoly__ is a simple collection of tools to help deal with multivariate polynomials _symbolically_ and functionally in R.  Polynomials are defined with the `mp()` function:

```{r mp, echo = TRUE}
library("mpoly")
mp("x + y")

mp("(x + 4 y)^2 (x - .25)")
```

[Term orders](https://en.wikipedia.org/wiki/Lexicographical_order#Monomials) are available with the reorder function:
```{r reordering}
(p <- mp("(x + y)^2 (1 + x)"))

reorder(p, varorder = c("y","x"), order = "lex")

reorder(p, varorder = c("x","y"), order = "glex")
```

Vectors of polynomials (`mpolyList`'s) can be specified in the same way:
```{r mpolyLists}
mp(c("(x+y)^2", "z"))
```

## Polynomial parts

You can extract pieces of polynoimals using the standard `[` operator, which works on its terms:
```{r subsetting}
p[1]

p[1:3]

p[-1]
```

There are also many other functions that can be used to piece apart polynomials, for example the leading term (default lex order):
```{r lt_lc_lm}
LT(p)

LC(p)

LM(p)
```

You can also extract information about exponents:
```{r degs}
exponents(p)

multideg(p)

totaldeg(p)

monomials(p)
```

## Polynomial arithmetic

Arithmetic is defined for both polynomials (`+`, `-`, `*` and `^`)...
```{r arithmetic}
p1 <- mp("x + y")

p2 <- mp("x - y")

p1 + p2

p1 - p2

p1 * p2

p1^2
```

... and vectors of polynomials:
```{r vector-arithmetic}
(ps1 <- mp(c("x", "y")))

(ps2 <- mp(c("2 x", "y + z")))

ps1 + ps2

ps1 - ps2

ps1 * ps2 
```


## Some calculus

You can compute derivatives easily:
```{r calculus}
p <- mp("x + x y + x y^2")

deriv(p, "y")

gradient(p)
```


## Function coercion

You can turn polynomials and vectors of polynomials into functions you can evaluate with `as.function()`.  Here's a basic example using a single multivariate polynomial:
```{r as-function-basic-single}
f <- as.function(mp("x + 2 y")) # makes a function with a vector argument

f(c(1,1))

f <- as.function(mp("x + 2 y"), vector = FALSE) # makes a function with all arguments

f(1, 1)
```

Here's a basic example with a vector of multivariate polynomials:
```{r as-function-basic-vector}
(p <- mp(c("x", "2 y")))

f <- as.function(p) 

f(c(1,1))

f <- as.function(p, vector = FALSE) 

f(1, 1)
```


Whether you're working with a single multivariate polynomial or a vector of them (`mpolyList`), if it/they are actually univariate polynomial(s), the resulting function is vectorized.  Here's an example with a single univariate polynomial.
```{r as-function-vectorized}
f <- as.function(mp("x^2"))

f(1:3)

(mat <- matrix(1:4, 2))

f(mat) # it's vectorized properly over arrays
```

Here's an example with a vector of univariate polynomials:
```{r as-function-vectorized2}
(p <- mp(c("t", "t^2")))

f <- as.function(p)
f(1)

f(1:3)
```

You can use this to visualize a univariate polynomials like this:
```{r packages, message=FALSE}
library("tidyverse"); theme_set(theme_classic())
```
```{r as-function, fig.height=3}
f <- as.function(mp("(x-2) x (x+2)"))
x <- seq(-2.5, 2.5, .1)

qplot(x, f(x), geom = "line")
```


For multivariate polynomials, it's a little more complicated:
```{r as-function-multi}
f <- as.function(mp("x^2 - y^2")) 
s <- seq(-2.5, 2.5, .1)
df <- expand.grid(x = s, y = s)
df$f <- apply(df, 1, f)
qplot(x, y, data = df, geom = "raster", fill = f)
```

Using [tidyverse-style coding](https://www.tidyverse.org) (install tidyverse packages with `install.packages("tidyverse")`), this looks a bit cleaner:
```{r as-function-multi-tidy}
f <- as.function(mp("x^2 - y^2"), vector = FALSE)
seq(-2.5, 2.5, .1) %>% 
  list("x" = ., "y" = .) %>% 
  cross_df() %>% 
  mutate(f = f(x, y)) %>% 
  ggplot(aes(x, y, fill = f)) + 
    geom_raster()
```

## Algebraic geometry

__Grobner bases are no longer implemented in mpoly; they're now in [m2r](https://github.com/coneill-math/m2r).__

```{r grobner}
# polys <- mp(c("t^4 - x", "t^3 - y", "t^2 - z"))
# grobner(polys)
```

Homogenization and dehomogenization:
```{r homogenization}
(p <- mp("x + 2 x y + y - z^3"))

(hp <- homogenize(p))

dehomogenize(hp, "t")

homogeneous_components(p)
```


## Special polynomials

__mpoly__ can make use of many pieces of the __polynom__ and __orthopolynom__ packages with `as.mpoly()` methods.  In particular, many special polynomials are available.

#### [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials)

You can construct [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) as follows:
```{r chebyshev1}
chebyshev(1)

chebyshev(2)

chebyshev(0:5)
```

And you can visualize them:
```{r chebyshev, fig.height=3}
s <- seq(-1, 1, length.out = 201); N <- 5
(chebPolys <- chebyshev(0:N))

df <- as.function(chebPolys)(s) %>% cbind(s, .) %>% as.data.frame()
names(df) <- c("x", paste0("T_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
```

#### [Jacobi polynomials](https://en.wikipedia.org/wiki/Jacobi_polynomials)

```{r jacobi, fig.height=3}
s <- seq(-1, 1, length.out = 201); N <- 5
(jacPolys <- jacobi(0:N, 2, 2))
 
df <- as.function(jacPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("P_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree) +
  coord_cartesian(ylim = c(-25, 25))
```

#### [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials)

```{r legendre, fig.height=3}
s <- seq(-1, 1, length.out = 201); N <- 5
(legPolys <- legendre(0:N))
 
df <- as.function(legPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("P_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
```

#### [Hermite polynomials](https://en.wikipedia.org/wiki/Hermite_polynomials)

```{r hermite, fig.height=3}
s <- seq(-3, 3, length.out = 201); N <- 5
(hermPolys <- hermite(0:N))

df <- as.function(hermPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("He_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
```

#### [(Generalized) Laguerre polynomials](https://en.wikipedia.org/wiki/Laguerre_polynomials)

```{r laguerre, fig.height=3}
s <- seq(-5, 20, length.out = 201); N <- 5
(lagPolys <- laguerre(0:N))

df <- as.function(lagPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("L_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree) +
  coord_cartesian(ylim = c(-25, 25))
```

#### [Bernstein polynomials](https://en.wikipedia.org/wiki/Bernstein_polynomial)

[Bernstein polynomials](https://en.wikipedia.org/wiki/Bernstein_polynomial) are not in __polynom__ or __orthopolynom__ but are available in __mpoly__ with `bernstein()`:
```{r bernstein, fig.height=3}
bernstein(0:4, 4)

s <- seq(0, 1, length.out = 101)
N <- 5 # number of bernstein polynomials to plot
(bernPolys <- bernstein(0:N, N))

df <- as.function(bernPolys)(s) %>% cbind(s, .) %>% as.data.frame
names(df) <- c("x", paste0("B_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)
```

You can use the `bernstein_approx()` function to compute the Bernstein polynomial approximation to a function.  Here's an approximation to the standard normal density:

```{r bernstein-approx, fig.height=3}
p <- bernstein_approx(dnorm, 15, -1.25, 1.25)
round(p, 4)

x <- seq(-3, 3, length.out = 101)
df <- data.frame(
  x = rep(x, 2),
  y = c(dnorm(x), as.function(p)(x)),
  which = rep(c("actual", "approx"), each = 101)
)
qplot(x, y, data = df, geom = "path", color = which)
```



## [Bezier polynomials and curves](https://en.wikipedia.org/wiki/Bézier_curve)

You can construct [Bezier polynomials](https://en.wikipedia.org/wiki/Bézier_curve) for a given collection of points with `bezier()`:
```{r bezier, fig.height=3}
points <- data.frame(x = c(-1,-2,2,1), y = c(0,1,1,0))
(bezPolys <- bezier(points))
```

And viewing them is just as easy:
```{r bezier-plot, fig.height = 3}
df <- as.function(bezPolys)(s) %>% as.data.frame

ggplot(aes(x = x, y = y), data = df) + 
  geom_point(data = points, color = "red", size = 4) +
  geom_path(data = points, color = "red", linetype = 2) +
  geom_path(size = 2)
```

Weighting is available also:

```{r bezier-weighting, fig.height = 3}
points <- data.frame(x = c(1,-2,2,-1), y = c(0,1,1,0))
(bezPolys <- bezier(points))
df <- as.function(bezPolys, weights = c(1,5,5,1))(s) %>% as.data.frame

ggplot(aes(x = x, y = y), data = df) + 
  geom_point(data = points, color = "red", size = 4) +
  geom_path(data = points, color = "red", linetype = 2) +
  geom_path(size = 2)
```

To make the evaluation of the Bezier polynomials stable, `as.function()` has a special method for Bezier polynomials that makes use of [de Casteljau's algorithm](https://en.wikipedia.org/wiki/De_Casteljau%27s_algorithm).  This allows `bezier()` to be used as a smoother:
```{r bezier-smooth, fig.height=3}
s <- seq(0, 1, length.out = 201) 
df <- as.function(bezier(cars))(s) %>% as.data.frame
qplot(speed, dist, data = cars) +
  geom_path(data = df, color = "red")
```



## Other stuff

I'm starting to put in methods for some other R functions:

```{r lm, fig.height=3}
set.seed(1)
n <- 101
df <- data.frame(x = seq(-5, 5, length.out = n))
df$y <- with(df, -x^2 + 2*x - 3 + rnorm(n, 0, 2))

mod <- lm(y ~ x + I(x^2), data = df)
(p <- mod %>% as.mpoly %>% round)
qplot(x, y, data = df) +
  stat_function(fun = as.function(p), colour = 'red')
```

```{r lm-poly, fig.height=3}
s <- seq(-5, 5, length.out = n)
df <- expand.grid(x = s, y = s) %>% 
  mutate(z = x^2 - y^2 + 3*x*y + rnorm(n^2, 0, 3))

(mod <- lm(z ~ poly(x, y, degree = 2, raw = TRUE), data = df))
as.mpoly(mod)
```



## Installation

* From CRAN: `install.packages("mpoly")`

* From Github (dev version): 
```R
# install.packages("devtools")
devtools::install_github("dkahle/mpoly")
```


## Acknowledgements

This material is based upon work partially supported by the National Science Foundation under Grant No. [1622449](https://www.nsf.gov/awardsearch/showAward?AWD_ID=1622449).