yasa: Spindles detection failing for some users

Wow, this was really some awesome detective work @remrama

Setting nopython=False does not help, but I think it makes sense because the compilation isn’t failing and so that isn’t changing behavior in our case. I think that only changes the fallback behavior if the initial compilation fails, which would’ve raised earlier errors if it happened. What we want to do is force the use of the fallback object mode. To do that, we can add forceobj=True, and that makes it work.

But the problem is, there are 100 other things that also fix the problem, and they all fix the problem independently. I’ve been racking my brain to try and figure out what common underlying change is going on in Numba when any of these fixes are used:

Things that fix _corr()

• Set forceobj=True - I didn’t find this too useful once I realized there were other fixes. It’s basically saying it works when treated as normal Python code, but in these other patches it doesn’t have to be.
• Set fastmath=True - This felt like it was going to be informative because it changes the underlying computations, so this would’ve made total sense. But alas, you can revert back to default fastmath=False and all the next things will will fix it.
• Set signature to "float64(float64[::1], float64[::1])" or just remove it entirely - We’re using 1d arrays that should comply with both C- and F-contiguous expectations (check with the ndarray.flags), so I think the current signature specification of A-contiguous is supposed to work. And it does on most systems. But not here, though you can fix the OG problem by telling numba to expect C-contiguous, or don’t tell it anything at all and it will infer it as C.
• Set error_model="numpy" - This also seemed promising. The error_model parameter changes how numba handles zero-division. But the default is python, which should raise an error when any zero-division attempt is made, but it isn’t doing that even when I make a simplified example. Despite not having errors in test cases, when switching to error_model="numpy", OG problem is fixed.
• Set boundscheck=True - Works with just this change.
• Use np.arange instead of range - Yes. This change – all alone by itself – fixes the problem.
• Pass xm2s, ym2s, r_num as keyword arguments - Do this instead of setting them within the function and it works. Between this and the last one, I thought maybe this was a sneaky object-mode thing, but numba tutorials are constantly setting variables and using range in their functions.

Again, all of these are patches by themselves. I have a jupyter notebook with printouts of each separate fix. I didn’t link it here bc it’s a bit extra but I could if it would help. While it seems great (easy fix!), it’s a bit unsettling to me because I have no idea what the root cause is. I’m hesitant to implement any of these because to me it all seems like completely unexpected behavior and has the potential for more downstream disruptions. I feel like these are just a ton of red herrings. I also feel like I’m losing my mind.

@remrama really confusing indeed. Btw I realized that we can even simplify further by using zip:

for xi, yi in zip(x, y):
    xm = xi - mx
    ym = yi - my
    r_num += xm * ym
    xm2s += xm**2
    ym2s += ym**2

It seems to be even a bit faster than the previous implementation. Can you check whether you have the bug when using zip?

Thanks so much @remrama for your in-depth investigation! So I think the error was introduced in #86. Interestingly, someone opened a very similar PR on antropy recently, and their solution was not to do an if..else statement, but instead to add a very small value to the denominator to avoid division by zero error.

Can you please try to run the code with:

import sys
eps = sys.float_info.epsilon

@jit("float64(float64[:], float64[:])", nopython=True)
def _corr(x, y):
    """Fast Pearson correlation."""
    n = x.size
    mx, my = x.mean(), y.mean()
    xm2s, ym2s, r_num = 0, 0, 0
    for i in range(n):
        xm = x[i] - mx
        ym = y[i] - my
        r_num += xm * ym
        xm2s += xm**2
        ym2s += ym**2
    r_d1 = np.sqrt(xm2s)
    r_d2 = np.sqrt(ym2s)
    r_den = r_d1 * r_d2
    return r_num / (r_den + eps)

If this is not working you can just do:

return r_num / (r_den + 10e-9)

Ok I think I found the problem. I was able to reproduce it on my dusty Desktop. It has to do with the moving correlation thresholding, and I know this sounds insane but I think it depends on a tiny difference in Windows version.

Before I go into detail can some people help verify something for me? Please show your output from this slightly modified version of the original sample code.

import sys
import platform
import numpy as np
import pandas as pd
import mne
import numba
import scipy
import yasa
import antropy

print("System", platform.platform())
print("Python", sys.version)
print(f"YASA {yasa.__version__}")
print(f"MNE {mne.__version__}")
print(f"NumPy {np.__version__}")
print(f"Pandas {pd.__version__}")
print(f"Numba {numba.__version__}")
print(f"SciPy {scipy.__version__}")
print(f"Antropy {antropy.__version__}")

sf = 200
data = np.loadtxt('data_N2_spindles_15sec_200Hz.txt')

# I expect this to work for everyone.
print("=" * 10, "Run with no moving correlation threshold.", "=" * 10)
sp = yasa.spindles_detect(data, sf, thresh={"corr": None}, verbose=True)

# I expect this to still break (no spindles found) for a subset of users.
print("=" * 10, "Run with moving correlation threshold set to default 0.65.", "=" * 10)
sp = yasa.spindles_detect(data, sf, verbose=True)

I wonder if someone who is getting the error could mess around with some of the detection parameters to see if that helps identify where the detection is failing. Sorry I can’t try it myself, I’m just not able to reproduce the initial problem.