Note, the following is partly to document my explorations in computing lat/on grids in numpy lately, in case others come across these issues in the future. It is partly a plea for some development of numerically accurate functions for computing lat/lon grids from a combination of inputs: bounds, counts, and resolutions.

I have been playing around with the decimal package a bit lately, and I discovered the concept of "fused multiply-add" operations for improved accuracy. I have come to realize that fma operations could be used to greatly improve the accuracy of linspace() and arange(). In particular, I have been needing improved results for computing latitude/longitude grids, which tend to be done in float32's to save memory (at least, this is true in data I come across).

Since fma is not available yet in python, please consider the following C snippet:

```
$ cat test_fma.c
#include<stdlib.h>
#include<math.h>
#include<stdio.h>

int main(){
    float res = 0.01;
    float x0 = -115.0;
    int cnt = 7001;
    float x1 = res * cnt + x0;
    float x2 = fma(res, cnt, x0);
    printf("x1 %.16f  x2 %.16f\n", x1, x2);
    float err1 = fabs(x1 - -44.990000000000000);
    float err2 = fabs(x2 - -44.990000000000000);
    printf("err1 %.5e  err2 %.5e\n", err1, err2);
    return 0;
}
$ gcc test_fma.c -lm -o test_fma

$ ./test_fma
x1 -44.9899978637695312  x2 -44.9900016784667969
err1 2.13623e-06  err2 1.67847e-06

```

And if you do similar for double-precision, the fma still yields significant accuracy over double precision explicit multiply-add.

Now, to the crux of my problem. It is next to impossible to generate a non-trivial numpy array of coordinates, even in double precision, without hitting significant numerical errors. Which has lead me down the path of using the decimal package (which doesn't play very nicely with numpy because of the lack of casting rules for it). Consider the following:

```

$ cat test_fma.py
from __future__ import print_function
import numpy as np
res = np.float32(0.01)
cnt = 7001
x0 = np.float32(-115.0)
x1 = res * cnt + x0
print("res * cnt + x0 = %.16f" % x1)
x = np.arange(-115.0, -44.99 + (res / 2), 0.01, dtype='float32')
print("len(arange()): %d  arange()[-1]: %16f" % (len(x), x[-1]))
x = np.linspace(-115.0, -44.99, cnt, dtype='float32')
print("linspace()[-1]: %.16f" % x[-1])

$ python test_fma.py
res * cnt + x0 = -44.9900015648454428
len(arange()): 7002  arange()[-1]:       -44.975044
linspace()[-1]: -44.9900016784667969

```

arange just produces silly results (puts out an extra element... adding half of the resolution is typically mentioned as a solution on mailing lists to get around arange()'s limitations -- I personally don't do this).

linspace() is a tad bit better than arange(), but still a little bit worse than just computing the value straight out. It also produces a result that is on par with fma(), so that is reassuring that it is about as accurate as can be at the moment. This also matches up almost exactly to what I would get if I used Decimal()s and then casted to float32's for final storage, so that's not too bad as well.

So, does it make any sense to improve arange by utilizing fma() under the hood? Also, any plans for making fma() available as a ufunc?

Notice that most of my examples required knowing the number of grid points ahead of time. But what if I didn't know that? What if I just have the bounds and the resolution? Then arange() is the natural fit, but as I showed, its accuracy is lacking, and you have to do some sort of hack to do a closed interval.

Thoughts? Comments?

Ben Root

(donning flame suit)

_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion