Skip to content

Conversation

@mscheltienne
Copy link
Member

@mscheltienne mscheltienne commented Nov 20, 2022

Improvement to make sure the arguments n_fft, n_overlap and n_per_seg are integers in Welch's method, to avoid this crpytic traceback when a float is provided (which is common if you do duration * raw.info["sfreq"] and forget to convert to int):

Effective window size : 4.000 (s)
Traceback (most recent call last):
  File "/Users/rellis/Library/Application Support/JetBrains/IntelliJIdea2021.2/plugins/python/helpers/pydev/pydevd.py", line 1483, in _exec
    pydev_imports.execfile(file, globals, locals)  # execute the script
  File "/Users/rellis/Library/Application Support/JetBrains/IntelliJIdea2021.2/plugins/python/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "/Users/rellis/dev/git/[project]/analysis/power-spectrum-measures.py", line 152, in <module>
    epochPs = epochs.compute_psd(method='welch', picks=picks, fmin=FMIN, fmax=FMAX, n_fft=4*epochs.info["sfreq"])
  File "<decorator-gen-277>", line 12, in compute_psd
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/epochs.py", line 2031, in compute_psd
    return EpochsSpectrum(
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/spectrum.py", line 1042, in __init__
    self._compute_spectra(data, fmin, fmax, n_jobs, method_kw, verbose)
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/spectrum.py", line 380, in _compute_spectra
    result = self._psd_func(
  File "<decorator-gen-199>", line 12, in psd_array_welch
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 184, in psd_array_welch
    f_spect = parallel(my_spect_func(d, func=func, freq_sl=freq_sl,
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 184, in <genexpr>
    f_spect = parallel(my_spect_func(d, func=func, freq_sl=freq_sl,
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 73, in _spect_func
    spect = np.apply_along_axis(
  File "<__array_function__ internals>", line 5, in apply_along_axis
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/numpy/lib/shape_base.py", line 379, in apply_along_axis
    res = asanyarray(func1d(inarr_view[ind0], *args, **kwargs))
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/mne/time_frequency/psd.py", line 52, in _decomp_aggregate_mask
    _, _, spect = func(epoch)
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/spectral.py", line 737, in spectrogram
    window, nperseg = _triage_segments(window, nperseg,
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/spectral.py", line 1965, in _triage_segments
    win = get_window(window, nperseg)
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 2153, in get_window
    return winfunc(*params)
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 1095, in hamming
    return general_hamming(M, 0.54, sym)
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 1017, in general_hamming
    return general_cosine(M, [alpha, 1. - alpha], sym)
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/scipy/signal/windows/windows.py", line 113, in general_cosine
    fac = np.linspace(-np.pi, np.pi, M)
  File "<__array_function__ internals>", line 5, in linspace
  File "/Users/rellis/opt/anaconda3/envs/rob-nme/lib/python3.9/site-packages/numpy/core/function_base.py", line 120, in linspace
    num = operator.index(num)
TypeError: 'float' object cannot be interpreted as an integer
python-BaseException

Process finished with exit code 1

From the forum post https://mne.discourse.group/t/set-window-size-for-welch-via-compute-psd/5958/6

@mscheltienne mscheltienne changed the title nsure n_fft, n_overlap and n_per_seg are integers Ensure n_fft, n_overlap and n_per_seg are integers Nov 20, 2022
@mscheltienne mscheltienne changed the title Ensure n_fft, n_overlap and n_per_seg are integers MRG: Ensure n_fft, n_overlap and n_per_seg are integers Nov 20, 2022
Comment on lines 150 to 151
if n_per_seg is not None:
_validate_type(n_per_seg, "int", "n_per_seg")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if n_per_seg is not None:
_validate_type(n_per_seg, "int", "n_per_seg")
_validate_type(n_per_seg, ("int", None), "n_per_seg")

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, it does not support multiple types if one of them is "int"

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if types == "int":
_ensure_int(item, name=item_name, extra=extra)
return # terminate prematurely

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure? To me it’s mostly to speed up the check but it should still work

Copy link
Member Author

@mscheltienne mscheltienne Nov 20, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, the function fails with a KeyError for "int" a bit further, so I guess the early exit is also due to this

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That looks like correct/useful behavior

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

... @mscheltienne I think we are digressing a bit here. In MNE-Python we try to use _ensure_int when we want to ensure a user is passing an int-ike, e.g.:

$ git grep " = _ensure_int"
mne/_ola.py:        n_pts = _ensure_int(n_pts, 'n_pts')
mne/_ola.py:        n_samples = _ensure_int(n_samples, 'n_samples')
mne/_ola.py:        n_overlap = _ensure_int(n_overlap, 'n_overlap')
mne/_ola.py:        n_total = _ensure_int(n_total, 'n_total')
mne/bem.py:    threshold = _ensure_int(threshold, 'threshold')
mne/filter.py:    filter_length = _ensure_int(filter_length, 'filter_length')
mne/io/pick.py:        info = _ensure_int(info, 'info', 'an int or Info')
...

It sounds like that is your goal here, yes? If so, I think we should continue using _ensure_int here, and @mscheltienne if you are still opposed to this function / don't see the point / think it's useless, we should debate its merits in a separate issue that addresses that.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not opposed in the slightest, especially on a minor point like this one, since both functions fulfill the same role. I am however curious, am I missing something about _ensure_int? _validate_type calls it internally, and the casting seems to occur only from int to int (otherwise it raises).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the casting seems to occur only from int to int

Based on this phrasing perhaps what you miss is that more specifically "the casting occurs from any int-like that we want to support, to a native Python int".

Why is this useful? Well here is one problematic example case:

>>> _validate_type(np.uint32(1), 'int')
>>> np.uint32(1) - np.uint32(2)
<stdin>:1: RuntimeWarning: overflow encountered in uint_scalars
4294967295
>>> _ensure_int(np.uint32(1)) - np.uint32(2)
-1

This may seem contrived, but once users start working with / saving data to different formats (and we use I/O internally, too, internally), this sort of stuff can happen, and motivated implementing this function. So can this:

>>> type(np.array(1))
<class 'numpy.ndarray'>
>>> type(_ensure_int(np.array(1)))
<class 'int'>

I can't remember offhand why having this 0-d ndarray was problematic downstream, but I think it was at some point...

And even for less problematic cases (e.g., np.int32 vs np.int64 for small values), having different variants of int-like variables being passed can end up being inefficient. For example, if you write a numba function that consumes it, numba will now have to JIT-compile a new variant for every dtype (or raise if it doesn't know how to handle the inputs).

So anything that requires, is made safer, or is made more consistent by knowing that what you work with is a native Python int will benefit from this function. And given that it's hard to see ahead of time how potential failures will come up and account for them, it's useful to use it most times you need an int.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for the details, I'll think about it and about how to include this in my future codes. I did not realize this difference between ints could cause issues.

Co-authored-by: Alexandre Gramfort <alexandre.gramfort@m4x.org>
@mscheltienne mscheltienne changed the title MRG: Ensure n_fft, n_overlap and n_per_seg are integers Ensure n_fft, n_overlap and n_per_seg are integers Nov 20, 2022
Copy link
Member

@agramfort agramfort left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@larsoner merge if you think it is the expected fix

@larsoner larsoner merged commit e565817 into mne-tools:main Nov 21, 2022
@larsoner
Copy link
Member

Thanks @mscheltienne !

@mscheltienne mscheltienne deleted the types branch November 21, 2022 16:36
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants