---
title: "The `disordR` package: design philosophy and a use-case in multivariate polynomials" 
author: "Robin K. S. Hankin"
date: "`r Sys.Date()`"
bibliography: disordR_arxiv.bib
vignette: >
  %\VignetteIndexEntry{The disordR package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

### Abstract

To cite the `disordR` package in publications please use
@hankin2022_disordR.  This document motivates the concept of
"disordered vector" using coefficients of multivariate polynomials as
represented by associative maps such as the `STL` map class as used by
the `mvp` package.  Values and keys of a map are stored in an
implementation-specific way so certain extraction and replacement
operations should be forbidden.  For example, consider a finite set
`S` of real numbers.  The "first" element of `S` is implementation
specific (because a set does not specify the order of its
elements)...but the maximum value of `S` has a well-defined value, as
does the sum.  The `disordR` package makes it impossible to execute
forbidden operations (such as finding the "first" element), while
allowing transparent R idiom for permitted operations such as finding
the maximum or sum.  Note that one cannot dispense with order
entirely, as we wish to consider keys and values of `STL` objects
separately, but we need to retain the ability to match individual keys
with their corresponding values.

An illustrative R session is given in which the `disordR` package is
used abstractly, without reference to any particular application, and
then shows how the `disordR` package is used in the `mvp` package.
`disordR` is used in the `clifford`, `freealg`, `hyper2`, `mvp`,
`spray`, `stokes`, and `weyl` packages.

### Introduction

In `C++` [@ISO14882], the `STL map` class [@josuttis1999] is an object
that associates a value to each of a set of keys.  Accessing values or
keys of a `map` object is problematic because the value-key pairs are
not stored in a well-defined order.  The situation is applicable to
any package which uses the `map` objects.  Consider, for example, the
`mvp` package which deals with multinomials using `STL` maps.  An
`mvp` object is a map from terms to coefficients, and a map has no
intrinsic ordering: the maps


```
x -> 1, y -> 3, xy -> 3, xy^3 -> 4
```

and

```
xy^3 -> 4, xy -> 3, x -> 1, y -> 3
```  

are the same map and correspond to the same multinomial (symbolically,
$x+3y+3xy+4xy^3=4xy^3+3xy+x+3y$).  Thus the coefficients of the
multinomial might be `c(1,3,3,4)` or `c(4,3,1,3)`, or indeed any
ordering.  Internally, the elements are stored in some order but the
order used is implementation-specific.  Quite often, I am interested
in the coefficients _per se_, without consideration of their meaning
in the context of a multivariate polynomial.  I might ask:

* "How many coefficients are there?"
* "What is the largest coefficient?"
* "Are any coefficients exactly equal to one?"
* "How many coefficients are greater than 2?"

These are reasonable and mathematically meaningful questions because
they have an answer that is independent of the order in which the
coefficients are stored.  Compare a meaningless question: "what is the
second coefficient?".  This is meaningless because of the order
ambiguity discussed above: the answer is at best
implementation-specific, but fundamentally it is a question that one
should not be allowed to ask.

To deal with the coefficients in isolation in R, one might be tempted
to use a multiset.  However, this approach does not allow one to link
the coefficients with the terms.  Suppose I coerce the coefficients to
a multiset object (as per the `sets` package, for example): then it is
impossible to extract the terms with coefficient greater than 2 (which
would be the polynomial $3y+3xy+4xy^3$) because the link between the
coefficients and the terms is not included in the multiset object.
Sensible questions involving this aspect of `mvp` objects might be:

* Give me all terms with coefficients greater than 2
* Give me all terms with positive coefficients
* Give me all terms with integer coefficients

and these questions cannot be answered if the the coefficients are
stored as a multiset (compare inadmissible questions such as "give me
the first three terms").  Further note that replacement methods are
mathematically meaningful, for example:

* Set any term with a negative coefficient to zero
* Add 100 to any coefficient less than 30

Again these operations are reasonable but precluded by multiset
formalism (compare inadmissible replacements: "replace the first two
terms with zero", or "double the last term" would be inadmissible).

***What we need is a system that forbids stupid questions and stupid
operations, while permitting sensible questions and operations***

The `disord` class of the `disordR` package is specificially designed
for this situation.  This class of object has a slot for the
coefficients in the form of a numeric R vector, but also another slot
which uses hash codes to prevent users from misusing the ordering of
the numeric vector.
   
For example, a multinomial `x+2y+3z` might have coefficients
`c(1,2,3)` or `c(3,1,2)`.  Package idiom to extract the coefficients
of a multivariate polynomial `a` is `coeffs(a)`; but this cannot
return a standard numeric vector.  If stored as a numeric vector, the
user might ask "what is the first element?" and this question should
not be asked [and certainly not answered!], because the elements are
stored in an implementation-specific order.  The `disordR` package
uses `disord` objects which are designed to return an error if such
inadmissible questions are asked.  But `disord` objects can answer
admissible questions and perform admissible operations.

Suppose we have two multivariate polynomials, `a` as defined as above
with `a=x+2y+3z` and `b=x+3y+4z`.  Even though the sum `a+b` is
well-defined algebraically, idiom such as `coeffs(a) + coeffs(b)` is
not defined because there is no guarantee that the coefficients of the
two multivariate polynomials are stored in the same order.  We might
have `c(1,2,3)+c(1,3,4)=c(2,5,7)` or `c(1,2,3)+c(1,4,3)=c(2,6,6)`,
with neither being more "correct" than the other.  In the package,
this ambiguity is rendered void: `coeffs(a) + coeffs(b)` will return
an error.  Note carefully that `a+b`---and hence `coeffs(a+b)`---is
perfectly well defined, although the result is subject to the same
ambiguity as `coeffs(a)`.

In the same way, `coeffs(a) + 1:3` is not defined and will return an
error.  Further, idiom such as `coeffs(a) <- 1:3`and `coeffs(a) <-
coeffs(b)` are not defined and will also return an error.  However,
note that

```
coeffs(a) + coeffs(a)
coeffs(a) + coeffs(a)^2
coeffs(a) <- coeffs(a)^2
coeffs(a) <- coeffs(a)^2 + 7
```

are perfectly well defined, with package idiom behaving as expected.

Idiom such as `disord(a) <- disord(a)^2` is OK: one does not need to
know the order of the coefficients on either side, so long as the
order is the same on both sides.  The idiomatic English equivalent
would be: "the coefficient of each term of `a` becomes its square";
note that this operation is insensitive to the order of coefficients.
The whole shebang is intended to make idiom such as `coeffs(a) <-
coeffs(a)%%2` possible, so we can manipulate polynomials over finite
rings, here $Z/2Z$.

The replacement methods are defined so that an expression like
`coeffs(a)[coeffs(a) < 5] <- 0` works as expected; the English idiom
would be "replace any coefficient less than 5 with 0".

To fix ideas, consider a fixed small mvp object:

```{r}
library("mvp")
a <- as.mvp("5 a c^3 + a^2 d^2 f^2 + 4 a^3 b e^3 + 3 b c f + 2 b^2 e^3")
a
```

Extraction presents issues; consider `coeffs(a)<3`.  This object has
Boolean elements but has the same ordering ambiguity as `coeffs(a)`.
One might expect that we could use this to extract elements of
`coeffs(a)`: specifically, those elements less than 5.  We may use
replace methods for coefficients if this makes sense.  Idiom such as

```
coeffs(a)[coeffs(a)<5] <- 4 + coeffs(a)[coeffs(a)<5]
coeffs(a) <- pmax(a,3)
```

is algebraically meaningful and allowed in the package.
Idiomatically: "Add 4 to any element less than 5"; "coefficients
become the parallel maximum of themselves and 3" respectively.
Further note that `coeffs(a) <- rev(coeffs(a))` is disallowed
(although `coeffs(a) <- rev(rev(coeffs(a)))` is meaningful and
admissible).

So the output of `coeffs(x)` is defined only up to an unknown
rearrangement.  The same considerations apply to the output of
`vars()`, which returns a list of character vectors in an undefined
order, and the output of `powers()`, which returns a numeric list
whose elements are in an undefined order.  However, even though the
order of these three objects is undefined individually, their ordering
is jointly consistent in the sense that the first element of
`coeffs(x)` corresponds to the first element of `vars(x)` and the
first element of `powers(x)`.  The identity of this element is not
defined---but whatever it is, the first element of all three accessor
methods refers to it.

Note also that a single term (something like `4a3*b*c^6`) has the same
issue: the variables are not stored in a well-defined order.  This
does not matter because the algebraic value of the term does not
depend on the order in which the variables appear and this term would
be equivalent to `4bc^6*a^3`.

### An R session with the `disordR` package

We will use the `disordR` package to show how the idiom works.

```{r}
library("disordR")
set.seed(0)
a <- rdis()
a
```

Object `a` is a `disord` object but it behaves similarly to a regular
numeric vector in many ways:

```{r}
a^2
a+1/a
```

Above, note how the result has the same hash code as `a`.  Other
operations that make sense are `max()` and `sort()`:

```{r}
max(a)
sort(a)
```

Above, see how the result is a standard numeric vector.  However,
inadmissible operations give an error:

```{r error=TRUE}
a[1]  # asking for the first element is inadmissible
a[1] <- 1000 # also cannot replace the first element
```

Standard R semantics generally work as expected:

```{r}
x <- a + 1/a
x
y <- a*2-9
y
x+y
```

Above, observe that objects `a`, `x` and `y` have the same hash code:
they are "compatible", in `disordR` idiom.  However, if we try to
combine object `a` with another object with different hash, we get
errors:

```{r,error=TRUE}
b <- rdis()
b
a
a+b
```

The error is given because objects `a` and `b` are stored in an
implementation-specific order (we say that `a` and `b` are
_incompatible_).  In the package, many extract and replace methods are
implemented whenever this is admissible:

```{r}
a[a<0.5] <- 0  # round down
a
b[b>0.6] <- b[b>0.6] + 3  # add 3 to every element greater than 0.6
b
```

Usual semantics follow, provided one is careful to maintain the hash
code:

```{r}
d <- disord(1:10)
d
e <- 10 + 3*d - d^2
e
e<4
d[e<4] <- e[e<4]
d
```

Above, the replacement command works because `d` and `e` _and_ `e<4`
[which is a Boolean `disord` object] all have the same hash code.

# An R session with the `mvp` package

The `mvp` package implements multivariate polynomials using the `STL`
map class.  Following commands only work as intended here with `mvp >=
1.0-12`.  Below we see how `disordR` idiom allows mathematically
meaningful operation while suppressing inadmissible ones:

```{r}
library("mvp")
set.seed(0)
a <- rmvp()
b <- rmvp()
a
b
```

Observe that standard multivariate polynomial algebra works:

```{r}
a + 2*b
(a+b)*(a-b) == a^2-b^2   # should be TRUE (expression is quite long)
```

We can extract the coefficients of these polynomials using the
`coeffs()` function:

```{r}
coeffs(a)
coeffs(b)
```

observe that the coefficients are returned as a `disord` object.  We
may manipulate the coefficients of a polynomial in many ways.  We may
do the following things:

```{r}
coeffs(a)[coeffs(a) < 4] <- 0   # set any coefficient of a that is <4 to zero
a
coeffs(b) <- coeffs(b)%%2       # consider coefficients of b modulo 2
b
```

However, many operations which have reasonable idiom are in fact
meaningless and are implicitly prohibited.  For example:


```{r}
x <- rmvp()     # set up new mvp objects x and y
y <- rmvp()
```

Then the following should all produce errors:

```{r error=TRUE}
coeffs(x) + coeffs(y)  # order implementation specific
coeffs(x) <- coeffs(y) # ditto
coeffs(x) <- 1:2       # replacement value not length 1
coeffs(x)[coeffs(x) < 3] <- coeffs(x)[coeffs(y) < 3]
```

## `vars()` and `powers()` return `disord` objects


The `disord()` function takes a list argument, and this is useful for
working with `mvp` objects:

```{r}
(a <- as.mvp("x^2 + 4 - 3*x*y*z"))
vars(a)
powers(a)
coeffs(a)
```

Note that the hash of all three objects is identical, generated from
the polynomial itself (not just the relevant element of the
three-element list that is an `mvp` object).  This allows us to do
some rather interesting things:

```{r}
double <- function(x){2*x}
(a <- rmvp())
pa <- powers(a)
va <- vars(a)
ca <- coeffs(a)
pa[ca<4] <- sapply(pa,double)[ca<4]
mvp(va,pa,ca)
```

Above, `a` was a multivariate polynomial and we doubled the powers of
all variables in terms with coefficients less than 4. Or even:

```{r}
a <- as.mvp("3 + 5*a*b - 7*a*b*x^2 + 2*a*b^2*c*d*x*y -6*x*y + 8*a*b*c*d*x")
a
pa <- powers(a)
va <- vars(a)
ca <- coeffs(a)
va[sapply(pa,length) > 4] <- sapply(va,toupper)[sapply(pa,length) > 4]
mvp(va,pa,ca)
```

Above, we took multivariate polynomial `a` and replaced the variable
names in every term with more than four variables with their uppercase
equivalents.

# References