# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:

import numpy as np
from scipy.stats import gamma

def biGammaHRF(times):
    """Return values for HRF at given times. Unit of time: second.
    This HRF is derived from the sum of two Gamma function, it starts at
    zero, and gets to zero sometime brfore 35 seconds.
    """
    # Gamma pdf for the peak
    peak_values = gamma.pdf(times, 6)
    # Gamma pdf for the undershoot
    undershoot_values = gamma.pdf(times, 12)
    # Combine them
    values = peak_values - 0.35 * undershoot_values
    # Scale max to 0.6
    return values / np.max(values) * 0.6


@@ -667,8 +667,405 @@
  * 0.6%0A%0A
+def vistaBiGammaHRF(times):%0A    %22%22%22Return values for HRF at given times. Unit of time: second.%0A    This HRF is derived from VISTA software.%0A    %22%22%22%0A    d1 = 5.4%0A    alpha1 = 5.98%0A    beta1 = 0.9%0A    c = 0.35%0A    d2 = 10.8%0A    alpha2 = 11.97%0A    beta2 = 0.9%0A    g = lambda t: ((t/d1)**alpha1)*np.exp((d1-t)/beta1)-c*((t/d2)**alpha2)*np.exp((d2-t)/beta2)%0A    return np.array(%5Bg(i) for i in times%5D)%0A%0A

add vistaHRF