pennmem/ptsa

MorletWaveletFilter timeseries dtype bug

Closed this issue · 3 comments

The following code doesn't work with the latest version of cmlreaders (0.7.2) and ptsa (2.0.0):

import numpy as np
from cmlreaders import CMLReader, get_data_index
df = get_data_index("r1")
s = 'R1111M'
exp = 'FR1'
sess = 0
sessions = df[np.logical_and(df["subject"] == s, df['experiment']==exp)]['session'].unique()
pairs = CMLReader(s, exp).load("pairs")
reader = CMLReader(s, exp, sess)
evs = reader.load("events")
evs = evs[evs.type=='COUNTDOWN_START']
eeg = reader.load_eeg(events=evs, rel_start=0, rel_stop=10000, scheme=pairs)
eeg = eeg.to_ptsa()
#Extract spectral power
from ptsa.data.filters import MorletWaveletFilter
wf = MorletWaveletFilter(timeseries=eeg, freqs=np.array(range(70, 171, 10)), width=5, output='power')
powers = wf.filter()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-30-67a550edab27> in <module>()
     21 from ptsa.data.filters import MorletWaveletFilter
     22 wf = MorletWaveletFilter(timeseries=eeg, freqs=np.array(range(70, 171, 10)), width=5, output='power')
---> 23 powers = wf.filter()

~/anaconda3/envs/CML/lib/python3.6/site-packages/ptsa/data/filters/morlet.py in filter(self)
    116             mt.set_output_type(morlet.COMPLEX)
    117 
--> 118         mt.set_signal_array(timeseries_reshaped)
    119         mt.set_wavelet_pow_array(powers_reshaped)
    120         mt.set_wavelet_phase_array(phases_reshaped)

~/anaconda3/envs/CML/lib/python3.6/site-packages/ptsa/extensions/morlet/morlet.py in set_signal_array(self, signal_array)
    259 
    260     def set_signal_array(self, signal_array):
--> 261         return _morlet.MorletWaveletTransformMP_set_signal_array(self, signal_array)
    262 
    263     def set_wavelet_pow_array(self, wavelet_pow_array):

TypeError: Array of type 'double' required.  Array of type 'short' given

Looks like the timeseries dtype either needs to be output by cmlreaders to_ptsa() as float64 or reset to float64 by PTSA.

Workaround for the time being:

eeg = eeg.to_ptsa().astype(np.double)

https://github.com/pennmem/ptsa_new/blob/b7fb1cbeec28c2acf5f1c9cea06dedc70bf6df36/ptsa/data/filters/morlet.py#L48-L51

The easiest solution is to call .astype(np.double) here when setting the timeseries on the filter. Otherwise, the underlying C++ code would have to be reworked significantly to explicitly handle typecasting or use templates to provide for different types.

On second thought, it looks like the old PTSA EEG readers always converted to double, so maybe it would make more sense for to_ptsa to do this to avoid other issues.