blog-static/content/blog/linear_multistep.md

---
title: "Generic Linear Multistep Method Evaluator using Chapel"
date: 2022-07-31T11:40:05-07:00
draft: true
tags: ["Chapel", "Mathematics"]
---

The Chapel programming language has a very pleasant mix of features for
general-purpose programming. Having helped review the recent change
adding [variadic arguments](https://chapel-lang.org/docs/1.14/primers/primers/varargs.html)
to Dyno, the work-in-progress "new compiler" for the language, I've been
excited in giving them a test drive (in the "old" compiler, for the time being).
Chapel can naturally express patterns like "as many variadic arguments as the
first argument specifies". In combination with this, Chapel lets you define
functions that perform computations on _types_ at compile time. This is a lot
of power -- Haskell programmers have done much more, with less!

All that's left is a nail for this unusually large hammer. Ideally, our
application would a family of related algorithms that differ in
predictable ways. Then, we could use type-level computation in Chapel to
pick particular algorithms from this family, without writing much extra
code. Fortunately, I know just the "family of related algorithms":
[linear multistep methods](https://en.wikipedia.org/wiki/Linear_multistep_method)
for computing approximate solutions to ordinary differential equations!
Aside from being a collection of related-but-different algorithms, like
we would want, LMMs are also right up Chapel's alley as a programming language
for high performance computing.

# Euler's Method
Suppose we have a (first-order) differential equation. Such a differential
equation has the form:
{{< latex >}}
y' = f(t, y)
{{< /latex >}}
That is, the value of the derivative of \\(y\\) (aka \\(y'\\)) at any point \\(t\\)
depends also on the value of \\(y\\) itself. You want to figure out what the
value of \\(y\\) is at any point \\(t\\). However, there is no general solution to
such an equation. In the general case, we'll have to settle for approximations.

Here is a simple approximation. Suppose that we knew a starting point for our
\\(y\\). For instance, we could know that at \\(t=0\\), the value of \\(y\\) is
\\(y(t) = y(0) = 1\\). With this, we have both \\(t\\) and \\(y\\), just
what we need to figure out the value of \\(y'(0)\\):
{{< latex >}}
y'(0) = f(0, y(0)) = f(0,1)
{{< /latex >}}

Okay, but how do we get \\(y(t)\\) for _any_ \\(t\\)? All we have is a point
and a derivative at that point. To get any further, we have to extrapolate.
At our point, the value of \\(y\\) is changing at a rate of \\(y'(0)\\). It's
not unreasonable to guess, then, that after some distance \\(h\\), the
value of \\(y\\) would be:
{{< latex >}}
y(0+h) \approx y(0) + hf(0,y(0))
{{< /latex >}}
Importantly, \\(h\\) shouldn't be large, since the further away you
extrapolate, the less accorate your approximation becomes. If we need
to get further than \\(0+h\\), we can use the new point we just computed,
and extrapolate again:
{{< latex >}}
y(0+2h) \approx y(0+h) + hf(0+h,y(0+h))
{{< /latex >}}
And so on and so forth. The last thing I'm going to do is clean up the above
equations. First, let's stop assuming that our initial \\(t=0\\); that doesn't work
in the general case. Instead, we'll use \\(t_0\\) to denote the initial value of
\\(t\\). Similarly, we can use \\(y_0\\) to denote our (known _a priori_)
value of \\(y(t_0)\\). From then, observing that we are moving in discrete
steps, estimating one point after another. We can denote the first point
we estimate as \\((t_1, y_1)\\), the second as \\((t_2, y_2)\\), and so on.
Since we're always moving our \\(t\\) by adding \\(h\\) to it, the formula for
a general \\(t_n\\) is easy enough:
{{< latex >}}
t_n = t_0 + nh
{{< /latex >}}
We can get the formula for \\(y_n\\) by looking at our aproximation steps above.
{{< latex >}}
\begin{gather}
  y_1 = y_0 + hf(t_0,y_0) \\
  y_2 = y_1 + hf(t_1,y_1) \\
  ... \\
  y_{n+1} = y_n + hf(t_n,y_n)
\end{gather}
{{< /latex >}}

What we have arrived at is [Euler's method](https://mathworld.wolfram.com/EulerForwardMethod.html).
It works pretty well, especially considering how easy it is to implement.
We could write it in Chapel like so:

```Chapel
proc runEulerMethod(step: real, count: int, t0: real, y0: real) {
    var y = y0;
    var t = t0;
    for i in 1..count {
        y += step*f(t,y);
        t += step;
    }
    return y;
}
```

Just eight lines of code!

# Linear Multistep Methods
In Euler's method, we use only the derivative at the _current_ point to drive
our subsequent approximations. Why does that have to be the case, though?
After a few steps, we've computed a whole history of points, each of which
gives us information about \\(y\\)'s derivatives. Of course, the more recent
derivative are more useful to us (they're closer, so using them would reduce
the error we get from extrapolating). Nevertheless, we can still try to
mix in a little bit of "older" information into our estimate of the derivative,
and try to get a better result that way. Take, for example, the two-step
Adams-Bashforth method:

{{< latex >}}
y_{n+2} = y_{n+1} + h\left[\frac{3}{2}f(t_{n+1}, y_{n+1})-\frac{1}{2}f(t_n,y_n)\right]
{{< /latex >}}

In this method, the _two_ last known derivatives are used for computing the next
point. More generally, an \\(s\\)-step linear multistep method uses the last
\\(s\\) points to compute the next, and always takes a linear combination
of the derivatives and the $y_n$ values themselves. The general form is:

{{< latex >}}
\sum_{j=0}^s a_jy_{n+j} = h\sum_{j=0}^s b_jf(t_{n+j}, y_{n+j})
{{< /latex >}}

Let's focus for a second only on methods where the very last \\(y\\) value
is used to compute the next. Furthermore, we won't allow our specifications
to depend on the slope at the point we're trying to compute (i.e., when computing
\\(y_{n+s}\\), we do not look at \\(f(t_{n+s}, y_{n+s})\\). The methods
we will consider will consequently have the following form:

{{< latex >}}
y_{n+s} = y_{n+s-1} + h\sum_{j=0}^{s-1} b_jf(t_{n+j}, y_{n+j})
{{< /latex >}}

Importantly, notice that a method of the above form is pretty
much uniquely determined by its list of coefficients. The
number of steps \\(s\\) is equal to the number of coefficients
\\(b_j\\), and rest of the equation simply "zips together"
a coefficient with a fixed expression representing the derivative
at a particular point in the past.

What we would like to do now is to be able to represent this
information entirely at compile-time. That is, we would like
to be able to fully describe a method using some Chapel code,
feed it to a function that implements the method, and pay
little-to-nothing in terms of runtime overhead. By encoding
information about a numerical method at compile-time,
we can also ensure that the function is correctly invoked. For
instance, Euler's method needs a single initial coordinate to
get started, \\(y_{0}\\), but the two-steps Adams-Bashforth method
above requires both \\(y_{0}\\) and another coordinate \\(y_{1}\\).
If we were to describe methods using a runtime construct, such
as perhaps a list, it would be entirely possible for a user to
feed in a list with not enough points, or too many. We can be
smarter than that.

Okay, but how do we represent a list at compile-time? Chapel
does have `param`s, which are integer, boolean, or string
values that can occur at the type-level and are known when
a program is being compiled. However, `param`s only support
primitive values; a list is far from a primitive value!

What is necessary is a trick used quite frequently in
Haskell land, a place where I suppose I am now a permanent
resident if not already a citizen. We can create two
completely unrelated records, just to tell the compiler
about two new types:

```Chapel
record empty {}
record cons {
    param head: real;
    type tail;
}
```

The new `cons` type is generic. It needs two
things: a `real` for the `head` field, and
a `type` for the `tail` field. Coming up
with a `real` is not that hard; we can, for instance,
write:

```Chapel
type t1 = cons(1.0, ?);
```

Alright, but what goes in the place of the question mark?
The second parameter of `cons` is a type, and a fully-instantiated
`cons` is itself a type. Thus, perhaps we could try nesting
two `cons`s:

```Chapel
type t2 = cons(1.0, cons(?, ?));
```

Once again, the first parameter to `cons` is a real number, and
coming up with one of those is easy.

```Chapel
type t3 = cons(1.0, cons(2.0, ?));
```

I hope it's obvious enough that we can keep
repeating this process to build up a longer
and longer chain of numbers. However, we need to
do something different if we are to get out of this
`cons`-loop. That's where `empty` comes in.

```
type t4 = cons(1.0, cons(2.0, empty));
```

Now, `t4` is a fully-concrete type. Note that we
haven't yet -- nor will we -- actually allocate
any records of type `cons` or `empty`; the purpose
of these two is solely to live at the type level.
Specifically, if we encode our special types
of linear multistep methods as lists of coefficients,
we can embed these lists of coefficients into Chapel
types by using `cons` and `empty`. Check out the
table below.

|Method | List | Chapel Type Expression |
|-------|----------------------|-------------------|
| Degenerate case | \\(\\varnothing\\) | `empty` |
| Euler's method  | \\(1\\)                  | `cons(1.0, empty)` |
| Two-step Adams-Bashforth  | \\(\frac{3}{2}, -\\frac{1}{2}\\)                  | `cons(3.0/2.0, cons(-0.5, empty))` |

As I mentioned earlier, the number of steps \\(s\\) a method takes
(which is also the number of initial data points required to
get the method started) is the number of coefficients. Thus,
we need a way to take a type like `t4` and count how many coefficients
it contains. Fortunately, Chapel allows us to write functions on
types, as well as to pattern match on types' arguments. An
important caveat is that matching occurs when a function is called,
so we can't somehow write `myType == cons(?h, ?t)`. Instead, we need
to write a function that _definitely_ accepts a type `cons` with
some arguments:

{{< codelines "chapel" "linear-multistep/lmm.chpl" 9 9 >}}

The above function works roughly as follows:
```Chapel
// initial call
initial(cons(1.0, cons(2.0, empty)))
// ?w gets matched with 1.0
// ?t gets matched with cons(2.0, empty)
// the function return is evaluated
1 + initial(t)
// Since t was matched with cons(2.0, empty), this is the same as
1 + initial(cons(2.0, empty))
// This is another call to initial.
// ?w gets matched with 2.0
// ?t gets matched with empty
// evaluation proceeds like last time
1 + 1 + initial(empty)
```
Here, we could we get stuck; what's `initial(empty)`? We need
another overload for `initial` that handles this case.

{{< codelines "chapel" "linear-multistep/lmm.chpl" 8 8 >}}

With this, the final expression will effectively become `1 + 1 + 0 = 2`,
which is indeed the length of the given list.

Another thing we're going to need is the ability to get a coefficient at a particular
location in the list. Such a function is pretty simple, though it
does also use pattern matching via `?w` and `?t`:

{{< codelines "chapel" "linear-multistep/lmm.chpl" 10 16 >}}

Alright, now's the time to get working on the actual implementation
of linear multistep methods. A function implementing them will
need a few parameters:

* A `type` parameter representing the linear multi-step method
  being implemented. This method will be encoded using our
  `cons` lists.
* A `real` representing the step size \\(h\\).
* An `int` representing how many iterations we should run.
* A `real` representing the starting x-coordinate of the method
  \\(t_0\\). We don't need to feed in a whole list \\(t_0, \\ldots, t_{s-1}\\)
  since each \\(t_{i}\\) is within a fixed distance \\(h\\) of
  its neighbors.
* As many `real` coordinates as required by the method,
  \\(y_0, \\ldots, y_{s-1}\\). 

Here's the corresponding Chapel code:

{{< codelines "chapel" "linear-multistep/lmm.chpl" 18 19 >}}

That last little bit, `n : real ... initial(method)`, is an
instance of Chapel's variable arguments functionality. Specifically,
this says that `runMethod` takes multiple `real` arguments,
the exact number of which is given by the expression `initial(method)`.
We've just seen that `initial` computes the number of coefficients
in a list, so the function will expect exactly as many \\(y_i\\) 
y-coordinates as there are coefficients.

The actual code is pretty boring.
{{< codelines "chapel" "linear-multistep/lmm.chpl" 18 39 >}}

That completes our implementation. All that's left is to actually
use our API! First, we need to encode the various methods:

{{< codelines "chapel" "linear-multistep/lmm.chpl" 46 47 >}}

Finally, let's try this on a function. Wikipedia runs
all the methods with the differential equation \\(y' = y\\),
for which the corresponding function `f` is:

{{< codelines "chapel" "linear-multistep/lmm.chpl" 42 42 >}}

Then, we can run the methods as follows (with initial coordinate
\\((0, 1)\\), as Wikipedia does).

{{< codelines "chapel" "linear-multistep/lmm.chpl" 50 56 >}}