rfftn axes specification no longer works

[Jeitan]

I ran into this issue with numpy 1.15.4 when some existing code (that was very robust and working fine) started to give crazy results, and I isolated it to np.fft.rfftn. I'm actually hoping that I'm just doing something stupid, but I can't find an error in my logic so far.

Since it's non-trivial to create on-the-fly example data for this, I've made some test data available in a npy file here (simulated ocean waves). The following code illustrates the issue - I have a cube of data which represents a simulated wave field in 3D (Time, X, Y). I am trying to make a PSD. Because the input data are real, I'm using rfftn. So, I want to perform the full complex FFT along axes 1 and 2 and the conjugate symmetric FFT along axis 0.

Reproducing code example:

import numpy as np
import numpy.fft as nfft
import matplotlib.pyplot as plt

with open('waves.npy', 'rb') as fid:
        waves = np.load(fid)

goodn = nfft.rfftn(waves.T)
good1 = nfft.rfft(waves, axis=0)
goodn1 = nfft.rfftn(waves, axes=[0])
goodn3 = nfft.rfftn(waves.T, axes=[0,1,2])

fig, axs = plt.subplots(figsize=(8,8), ncols=2, nrows=2)
axs[0][0].imshow(np.ma.log10(np.abs(goodn[0,:,:])**2), aspect='auto')
axs[1][0].imshow(np.ma.log10(np.abs(good1[0,:,:])**2), aspect='auto')
axs[1][1].imshow(np.ma.log10(np.abs(goodn1[0,:,:])**2), aspect='auto')
axs[0][1].imshow(np.ma.log10(np.abs(goodn3[0,:,:])**2), aspect='auto')

All of the above actually produces plots that are correct. In the top two plots in the image below, you can see the deep water dispersion curve.

If I transpose my data so the 'real' axis that I want to use is the last axis and use rfftn out of the box, it works - "goodn", upper left. If I do the same thing but specify all 3 axes (without changing the order), it also works - "goodn3", upper right. If I use rfft and specify axis=0, instead of using rfftn, I get junk because I haven't taken the transform of the X and Y dimensions - "good1", lower left. If I do essentially the equivalent by using rfftn but only specifying axes=[0], I get the same exact junk - "goodn1", lower right.

Error message:

However, what I had been doing, and what I want to do, and from the docs I think I should be able to do, is this:

realgood = nfft.rfftn(waves, axes=[1,2,0])

In an older version of numpy (at least up to 1.13), that worked perfectly. In numpy 1.15.4, that gave me the kind of junk described above. It looked like the X and Y axes were not actually being transformed (full complex transform) like they should be. I thought I might get around it by upgrading to 1.16.2 ... however now what that does is give me this error:

Traceback (most recent call last):

  File "<ipython-input-11-3dad548cb9e8>", line 1, in <module>
    runfile('C:/Users/asmith/Documents/Git/python_utils/scripts/gist.py', wdir='C:/Users/asmith/Documents/Git/python_utils/scripts')

  File "C:\Users\asmith\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 786, in runfile
    execfile(filename, namespace)

  File "C:\Users\asmith\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 110, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/asmith/Documents/Git/python_utils/scripts/gist.py", line 57, in <module>
    realgood = nfft.rfftn(waves, axes=[1,2,0])

  File "C:\Users\asmith\AppData\Local\Continuum\anaconda3\lib\site-packages\mkl_fft\_numpy_fft.py", line 1043, in rfftn
    output = mkl_fft.rfftn_numpy(a, s, axes)

  File "mkl_fft\_pydfti.pyx", line 947, in mkl_fft._pydfti.rfftn_numpy

  File "mkl_fft\_pydfti.pyx", line 846, in mkl_fft._pydfti._fftnd_impl

  File "mkl_fft\_pydfti.pyx", line 724, in mkl_fft._pydfti._iter_fftnd

  File "mkl_fft\_pydfti.pyx", line 667, in mkl_fft._pydfti._init_nd_shape_and_axes

ValueError: all axes must be unique

If I switch the two spatial axes (which should be the same, yes? Axes 1 and 2 are the same size), I get a different error:

realgood = nfft.rfftn(waves, axes=[2,1,0])

gives me this:

Traceback (most recent call last):

  File "<ipython-input-15-3dad548cb9e8>", line 1, in <module>
    runfile('C:/Users/asmith/Documents/Git/python_utils/scripts/gist.py', wdir='C:/Users/asmith/Documents/Git/python_utils/scripts')

  File "C:\Users\asmith\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 786, in runfile
    execfile(filename, namespace)

  File "C:\Users\asmith\AppData\Local\Continuum\anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 110, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/asmith/Documents/Git/python_utils/scripts/gist.py", line 57, in <module>
    realgood = nfft.rfftn(waves, axes=[2,1,0])

  File "C:\Users\asmith\AppData\Local\Continuum\anaconda3\lib\site-packages\mkl_fft\_numpy_fft.py", line 1043, in rfftn
    output = mkl_fft.rfftn_numpy(a, s, axes)

  File "mkl_fft\_pydfti.pyx", line 951, in mkl_fft._pydfti.rfftn_numpy

ValueError: could not broadcast input array from shape (512,64) into shape (512,512)

Numpy/Python version information:

1.16.2 3.7.3 (default, Mar 27 2019, 17:13:21) [MSC v.1915 64 bit (AMD64)]

[charris]

Did you get numpy from anaconda? I ask, because the mkl_fft code is not ours.

[Jeitan]

I did. I didn't realize that getting numpy through anaconda was somehow different? Should I try getting it through pip instead?

Although actually even so ... shouldn't numpy work when installed with anaconda?

[charris]

Mixing pip and anaconda can be tricky and cause problems, there are several issues related to that here. I don't know what the current status is, but I think it is not recommended. Might open an issue at https://github.com/ContinuumIO/anaconda-issues.

Although actually even so ... shouldn't numpy work when installed with anaconda?

Different libraries are linked.

[tylerjereddy]

I'd suggest perhaps trying conda install nomkl and then reinstalling the package if conda doesn't adjust it for you to remove the mkl patching.

See here: https://docs.anaconda.com/mkl-optimizations/#uninstalling-mkl

[Jeitan]

Hmmm, I might need to investigate ... I did "conda update conda" when I updated numpy, so that thinks it's up to date, but "conda install nomkl" can't find anything called nomkl. I will see what I can straighten out there and update in a couple days.

[charris]

I assume that anaconda is running tests, so perhaps we need a test that exposes this problem. @Jeitan Can you try running the tests? Something like

In [2]: np.fft.test('full')
NumPy version 1.15.4
NumPy relaxed strides checking option: True
..............................                                           [100%]
30 passed in 0.92 seconds

[Jeitan]

@charris Okay I did that and here it is:

In [3]: np.fft.test('full')
NumPy version 1.16.2
NumPy relaxed strides checking option: True
..............................                                           [100%]
30 passed in 1.17 seconds

This was from the IPython console inside of Spyder (and I verified that the errors above still show in the current session).

[Jeitan]

I also learned that nomkl is not available on Windows (this is a Windows system I'm having issues on). I don't have an easy way to test anything on a Linux system ...

I did try to create a clean conda environment with the same results - I just installed python 3.7, conda, jupyter, numpy, and matplotlib. Python 3.7.3, conda 4.6.12, numpy 1.16.2, matplotlib 3.0.3, and various things associated with jupyter. The mkl that got installed was mkl 2019.3, mkl_fft 1.0.10, mkl_random 1.0.2.

[Jeitan]

It's looking more and more like some kind of package/library alignment. I tried making clean environments based on defaults for anaconda 5.0.1, 5.3, 2018.12, and 2019.03. The 5.0.1 worked fine, the 5.3 worked fine. As soon as it hit Anaconda 2018.12, the broadcasting error showed up. All were Python 3.x, except 5.3, which for some reason had conflicts so I caved and used Python 2.7. The versions included therein were numpy=1.15.1, mkl=2018.0.3, mkl_fft=1.0.4. That got me thinking and I made two non-full-anaconda environments, with the following requirements:

1. python=2 numpy=1.15.1 mkl=2018.0.3 mkl_fft=1.0.4 matplotlib
2. python=3 numpy=1.15.1 mkl=2018.0.3 mkl_fft=1.0.4 matplotlib

The same code worked in in environment 1 but hit the broadcasting error in environment 2. So I tried the most recent versions of everything, except using python 2, and voila it works fine. I'll open a ticket on the Anaconda issues page, thanks for all the pointers!

[charris]

@Jeitan Would it be correct to say that NumPy needs a test that fails for this particulat issue? Ideally this sort of problem should be caught early rather than after release.

[Jeitan]

That would probably be a good thing, I would say ... but right now I really have no idea what is actually wrong, so I wouldn't know how to test it. Although I suppose any test for rfftn that performs the conjugate symmetric transform over an axis that isn't the last one would have caught this, whatever it may be that causes it.