pennmem / ptsa

PTSA - Python Time Series Analysis
https://pennmem.github.io/ptsa/
GNU General Public License v3.0
20 stars 9 forks source link

MorletWaveletFilter timeseries dtype bug #203

Closed esolomon closed 6 years ago

esolomon commented 6 years ago

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.

mivade commented 6 years ago

Workaround for the time being:

eeg = eeg.to_ptsa().astype(np.double)
mivade commented 6 years ago

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.

mivade commented 6 years ago

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.