<commit_before>"""
This module provides definitions of various hemodynamic response functions (hrf).

In particular, it provides <NAME>'s canonical HRF, AFNI's default HRF, and
a spectral HRF.
"""
import numpy as N
import numpy.linalg as L

from neuroimaging import traits

from neuroimaging.modalities.fmri import filters
from neuroimaging.modalities.fmri.utils import LinearInterpolant as interpolant

def glover2GammaDENS(peak_hrf, fwhm_hrf):
    alpha = N.power(peak_hrf / fwhm_hrf, 2) * 8 * N.log(2.0)
    beta = N.power(fwhm_hrf, 2) / peak_hrf / 8 / N.log(2.0)
    coef = peak_hrf**(-alpha) * N.exp(peak_hrf / beta)
    return coef, filters.GammaDENS(alpha + 1., 1. / beta)

def _glover(peak_hrf=(5.4, 10.8), fwhm_hrf=(5.2, 7.35), dip=0.35):
    coef1, gamma1 = glover2GammaDENS(peak_hrf[0], fwhm_hrf[0])
    coef2, gamma2 = glover2GammaDENS(peak_hrf[1], fwhm_hrf[1])
    f = filters.GammaCOMB([[coef1,gamma1],[-dip*coef2, gamma2]])
    dt = 0.02
    t = N.arange(0, 50 + dt, dt)
    c = (f(t) * dt).sum()
    return filters.GammaCOMB([[coef1/c,gamma1],[-dip*coef2/c, gamma2]])

# Glover, 'canonical HRF'

glover = filters.Filter(_glover(), names=['glover'])
glover_deriv = filters.Filter([_glover(), _glover().deriv(const=-1.)],
                              names=['glover', 'dglover'])
canonical = glover

# AFNI's default HRF (at least at some point in the past)

afni = filters.Filter(filters.GammaDENS(9.6, 1.0/0.547), ['gamma'])

class SpectralHRF(filters.Filter):
    '''
    Delay filter with spectral or Taylor series decomposition
    for estimating delays.

    Liao et al. (2002).
    '''

    dt = traits.Float(0.02)
    tmax = traits.Float(500.0)
    ncomp = traits.Int(2)
    names = traits.ListStr(['glover'])

    spectral = traits.true

    def __init__(self, input_hrf=canonical, **keywords):
        filters.Filter.__init__(self, input_hrf, **keywords)
        if self.n != 1:
            raise ValueError, 'expecting one HRF for spectral decomposition'
        self.deltaPCA()

    def deltaPCA(self, tmax=50., lower=-15.0, delta=N.arange(-4.5,4.6,0.1)):
        """
        Perform an expansion of fn, shifted over the values in delta.
        Effectively, a Taylor series approximation to fn(t+delta), in delta,
        with basis given by the filter elements. If fn is None, it assumes
        fn=IRF[0], that is the first filter.

        >>> import numpy as N
        >>> from pylab import plot, title, show
        >>> from neuroimaging.modalities.fmri.hrf import glover, glover_deriv, SpectralHRF
        >>>
        >>> ddelta = 0.25
        >>> delta = N.arange(-4.5,4.5+ddelta, ddelta)
        >>> time = N.arange(0,20,0.2)
        >>>
        >>> hrf = glover_deriv
        >>>
        >>> canonical = glover
        >>> taylor = hrf.deltaPCA(delta)
        >>> curplot = plot(time, taylor.components[1](time))
        >>> curplot = plot(time, taylor.components[0](time))
        >>> curtitle=title('Shift using Taylor series -- components')
        >>> show()
        >>>
        >>> curplot = plot(delta, taylor.coef[1](delta))<commit_msg>small trait cleanup<commit_after>"""
This module provides definitions of various hemodynamic response functions (hrf).

In particular, it provides <NAME>'s canonical HRF, AFNI's default HRF, and
a spectral HRF.
"""
import numpy as N
import numpy.linalg as L

from neuroimaging import traits

from neuroimaging.modalities.fmri import filters
from neuroimaging.modalities.fmri.utils import LinearInterpolant as interpolant

def glover2GammaDENS(peak_hrf, fwhm_hrf):
    alpha = N.power(peak_hrf / fwhm_hrf, 2) * 8 * N.log(2.0)
    beta = N.power(fwhm_hrf, 2) / peak_hrf / 8 / N.log(2.0)
    coef = peak_hrf**(-alpha) * N.exp(peak_hrf / beta)
    return coef, filters.GammaDENS(alpha + 1., 1. / beta)

def _glover(peak_hrf=(5.4, 10.8), fwhm_hrf=(5.2, 7.35), dip=0.35):
    coef1, gamma1 = glover2GammaDENS(peak_hrf[0], fwhm_hrf[0])
    coef2, gamma2 = glover2GammaDENS(peak_hrf[1], fwhm_hrf[1])
    f = filters.GammaCOMB([[coef1,gamma1],[-dip*coef2, gamma2]])
    dt = 0.02
    t = N.arange(0, 50 + dt, dt)
    c = (f(t) * dt).sum()
    return filters.GammaCOMB([[coef1/c,gamma1],[-dip*coef2/c, gamma2]])

# Glover, 'canonical HRF'

glover = filters.Filter(_glover(), names=['glover'])
glover_deriv = filters.Filter([_glover(), _glover().deriv(const=-1.)],
                              names=['glover', 'dglover'])
canonical = glover

# AFNI's default HRF (at least at some point in the past)

afni = filters.Filter(filters.GammaDENS(9.6, 1.0/0.547), ['gamma'])

class SpectralHRF(filters.Filter):
    '''
    Delay filter with spectral or Taylor series decomposition
    for estimating delays.

    Liao et al. (2002).
    '''

    dt = traits.Float(0.02)
    tmax = traits.Float(500.0)
    ncomp = traits.Int(2)
    names = traits.ListStr(['glover'])

    def __init__(self, input_hrf=canonical, spectral=True, **keywords):
        filters.Filter.__init__(self, input_hrf, **keywords)
        self.spectral = spectral
        if self.n != 1:
            raise ValueError, 'expecting one HRF for spectral decomposition'
        self.deltaPCA()

    def deltaPCA(self, tmax=50., lower=-15.0, delta=N.arange(-4.5,4.6,0.1)):
        """
        Perform an expansion of fn, shifted over the values in delta.
        Effectively, a Taylor series approximation to fn(t+delta), in delta,
        with basis given by the filter elements. If fn is None, it assumes
        fn=IRF[0], that is the first filter.

        >>> import numpy as N
        >>> from pylab import plot, title, show
        >>> from neuroimaging.modalities.fmri.hrf import glover, glover_deriv, SpectralHRF
        >>>
        >>> ddelta = 0.25
        >>> delta = N.arange(-4.5,4.5+ddelta, ddelta)
        >>> time = N.arange(0,20,0.2)
        >>>
        >>> hrf = glover_deriv
        >>>
        >>> canonical = glover
        >>> taylor = hrf.deltaPCA(delta)
        >>> curplot = plot(time, taylor.components[1](time))
        >>> curplot = plot(time, taylor.components[0](time))
        >>> curtitle=title('Shift using Taylor series -- components')
        >>> show()
        >>>
        >>> curplot = plot(delta, taylor.coef[1](delta))
{'fraction_non_alphanumeric': 0.14581312641448432, 'fraction_numerical': 0.03136113805366958, 'mean_word_length': 2.9584133077415227, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 113, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2351800', 'n_tokens_mistral': 2303, 'n_tokens_neox': 2136, 'n_words': 530}