---
title: "Brobdingnagian matrices"
author: "Robin K. S. Hankin"
output: html_vignette
bibliography: brob.bib
vignette: >
  %\VignetteIndexEntry{brobmat}
  %\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/Brobdingnag.png", package = "Brobdingnag"))
```

To cite the Brobdingnag package in publications, please use
@hankin2007_brob.  R package `Brobdingnag` allows for very large
numbers, and the package now has basic functionality for matrices with
Brobdingnagian entries: objects of class `brobmat`.  It includes
matrix multiplication and addition, but determinants and matrix
inverses are not implemented.  First load the package:

```{r loadpkg}
library("Brobdingnag")
```

The standard way to create a Brobdgingnagian matrix (a
`brobmat`) is to use function `brobmat()` which takes arguments
similar to `matrix()` and returns a matrix of entries created with
`brob()`:

```{r firstbrobmat}
M1 <- brobmat(-10:13,4,6)
colnames(M1) <- state.abb[1:6]
M1
```

Above, note that all entries of `M1` are greater than zero; `M1[1,1]`,
for example, is `exp(-10)`, or $e^{-10}\simeq 4.54\times 10^{-5}$.
For negative matrix entries, function `brobmat()` takes a Boolean
argument `positive` that specifies the sign:

```{r secondbrobmat}
M2 <- brobmat(
c(1,104,-66,45,1e40,-2e40,1e-200,232.2),2,4,
positive=c(T,F,T,T,T,F,T,T))
M2
```

Standard matrix arithmetic is implemented, thus:

```{r standardarith}
rownames(M2) <- c("a","b")
colnames(M2) <- month.abb[1:4]
M2
M2[2,3] <- 0
M2
M2+1000
```

We can also do matrix multiplication, although it is slow:


```{r matmult}
M2 %*% M1
```


## Numerical verification: matrix multiplication

We will verify matrix multiplication by carrying out the same
operation in two different ways.  First, create two largish
Brobdingnagian matrices:

```{r verifymatmult}
nrows <- 11
ncols <- 18
M3 <- as.brobmat(matrix(rnorm(nrows*ncols),nrows,ncols))
M4 <- as.brobmat(matrix(rnorm(nrows*ncols),ncols,nrows))
M3[1:3,1:3]
```

Now calculate the matrix product by coercing to numeric matrices and
using base R matrix multiplication:

```{r usebaseR}
p1 <- as.matrix(M3) %*% as.matrix(M4)
```

Secondly, we use Brobdingnagian matrix multiplication, and then
coercing to numeric:

```{r usebrob}
p2 <- as.matrix(M3 %*% M4)
```

The difference

```{r showdiff}
max(abs(p1-p2))
```

is small.  Now the other way:

```{r otherway}
q1 <- M3 %*% M4
q2 <- as.brobmat(as.matrix(M3) %*% as.matrix(M4))
max(abs(as.brob(q1-q2)))
```

Above we see that the difference of `exp(-30.267)` (about $7\times
10^{-14}$), is small.

## Numerical verification: integration with the `cubature` package

The matrix functionality of the `Brobdingnag` package was originally
written to leverage the functionality of the `cubature` package.  Here
I give some numerical verification for this.

Suppose we wish to evaluate

\[
\int_{x=0}^{x=4}(x^2-4)\,dx
\]

using numerical methods.  See how the integrand includes positive and
negative values; the theoretical value is $\frac{16}{3}=5.33\ldots$.
The `cubature` idiom for this would be

```{r,label = numericalintegration}
library("cubature")

f.numeric <- function(x){x^2 - 4}

out.num <- cubature::hcubature(f = f.numeric, lowerLimit = 0, upperLimit = 4, vectorInterface = TRUE)
out.num
```

and the Brobdingnagian equivalent would be


```{r,label = numericalintegrationbrob}

f.brob <- function(x) {
	x <- as.brob(x[1, ])
	as.matrix( brobmat(x^2 - 4, ncol = length(x)))
}

out.brob <- cubature::hcubature(f = f.brob, lowerLimit = 0, upperLimit = 4, vectorInterface = TRUE)
out.brob
```

We may compare the two methods:

```{r,label=comparebrobandnumeric}
out.brob$integral - out.num$integral
```

### References {-}