issue_comments

On second thought we should add this to a FAQ page.

Yes that sounds right. Thanks!

• https://github.com/pydata/xarray/pull/5560 introduced "use_bottleneck" option to disable/enable using bottleneck. can we close this issue or keep it open?

Oh yes, of course! I've underestimated the low precision of float32 values above 2**24. Thanks for the hint.

I think (!) xarray is not effected any longer, but pandas is. Bisecting the GIT history leads to commit 0b9ab2d1, which means that xarray >= v0.10.9 should not be affected. Uninstalling bottleneck is also a valid workaround.

<s>Bottleneck's documentation explicitly mentions that no error is raised in case of an overflow. But it seams to be very evil behavior, so it might be worth reporting upstream.</s> What do you think? (I think kwgoodman/bottleneck#164 is something different, isn't it?) Edit: this is not an overflow. It's a numerical error by not applying pairwise summation.

A couple of minimal examples:

```python

import numpy as np import pandas as pd import xarray as xr import bottleneck as bn bn.nanmean(np.ones(225, dtype=np.float32))
0.5 pd.Series(np.ones(225, dtype=np.float32)).mean()
0.5 xr.DataArray(np.ones(2**25, dtype=np.float32)).mean() # not affected for this version <xarray.DataArray ()> array(1., dtype=float32) ```

Done with the following versions: bash $ pip3 freeze Bottleneck==1.2.1 numpy==1.16.1 pandas==0.24.1 xarray==0.11.3 ...

The difference is that Bottleneck does the sum in the naive way, whereas NumPy uses the more numerically stable pairwise summation.

Oh hm, I think I didn't really understand what happens in bottleneck.nanmean(). I understand that integers can overflow and that float32 have varying absolute precision. The max float32 3.4E+38 is not hit here. So how can the mean of a list of ones be 0.5?

Isn't this what bottleneck is doing? Summing up a bunch of float32 values and then dividing by the length?

```

d = np.ones(2**25, dtype=np.float32) d.sum()/np.float32(len(d)) 1.0 ```

Ah ok, I suppose bottleneck is indeed now avoided for float32 xarray. Yeah that issue is for a different function, but the source of the problem and proposed solution in the thread is the same - use higher precision intermediates for float32 (double arithmetic); a small speed vs accuracy/precision trade off.

Is it worth changing bottleneck to use double for single precision reductions? AFAICT this is a matter of changing npy_DTYPE0 to double in the float{64,32} versions of functions in reduce_template.c.

Would it be worth adding a warning (until the right solution is found) if someone is doing .mean() on a DataArray which is float32?

I would rather pick option (1) above, that is, "Stop using bottleneck on float32 arrays"

Sorry to unearth this issue again, but I just got bitten by this quite badly. I'm looking at absolute temperature perturbations and bottleneck's implementation together with my data being loaded as float32 (correctly, as it's stored like that) causes an error on the size of the perturbations I'm looking for.

Example:

``` In [1]: import numpy as np ...: import bottleneck

In [2]: a = 300np.ones((800*2,), dtype=np.float32)

In [3]: np.mean(a) Out[3]: 300.0

In [4]: bottleneck.nanmean(a) Out[4]: 302.6018981933594 ```

Would it be worth adding a warning (until the right solution is found) if someone is doing .mean() on a DataArray which is float32?

Based a little experimentation (https://gist.github.com/leifdenby/8e874d3440a1ac96f96465a418f158ab) bottleneck's mean function builds up significant errors even with moderately sized arrays if they are float32, so I'm going to stop using .mean() as-is from now and always pass in dtype=np.float64.

@matteodefelice you didn't decide on float32, but your data is stored that way. It's really hard to make choices about numerical precision for computations automatically: if we converted automatically to float64, somebody else would be complaining about unexpected memory usage :).

Looking at our options, we could:

1. Stop using bottleneck on float32 arrays, or provide a flag or option to disable using bottleneck. This is not ideal, because bottleneck is much faster.
2. Automatically convert float32 arrays to float64 before doing aggregations. This is not ideal, because it could significant increase memory requirements.
3. Add a dtype option for aggregations (like NumPy) and consider defaulting to dype=np.float64 when doing aggregations on float32 arrays. I would generally be happy with this, but bottleneck currently doesn't provide the option currently.
4. Write a higher precision algorithm for bottleneck's mean.

Thanks all guys for the replies. @Aegaeon I get the same your results with bottleneck... @shoyer The point is that I haven't decided the use of float32 and — yes — using .astype(np.float64) solves the issue...the point is that is not an expected behaviour, with such standard dataset I would not expect any problem related to numerical precision...

Yes, this is probably related to the fact that .mean() in xarray uses bottleneck if available, and bottleneck has a slightly different mean implementation, quite possibly with a less numerically stable algorithm.

The fact that the dtype is float32 is a sign that this is probably a numerical precision issue. Try casting with .astype(np.float64) and see if the problem goes away.

If you really cared about performance using float32, the other thing to do to improve conditioning is to subtract and add a number close to the mean, e.g., (ds.var167 - 270).mean() + 270.

I think this might be a problem with bottleneck? My interpretation of _create_nan_agg_method in xarray/core/ops.py is that it may use bottleneck to get the mean unless you pass skipna=False or specify multiple axes. And,

```python In [2]: import bottleneck In [3]: bottleneck.version Out[3]: '1.2.0'

In [6]: bottleneck.nanmean(ds.var167.data) Out[6]: 261.6441345214844 ```

Forgive me if I'm wrong, I'm still a bit new.

Does it make a difference if you load the data first? (ds.var167.load().mean()) Or use python 3?

Also on macOS, and I can reproduce.

Using python 2.7.11, xarray 0.9.1, dask 0.14.1 installed through Anaconda. I get the same results with xarray 0.9.1-38-gc0178b7 from GitHub.

```python In [3]: ds = xarray.open_dataset('ERAIN-t2m-1983-2012.seasmean.nc')

In [4]: ds.var167.mean() Out[4]: <xarray.DataArray 'var167' ()> array(261.6441345214844, dtype=float32) ```

Curiously, I get the right results with skipna=False...

python In [10]: ds.var167.mean(skipna=False) Out[10]: <xarray.DataArray 'var167' ()> array(278.6246643066406, dtype=float32)

... or by specifying coordinates to average over:

python In [5]: ds.var167.mean(('time', 'lat', 'lon')) Out[5]: <xarray.DataArray 'var167' ()> array(278.6246643066406, dtype=float32)

Ok, I am on MacOS: - Python 2.7.13 from Macports - Dask 0.14.1 from Macports - xarray from GitHub

I can't reproduce this:

```python In [6]: ds = xr.open_dataset('./Downloads/ERAIN-t2m-1983-2012.seasmean.nc')

In [7]: ds.var167.mean() Out[7]: <xarray.DataArray 'var167' ()> array(278.6246643066406, dtype=float32)

In [8]: ds.var167.data.mean() Out[8]: 278.62466 ```

which version of xarray, dask, python are you using?