---
title: "Functions `contract()` and `contract_elementary()` in the `stokes` package"
author: "Robin K. S. Hankin"
output: html_vignette
bibliography: stokes.bib
link-citations: true
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{contract}
  %\usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
set.seed(0)
library("stokes")
library("spray")
library("disordR")
library("magrittr")
knitr::opts_chunk$set(echo = TRUE)
options(rmarkdown.html_vignette.check_title = FALSE)
knit_print.function <- function(x, ...){dput(x)}
registerS3method(
  "knit_print", "function", knit_print.function,
  envir = asNamespace("knitr")
)
```

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

```{r,label=showcontract,comment=""}
contract
contract_elementary
```

To cite the `stokes` package in publications, please use
@hankin2022_stokes.  Given a $k$-form $\phi\colon
V^k\longrightarrow\mathbb{R}$ and a vector $\mathbf{v}\in V$, the
_contraction_ $\phi_\mathbf{v}$ of $\phi$ and $\mathbf{v}$ is a
$k-1$-form with

\[
  \phi_\mathbf{v}\left(\mathbf{v}^1,\ldots,\mathbf{v}^{k-1}\right) =
  \phi\left(\mathbf{v},\mathbf{v}^1,\ldots,\mathbf{v}^{k-1}\right)
\]

provided $k>1$; if $k=1$ we specify
$\phi_\mathbf{v}=\phi(\mathbf{v})$.  If @spivak1965 does discuss this,
I have forgotten it.  Function `contract_elementary()` is a low-level
helper function that translates elementary $k$-forms with coefficient
1 (in the form of an integer vector corresponding to one row of an
index matrix) into its contraction with $\mathbf{v}$; function
`contract()` is the user-friendly front end.  We will start with some
simple examples.  I will use `phi` and $\phi$ to represent the same
object.

```{r label=simpleexample}
(phi <- as.kform(1:5))
```

Thus $k=5$ and we have $\phi=\mathrm{d}x^1\wedge\mathrm{d}x^2\wedge
\mathrm{d}x^3\wedge\mathrm{d}x^4\wedge\mathrm{d}x^5$.  We have that
$\phi$ is a linear alternating map with

$$\phi\left(\begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix},
\begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix},
\begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix} \right)=1$$.

The contraction of $\phi$ with any vector $\mathbf{v}$ is thus a
$4$-form mapping $V^4$ to the reals with
$\phi_\mathbf{v}\left(\mathbf{v}^1,\mathbf{v}^2,\mathbf{v}^3,\mathbf{v}^4\right)=\phi\left(\mathbf{v},\mathbf{v}^1,\mathbf{v}^2,\mathbf{v}^3,\mathbf{v}^4\right)$.
Taking the simplest case first, if $\mathbf{v}=(1,0,0,0,0)$ then

```{r label=contract10000}
v <- c(1,0,0,0,0)
contract(phi,v)
```

that is, a linear alternating map from $V^4$ to the reals with

$$\phi_\mathbf{v}\left(
\begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix},
\begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=1$$.

(the contraction has the first argument of $\phi$ understood to be
$\mathbf{v}=(1,0,0,0,0)$).  Now consider $\mathbf{w}=(0,1,0,0,0)$:

```{r label=contract01000}
w <- c(0,1,0,0,0)
contract(phi,w)
```

$$\phi_\mathbf{w}\left(
\begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix},
\begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=1
\qquad\mbox{or}\qquad
\phi_\mathbf{w}\left(
\begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix},
\begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix},
\begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=-1$$.

Contraction is linear, so we may use more complicated vectors:

```{r complicatedvectors}
contract(phi,c(1,3,0,0,0))
contract(phi,1:5)
```

We can check numerically that the contraction is linear in its vector
argument: $\phi_{a\mathbf{v}+b\mathbf{w}}=
a\phi_\mathbf{v}+b\phi_\mathbf{w}$.


```{r label=verifylinearityinv,cache=TRUE}
a <- 1.23
b <- -0.435
v <- 1:5
w <- c(-3, 2.2, 1.1, 2.1, 1.8)

contract(phi,a*v + b*w) == a*contract(phi,v) + b*contract(phi,w)
```

We also have linearity in the alternating form:
$(a\phi+b\psi)_\mathbf{v}=a\phi_\mathbf{v} + b\psi_\mathbf{v}$.

```{r label=verifylinearityinphi}
(phi <- rform(2,5))
(psi <- rform(2,5))
a <- 7
b <- 13
v <- 1:7
contract(a*phi + b*psi,v) == a*contract(phi,v) + b*contract(psi,v)
```

## Contraction of products

@weintraub2014 gives us the following theorem, for any $k$-form $\phi$
and $l$-form $\psi$:

\[
\left(\phi\wedge\psi\right)_\mathbf{v} = 
\phi_\mathbf{v}\psi + (-1)^k\phi\wedge\psi_\mathbf{v}.\]

We can verify this numerically with $k=4,l=5$:

```{r label=verifycross}
phi <- rform(terms=5,k=3,n=9)
psi <- rform(terms=9,k=4,n=9)
v <- sample(1:100,9)
contract(phi^psi,v) ==  contract(phi,v) ^ psi - phi ^ contract(psi,v)
```

The theorem is verified.  We note in passing that the object itself is quite complicated:

```{r,label=quitecomplicated}
summary(contract(phi^psi,v))
```

We may also switch $\phi$ and $\psi$, remembering to change the sign:

```{r, label=signswitch}
contract(psi^phi,v) ==  contract(psi,v) ^ phi + psi ^ contract(phi,v)
```


## Repeated contraction


It is of course possible to contract a contraction.  If $\phi$ is a
$k$-form, then $\left(\phi_\mathbf{v}\right)_\mathbf{w}$ is a $k-2$
form with

$$
\left(\phi_\mathbf{u}\right)_\mathbf{v}\left(\mathbf{w}^1,\ldots,\mathbf{w}^{k-2}\right)=\phi\left(\mathbf{u},\mathbf{v},\mathbf{w}^1,\ldots,\mathbf{w}^{k-2}\right)
$$

And this is straightforward to realise in the package:

```{r label=straight}
(phi <- rform(2,5))
u <- c(1,3,2,4,5,4,6)
v <- c(8,6,5,3,4,3,2)
contract(contract(phi,u),v)
```

But `contract()` allows us to perform both contractions in one
operation: if we pass a matrix $M$ to `contract()` then this is
interpreted as repeated contraction with the columns of $M$:

```{r bothatonce,cache=TRUE}
M <- cbind(u,v)
contract(contract(phi,u),v) == contract(phi,M)
```

We can verify directly that the system works as intended.  The lines
below strip successively more columns from argument `V` and contract
with them:

```{r verifydirect,cache=TRUE}
(o <- kform(spray(t(replicate(2, sample(9,4))), runif(2))))
V <- matrix(rnorm(36),ncol=4)
jj <- c(
   as.function(o)(V),
   as.function(contract(o,V[,1,drop=TRUE]))(V[,-1]), # scalar
   as.function(contract(o,V[,1:2]))(V[,-(1:2),drop=FALSE]),
   as.function(contract(o,V[,1:3]))(V[,-(1:3),drop=FALSE]),
   as.function(contract(o,V[,1:4],lose=FALSE))(V[,-(1:4),drop=FALSE])
)
print(jj)
max(jj) - min(jj) # zero to numerical precision
```

and above we see agreement to within numerical precision.  If we pass
three columns to `contract()` the result is a $0$-form:

```{r label=getazeroform}
contract(o,V)
```

In the above, the result is coerced to a scalar which is returned in
the form of a `disord` object; in order to work with a formal $0$-form
(which is represented in the package as a `spray` with a zero-column
index matrix) we can use the `lost=FALSE` argument:

```{r doanothercontractnolose}
contract(o,V,lose=FALSE)
```

thus returning a $0$-form.  If we iteratively contract a
$k$-dimensional $k$-form, we return the determinant, and this may be
verified as follows:

```{r,label=verifydeterminant}
o <- as.kform(1:5)
V <- matrix(rnorm(25),5,5)
LHS <- det(V)
RHS <- contract(o,V)
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)
```

Above we see agreement to within numerical error.

## Contraction from first principles

Suppose we wish to contract $\phi=dx^{i_1}\wedge\cdots\wedge dx^{i_k}$
with vector $\mathbf{v}=(v_1\mathbf{e}_1,\ldots,v_k\mathbf{e}_k)$.
Thus we seek $\phi_\mathbf{v}$ with
$\phi_\mathbf{v}\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right) =
dx^{i_1}\wedge\cdots\wedge
dx^{i_k}\left(\mathbf{v},\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)$.
Writing $\mathbf{v}=v_1\mathbf{e}_1+\cdots+\mathbf{e}_k$, we have

\begin{eqnarray}
\phi_\mathbf{v}\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right) &=&
dx^{i_1}\wedge\cdots\wedge
dx^{i_k}\left(\mathbf{v},\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)\\&=& 
dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(v_1\mathbf{e}_1+\cdots+v_k\mathbf{e}_k,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)\\&=& 
v_1 dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(\mathbf{e}_1,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)+\cdots+
v_k dx^{i_1}\wedge\cdots\wedge dx^{i_k}\left(\mathbf{e}_k,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right).
\end{eqnarray}

where we have exploited linearity.  To evaluate this it is easiest and
most efficient  to express  $\phi$ as  $\bigwedge_{j=1}^kdx^{i_j}$ and
cycle  through the  index $j$,  and  use various  properties of  wedge
products:

\begin{eqnarray}
dx^{i_1}\wedge\cdots\wedge dx^{i_k} &=&
(-1)^{j-1} dx^{i_j}\wedge\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{i_j}}\wedge\cdots\wedge dx^{i-k}\right)\\ &=&
(-1)^{j-1} k\operatorname{Alt}\left(dx^{i_j}\otimes\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{i_j}}\wedge\cdots\wedge dx^{i-k}\right)\right)
\end{eqnarray}

(above, a hat indicates a term's being omitted).  From this, we see
that $l\not\in L\longrightarrow\phi=0$ (where $L=\left\lbrace
i_1,\ldots i_k\right\rbrace$ is the index set of $\phi$), for all the
$dx^p$ terms kill $\mathbf{e}_l$.  On the other hand, if $l\in L$ we
have

\begin{eqnarray}
\phi_{\mathbf{e}_l}\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right) &=&
\left(dx^{l}\wedge\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{l}}\wedge\cdots\wedge
dx^{i_k}\right)\right)\left(\mathbf{e}_l,\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)\\ &=&
(-1)^{l-1}k\left(dx^{l}\otimes\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{l}}\wedge\cdots\wedge
dx^{i_k}\right)\right)\left(\mathbf{e}_l,\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)\right)\\ &=&
(-1)^{l-1}k\left(dx^{i_1}\wedge\cdots\wedge\widehat{dx^{l}}\wedge\cdots\wedge
dx^{i_k}\right)\left(\mathbf{v}_1,\ldots,\mathbf{v}_{k-1}\right)
\end{eqnarray}

## Worked example using `contract_elementary()`

Function `contract_elementary()` is a bare-bones low-level no-frills
helper function that returns $\phi_\mathbf{v}$ for $\phi$ an
elementary form of the form $dx^{i_1}\wedge\cdots\wedge dx^{i_k}$.
Suppose we wish to contract $\phi=dx^1\wedge dx^2\wedge dx^5$ with
vector $\mathbf{v}=(1,2,10,11,71)^T$.

Thus we seek $\phi_\mathbf{v}$ with
$\phi_\mathbf{v}\left(\mathbf{v}_1,\mathbf{v}_2 \right)=dx^1\wedge
dx^2\wedge dx^5\left(\mathbf{v},\mathbf{v}_1,\mathbf{v}_2\right)$.
Writing $\mathbf{v}=v_1\mathbf{e}_1+\cdots+v_5\mathbf{e}_5$ we have

\begin{eqnarray}
\phi_\mathbf{v}\left(\mathbf{v}_1,\mathbf{v}_2 \right)                         &=&
dx^1\wedge dx^2\wedge dx^5\left(\mathbf{v},\mathbf{v}_1,\mathbf{v}_2\right)\\  &=&
dx^1\wedge dx^2\wedge dx^5\left(v_1\mathbf{e}_1+\cdots+v_5\mathbf{e}_5,\mathbf{v}_1,\mathbf{v}_2\right)\\&=&
v_1 dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_1,\mathbf{v}_1,\mathbf{v}_2\right)+
v_2 dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_2,\mathbf{v}_1,\mathbf{v}_2\right)\\
&{}&\qquad +v_3dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_3,\mathbf{v}_1,\mathbf{v}_2\right)+
v_4 dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_4,\mathbf{v}_1,\mathbf{v}_2\right)\\
&{}&\qquad\qquad +v_5dx^1\wedge dx^2\wedge dx^5\left(\mathbf{e}_5,\mathbf{v}_1,\mathbf{v}_2\right)\\&=&
v_1 dx^2\wedge dx^5\left(\mathbf{v}_1,\mathbf{v}_2\right)-
v_2 dx^1\wedge dx^5\left(\mathbf{v}_1,\mathbf{v}_2\right)+0+0+
v_5 dx^1\wedge dx^2\left(\mathbf{v}_1,\mathbf{v}_2\right)
\end{eqnarray}

(above, the zero terms are because the vectors $\mathbf{e}_3$ and
$\mathbf{e}_4$ are killed by $dx^1\wedge dx^2\wedge dx^5$).  We can
see that the way to evaluate the contraction is to go through the
terms of $\phi$ [that is, $dx^1$, $dx^2$, and $dx^5$] in turn, and sum
the resulting expressions:


```{r,fromfirst}
o <- c(1,2,5)
v <- c(1,2,10,11,71)
(
(-1)^(1+1) * as.kform(o[-1])*v[o[1]] + 
(-1)^(2+1) * as.kform(o[-2])*v[o[2]] +
(-1)^(3+1) * as.kform(o[-3])*v[o[3]]
)
```

This is performed more succinctly by `contract_elementary()`:


```{r, label=ceex}
contract_elementary(o,v)
```


## The "meat" of `contract()`

Given a vector `v`, and `kform` object `K`, the meat of `contract()` would be


```
Reduce("+", Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K))))
```


I will show this in operation with simple but nontrivial arguments.


```{r,label=meatcontract}
(K <- as.kform(spray(matrix(c(1,2,3,6,2,4,5,7),2,4,byrow=TRUE),1:2)))
v <- 1:7
```

Then the inside bit would be

```{r,label=insidebitofmeat}
apply(index(K), 1, contract_elementary, v)
```

Above we see a two-element list, one for each elementary term of `K`.
We then use base R's `Map()` function to multiply each one by the appropriate coefficient:


```{r usemap}
Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K)))
```

And finally use `Reduce()` to sum the terms:

```{r usereduce}
Reduce("+",Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K))))
```

However, it might be conceptually easier to use `magrittr` pipes to
give an equivalent definition:

```{r usemagrittr}
K                                %>%
index                              %>%
apply(1,contract_elementary,v)       %>%
Map("*", ., K %>% coeffs %>% elements) %>%
Reduce("+",.)
```

Well it might be clearer to Hadley but frankly YMMV.

```{r tidyup, include=FALSE}
rm(phi)
```


## References