solve pyyawt compatibility

.gitignore

@@ -2,3 +2,5 @@
 **/*.pyc
 **/*.pyo
 *.egg-info/*
+build/
+

build/lib/rsHRF/VERSION

@@ -0,0 +1 @@
+1.5.8

build/lib/rsHRF/__about__.py

@@ -0,0 +1,12 @@
+"""Base module variables."""
+from os import path as op
+
+with open(op.join(op.dirname(op.realpath(__file__)), "VERSION"), "r") as fh:
+    __version__ = fh.read().strip('\n')
+
+__packagename__ = 'rsHRF'
+__url__ = 'https://github.com/BIDSapps/rsHRF'
+
+DOWNLOAD_URL = (
+    'https://github.com/BIDS-Apps/{name}/'.format(
+        name=__packagename__))

build/lib/rsHRF/__init__.py

@@ -0,0 +1,10 @@
+import rsHRF.spm_dep
+import rsHRF.processing
+import rsHRF.canon
+import rsHRF.sFIR
+import rsHRF.parameters
+import rsHRF.fourD_rsHRF
+import rsHRF.CLI
+__all__ = ["spm_dep", "processing", "canon", "utils",
+           "sFIR", "parameters",   "basis_functions",
+           "fourD_rsHRF.py", "CLI.py", "iterative_wiener_deconv"]

build/lib/rsHRF/__main__.py

@@ -0,0 +1,3 @@
+from rsHRF import CLI
+
+CLI.main()

build/lib/rsHRF/basis_functions/__init__.py

@@ -0,0 +1 @@
+from . import basis_functions

build/lib/rsHRF/basis_functions/basis_functions.py

@@ -0,0 +1,69 @@
+import math
+import numpy as np
+from scipy.stats import gamma
+from rsHRF import canon, spm_dep, sFIR
+
+import warnings
+warnings.filterwarnings("ignore")
+
+"""
+BASIS FUNCTION COMPUTATION
+"""
+
+def get_basis_function(bold_sig_shape, para):
+    N, nvar = bold_sig_shape
+    dt = para['dt'] # 'time bin for basis functions {secs}';
+    l  = para['len']
+    if "gamma" in para['estimation']:
+        pst = np.arange(0,l+0.01,dt) #the 0.01 is because the rounding is different between python and matlab
+        bf  = gamma_bf(pst,para['order'])
+    elif 'fourier' in para['estimation'] or 'hanning' in para['estimation']:
+        pst = np.arange(0,l+0.01,dt) #the 0.01 is because the rounding is different between python and matlab
+        pst = pst/max(pst)
+        bf  = fourier_bf(pst,para)
+    elif 'canon' in para['estimation']:
+        bf = canon2dd_bf(bold_sig_shape, para)
+    bf = spm_dep.spm.spm_orth(np.asarray(bf))
+    return bf
+
+def canon2dd_bf(data_shape, para):
+    """
+    Returns canon basis functions
+    """
+    N, nvar = data_shape
+    bf = canon.canon_hrf2dd.wgr_spm_get_canonhrf(para)
+    if 'Volterra' in para:
+        if para['Volterra'] == 2:
+            bf2 = np.einsum('i...,j...->ij...',bf.T,bf.T).reshape(-1, bf.shape[0]).T
+            bf  = np.column_stack((bf, bf2))
+    return bf
+
+def fourier_bf(pst, para):
+    """
+    Returns Fourier (Hanning) basis functions
+    """
+    if "hanning" in para['estimation']:
+        g = (1 - np.cos(2*math.pi*pst)) / 2
+    else:
+        g = np.ones(len(pst))
+    sin_ = lambda x : np.sin(x*2*math.pi)
+    cos_ = lambda x : np.cos(x*2*math.pi) 
+    sin_ = np.vectorize(sin_)
+    cos_ = np.vectorize(cos_)
+    arr = np.arange(1, para['order']+1)
+    s = sin_(np.einsum('i,j->ij',arr,pst))
+    c = cos_(np.einsum('i,j->ij',arr,pst))
+    s = np.multiply(g,s)
+    c = np.multiply(g,c)
+    g = np.expand_dims(g, axis=1).T
+    return np.concatenate((g, s, c), axis=0).T
+
+def gamma_bf(u,h):
+    """
+    Returns Gamma basis functions
+    """
+    arr = np.arange(2, h+2)
+    f = np.vectorize(gamma.pdf, signature='(n),()->(n)')
+    m = np.power(2, arr)
+    return f(u,m).T
+

build/lib/rsHRF/canon/__init__.py

@@ -0,0 +1 @@
+from . import canon_hrf2dd

build/lib/rsHRF/canon/canon_hrf2dd.py

@@ -0,0 +1,29 @@
+import numpy as np
+from ..spm_dep    import spm
+
+import warnings
+warnings.filterwarnings("ignore")
+
+def wgr_spm_get_canonhrf(xBF):
+    dt     = xBF['dt']
+    fMRI_T = xBF['T']
+    p      = np.array([6, 16, 1, 1, 6, 0, 32], dtype=float)
+    p[len(p) - 1] = xBF['len']
+    bf = spm.spm_hrf(dt, p, fMRI_T)
+    bf = bf[:, np.newaxis]
+    # time-derivative
+    if xBF['TD_DD']:
+        dp   = 1
+        p[5] = p[5] + dp
+        D    = (bf[:, 0] - spm.spm_hrf(dt, p, fMRI_T)) / dp
+        D    = D[:, np.newaxis]
+        bf   = np.append(bf, D, axis=1)
+        p[5] = p[5] - dp
+        # dispersion-derivative
+        if xBF['TD_DD'] == 2:
+            dp   = 0.01
+            p[2] = p[2] + dp
+            D    = (bf[:, 0] - spm.spm_hrf(dt, p, fMRI_T)) / dp
+            D    = D[:, np.newaxis]
+            bf   = np.append(bf, D, axis=1)
+    return bf

build/lib/rsHRF/fourD_rsHRF.py

@@ -0,0 +1,201 @@
+import os
+import matplotlib
+matplotlib.use('agg')
+import numpy             as np
+import nibabel           as nib
+import scipy.io          as sio
+import matplotlib.pyplot as plt
+from bids.layout  import BIDSLayout, parse_file_entities 
+from scipy        import stats, signal
+from scipy.sparse import lil_matrix
+from rsHRF        import spm_dep, processing, parameters, basis_functions, utils, iterative_wiener_deconv
+
+import warnings
+warnings.filterwarnings("ignore")
+
+def demo_rsHRF(input_file, mask_file, output_dir, para, p_jobs, file_type=".nii", mode="bids", wiener=False, temporal_mask=[]):
+    # book-keeping w.r.t parameter values
+    if 'localK' not in para or para['localK'] == None:
+        if para['TR']<=2:
+            para['localK'] = 1
+        else:
+            para['localK'] = 2
+    # creating the output-directory if not already present
+    if not os.path.isdir(output_dir):
+        os.mkdir(output_dir)
+    # for four-dimensional input
+    if mode != 'time-series':
+        if mode == 'bids' or mode == 'bids w/ atlas':
+            name = input_file.split('/')[-1].split('.')[0]
+            v1 = spm_dep.spm.spm_vol(input_file)
+        else:
+            name = input_file.split('/')[-1].split('.')[0]   
+            v1 = spm_dep.spm.spm_vol(input_file)
+        if mask_file != None:
+            if mode == 'bids':
+                mask_name = mask_file.split('/')[-1].split('.')[0]
+                v = spm_dep.spm.spm_vol(mask_file)
+            else:
+                mask_name = mask_file.split('/')[-1].split('.')[0]
+                v = spm_dep.spm.spm_vol(mask_file)
+            if file_type == ".nii" or file_type == ".nii.gz":
+                brain = spm_dep.spm.spm_read_vols(v)
+            else:
+                brain = v.agg_data().flatten(order='F')
+            if  ((file_type == ".nii" or file_type == ".nii.gz") and \
+                    v1.header.get_data_shape()[:-1] != v.header.get_data_shape()) or \
+                ((file_type == ".gii" or file_type == ".gii.gz") and \
+                    v1.agg_data().shape[0]!= v.agg_data().shape[0]):
+                raise ValueError ('Inconsistency in input-mask dimensions' + '\n\tinput_file == ' + name + file_type + '\n\tmask_file == ' + mask_name + file_type)
+            else:
+                if file_type == ".nii" or file_type == ".nii.gz" :
+                    data   = v1.get_fdata()      
+                else:
+                    data   = v1.agg_data()
+        else:
+            print('No atlas provided! Generating mask file...')
+            if file_type == ".nii" or file_type == ".nii.gz" :
+                data   = v1.get_fdata() 
+                brain  = np.nanvar(data.reshape(-1, data.shape[3]), -1, ddof=0)
+            else:
+                data   = v1.agg_data()
+                brain  = np.nanvar(data, -1, ddof=0) 
+            print('Done')
+        voxel_ind  = np.where(brain > 0)[0]
+        mask_shape = data.shape[:-1]
+        nobs       = data.shape[-1]
+        data1      = np.reshape(data, (-1, nobs), order='F').T
+        bold_sig = stats.zscore(data1[:, voxel_ind], ddof=1)
+   # for time-series input
+    else:
+        name = input_file.split('/')[-1].split('.')[0]
+        data1 = (np.loadtxt(input_file, delimiter=","))
+        if data1.ndim == 1:
+            data1 = np.expand_dims(data1, axis=1)
+        nobs = data1.shape[0]
+        bold_sig = stats.zscore(data1, ddof=1)
+    if len(temporal_mask) > 0 and len(temporal_mask) != nobs:
+            raise ValueError ('Inconsistency in temporal_mask dimensions.\n' + 'Size of mask: ' + str(len(temporal_mask)) + '\n' + 'Size of time-series: ' + str(nobs))
+    bold_sig = np.nan_to_num(bold_sig)
+    bold_sig_deconv = processing. \
+                      rest_filter. \
+                      rest_IdealFilter(bold_sig, para['TR'], para['passband_deconvolve'])   
+    bold_sig = processing. \
+               rest_filter. \
+               rest_IdealFilter(bold_sig, para['TR'], para['passband'])   
+    data_deconv  = np.zeros(bold_sig.shape)
+    event_number = np.zeros((1, bold_sig.shape[1]))
+    print('Retrieving HRF ...')
+    #Estimate HRF for the fourier / hanning / gamma / cannon basis functions
+    if not (para['estimation'] == 'sFIR' or para['estimation'] == 'FIR'):
+        bf = basis_functions.basis_functions.get_basis_function(bold_sig.shape, para)
+        beta_hrf, event_bold = utils.hrf_estimation.compute_hrf(bold_sig, para, temporal_mask, p_jobs, bf=bf)
+        hrfa = np.dot(bf, beta_hrf[np.arange(0, bf.shape[1]), :])
+    #Estimate HRF for FIR and sFIR
+    else:
+        para['T'] = 1        
+        beta_hrf, event_bold = utils.hrf_estimation.compute_hrf(bold_sig, para, temporal_mask, p_jobs)
+        hrfa = beta_hrf[:-1,:]
+    nvar = hrfa.shape[1]
+    PARA = np.zeros((3, nvar))
+    for voxel_id in range(nvar):
+        hrf1 = hrfa[:, voxel_id]
+        PARA[:, voxel_id] = \
+            parameters.wgr_get_parameters(hrf1, para['TR'] / para['T'])
+    print('Done')
+    print('Deconvolving HRF ...')
+    if para['T'] > 1:
+        hrfa_TR = signal.resample_poly(hrfa, 1, para['T'])
+    else:
+        hrfa_TR = hrfa
+    for voxel_id in range(nvar):
+        hrf = hrfa_TR[:, voxel_id]
+        if not wiener:
+            H = np.fft.fft(
+                np.append(hrf,
+                          np.zeros((nobs - max(hrf.shape), 1))), axis=0)
+            M = np.fft.fft(bold_sig_deconv[:, voxel_id])
+            data_deconv[:, voxel_id] = \
+                np.fft.ifft(H.conj() * M / (H * H.conj() + .1*np.mean((H * H.conj()))))
+        else:
+            data_deconv[:, voxel_id] = iterative_wiener_deconv.rsHRF_iterative_wiener_deconv(bold_sig_deconv[:, voxel_id], hrf)
+        event_number[:, voxel_id] = np.amax(event_bold[voxel_id].shape)
+    print('Done')
+    print('Saving Output ...')
+    # setting the output-path
+    if mode == 'bids' or mode == 'bids w/ atlas':
+        layout_output = BIDSLayout(output_dir)
+        entities = parse_file_entities(input_file)
+        sub_save_dir = layout_output.build_path(entities).rsplit('/',1)[0]            
+    else:
+        sub_save_dir = output_dir
+    if not os.path.isdir(sub_save_dir):
+        os.makedirs(sub_save_dir, exist_ok=True)
+    dic = {'para': para, 'hrfa': hrfa, 'event_bold': event_bold, 'PARA': PARA}
+    ext = '_hrf.mat'
+    if mode == "time-series":
+        dic["event_number"] = event_number
+        dic["data_deconv"]  = data_deconv
+        ext = '_hrf_deconv.mat'
+    name = name.rsplit('_bold', 1)[0]   
+    sio.savemat(os.path.join(sub_save_dir, name + ext), dic)
+    HRF_para_str = ['height', 'T2P', 'FWHM']
+    if mode != "time-series":
+        mask_data = np.zeros(mask_shape).flatten(order='F')
+        for i in range(3):
+            fname = os.path.join(sub_save_dir,
+                                 name + '_' + HRF_para_str[i])
+            mask_data[voxel_ind] = PARA[i, :]
+            mask_data = mask_data.reshape(mask_shape, order='F')
+            spm_dep.spm.spm_write_vol(v1, mask_data, fname, file_type)
+            mask_data = mask_data.flatten(order='F')
+        fname = os.path.join(sub_save_dir, name + '_eventnumber')
+        mask_data[voxel_ind] = event_number
+        mask_data = mask_data.reshape(mask_shape, order='F')
+        spm_dep.spm.spm_write_vol(v1, mask_data, fname, file_type)
+        mask_data = np.zeros(data.shape)
+        dat3 = np.zeros(data.shape[:-1]).flatten(order='F')
+        for i in range(nobs):
+            fname = os.path.join(sub_save_dir, name + '_deconv')
+            dat3[voxel_ind] = data_deconv[i, :]
+            dat3 = dat3.reshape(data.shape[:-1], order='F')
+            if file_type == ".nii" or file_type == ".nii.gz" :
+                mask_data[:, :, :, i] = dat3
+            else:
+                mask_data[:, i] = dat3
+            dat3 = dat3.flatten(order='F')
+        spm_dep.spm.spm_write_vol(v1, mask_data, fname, file_type)
+    pos = 0
+    while pos < hrfa_TR.shape[1]:
+        if np.any(hrfa_TR[:,pos]):
+            break 
+        pos += 1
+    event_plot = lil_matrix((1, nobs))
+    if event_bold.size:
+        event_plot[:, event_bold[pos]] = 1
+    else:
+        print("No Events Detected!")
+        return 0
+    event_plot = np.ravel(event_plot.toarray())
+    plt.figure()
+    plt.plot(para['TR'] * np.arange(1, np.amax(hrfa_TR[:, pos].shape) + 1),
+             hrfa_TR[:, pos], linewidth=1)
+    plt.xlabel('time (s)')
+    plt.savefig(os.path.join(sub_save_dir, name + '_hrf_plot.png'))
+    plt.figure()
+    plt.plot(para['TR'] * np.arange(1, nobs + 1),
+             np.nan_to_num(stats.zscore(bold_sig[:, pos], ddof=1)),
+             linewidth=1)
+    plt.plot(para['TR'] * np.arange(1, nobs + 1),
+             np.nan_to_num(stats.zscore(data_deconv[:, pos], ddof=1)),
+             color='r', linewidth=1)
+    markerline, stemlines, baseline = \
+        plt.stem(para['TR'] * np.arange(1, nobs + 1), event_plot)
+    plt.setp(baseline, 'color', 'k', 'markersize', 1)
+    plt.setp(stemlines, 'color', 'k')
+    plt.setp(markerline, 'color', 'k', 'markersize', 3, 'marker', 'd')
+    plt.legend(['BOLD', 'Deconvolved BOLD', 'Events'], loc='best')
+    plt.xlabel('time (s)')
+    plt.savefig(os.path.join(sub_save_dir, name + '_deconvolution_plot.png'))   
+    print('Done')
+    return 0

build/lib/rsHRF/iterative_wiener_deconv.py

@@ -0,0 +1,39 @@
+# import pyyawt
+import pywt
+import numpy as np 
+from rsHRF.processing import knee 
+
+def rsHRF_iterative_wiener_deconv(y, h, Iterations=1000):
+    N = y.shape[0]
+    nh = max(h.shape)
+    h = np.append(h, np.zeros((N - nh, 1)))
+    H = np.fft.fft(h, axis=0)
+    Y = np.fft.fft(y, axis=0)
+    coeffs = pywt.wavedec(abs(y), 'db2', level=1)
+    sigma = np.median(np.abs(coeffs[-1])) / 0.6745
+    Phh = np.square(abs(H))
+    sqrdtempnorm = ((((np.linalg.norm(y-np.mean(y), 2)**2) - (N-1)*(sigma**2))) / (np.linalg.norm(h,1))**2)
+    Nf = (sigma**2)*N
+    tempreg = Nf/sqrdtempnorm
+    Pxx0 = np.square(abs(np.multiply(Y, (np.divide(np.conj(H), (Phh + N*tempreg))))))
+    Pxx = Pxx0
+    Sf = Pxx.reshape(-1, 1)
+    for i in range(0, Iterations):
+        M = np.divide(np.multiply(np.multiply(np.conjugate(H), Pxx), Y), np.add(np.multiply(np.square(abs(H)), Pxx), Nf))
+        PxxY = np.divide(np.multiply(Pxx, Nf), np.add(np.multiply(np.square(abs(H)), Pxx), Nf))
+        Pxx = np.add(PxxY, np.square(abs(M)))
+        Sf = np.concatenate((Sf, Pxx.reshape(-1, 1)), axis=1)
+    dSf = np.diff(Sf, 1, 1)
+    dSfmse = np.mean(np.square(dSf), axis=1)
+    _, idx = knee.knee_pt(dSfmse)
+    idm = np.argmin(dSfmse)
+    ratio = np.abs(dSfmse[idx] - dSfmse[idm])/(np.abs(np.max(dSfmse) - np.min(dSfmse)))
+    if ratio > 0.5: 
+        id0 = idm 
+    else:
+        id0 = idx 
+    # Safe indexing to avoid IndexError
+    safe_index = min(id0 + 1, Sf.shape[1] - 1)
+    Pxx = Sf[:, safe_index]
+    WienerFilterEst = np.divide(np.multiply(np.conj(H), Pxx), np.add(np.multiply(np.square(abs(H)), Pxx), Nf))
+    return np.real(np.fft.ifft(np.multiply(WienerFilterEst, Y)))