Fixed issue with time-series convolution

skued/time_series/fitting.py

@@ -5,6 +5,7 @@
 
 import numpy as np
 from math import sqrt, log
+from functools import wraps
 
 
 def exponential(time, tzero, amp, tconst, offset=0):
@@ -43,7 +44,8 @@ def exponential(time, tzero, amp, tconst, offset=0):
     biexponential : bi-exponential curve with onset
     """
     return (
-        np.heaviside(time - tzero, 1 / 2) * (amp * (1 - np.exp(-(time - tzero) / tconst)))
+        np.heaviside(time - tzero, 1 / 2)
+        * (amp * (1 - np.exp(-(time - tzero) / tconst)))
         + offset
     )
 
@@ -93,36 +95,7 @@ def biexponential(time, tzero, amp1, amp2, tconst1, tconst2, offset=0):
     return arr
 
 
-def gauss_kernel(t, t0, fwhm):
-    """
-    Gaussian convolution kernel.
-
-    Parameters
-    ----------
-    x : array-like
-        Independant variable
-    t0 : array-like
-        t0 offset
-    fwhm : float
-        Full-width at half-maximum of the Gaussian kernel
-    """
-    std = fwhm / (2 * sqrt(2 * log(2)))
-    return (1/(np.sqrt(2*np.pi)*std)) * np.exp(-(1.0*t-t0)**2 / (2*std**2))
-
-
-def convolve(arr, kernel):
-    """ Convolution of array with kernel. """
-    #logger.debug("Convolving...")
-    npts = min(len(arr), len(kernel))
-    pad  = np.ones(npts)
-    tmp  = np.concatenate((pad*arr[0], arr, pad*arr[-1]))
-    norm = np.sum(kernel)
-    out  = np.convolve(tmp, kernel, mode='same')
-    noff = int((len(out) - npts)/2)
-    return out[noff:noff+npts]/norm
-
 # TODO: test with unevenly-spaced data points
-# TODO: figure out wth is going on with the width
 # TODO: add example
 # TODO: add example to user guide
 def with_irf(fwhm, f):
@@ -145,7 +118,42 @@ def with_irf(fwhm, f):
     --------
 
     """
+
+    @wraps(f)
     def f_(time, *args, **kwargs):
-        kernel = gauss_kernel(time, t0=0, fwhm=fwhm)
-        return convolve(f(time, *args, **kwargs), kernel)
-    return f_
+        kernel = _gauss_kernel(time, fwhm=fwhm)
+        return _convolve(f(time, *args, **kwargs), kernel)
+
+    return f_
+
+
+def _gauss_kernel(t, fwhm):
+    """
+    Gaussian convolution kernel.
+
+    Parameters
+    ----------
+    t : array-like
+        Independant variable
+    fwhm : float
+        Full-width at half-maximum of the Gaussian kernel
+    """
+    # It is important that the kernel be centered in the array `t`
+    # for the convolution operation to work properly
+    t0 = t[int(len(t) / 2)]
+
+    std = fwhm / (2 * sqrt(2 * log(2)))
+    return (1 / (np.sqrt(2 * np.pi) * std)) * np.exp(
+        -((1.0 * t - t0) ** 2) / (2 * std ** 2)
+    )
+
+
+def _convolve(arr, kernel):
+    """ Convolution of array with kernel. """
+    npts = min(len(arr), len(kernel))
+    pad = np.ones(npts)
+    tmp = np.concatenate((pad * arr[0], arr, pad * arr[-1]))
+    norm = np.sum(kernel)
+    out = np.convolve(tmp, kernel, mode="same")
+    noff = int((len(out) - npts) / 2)
+    return out[noff : noff + npts] / norm

tests/time_series/test_fitting.py

@@ -18,7 +18,9 @@ def setUp(self):
     def test_tzero_limits(self):
         """ Test that the output of ``exponential`` has the correct time-zero """
         t = np.arange(-10, 50, step=0.3)
-        I = exponential(t, tzero=self.tzero, amp=self.amp, tconst=self.tconst, offset=self.offset)
+        I = exponential(
+            t, tzero=self.tzero, amp=self.amp, tconst=self.tconst, offset=self.offset
+        )
 
         # Check that all values before time-zero are the amplitude
         self.assertTrue(np.all(np.equal(I[t < self.tzero], self.offset)))
@@ -27,15 +29,19 @@ def test_tzero_limits(self):
     def test_positivity(self):
         """ Test that the output of ``exponential`` is always positive. """
         t = np.arange(-10, 50, step=0.3)
-        I = exponential(t, tzero=self.tzero, amp=self.amp, tconst=self.tconst, offset=self.offset)
+        I = exponential(
+            t, tzero=self.tzero, amp=self.amp, tconst=self.tconst, offset=self.offset
+        )
 
         self.assertTrue(np.all(I >= 0))
 
     def test_amplitude(self):
         """ Test that the output of ``exponential`` is at most 0, 
         and at least ``amp`` (if ``amp`` is negative). """
         t = np.arange(-10, 50, step=0.3)
-        I = exponential(t, tzero=self.tzero, amp=self.amp, tconst=self.tconst, offset=self.offset)
+        I = exponential(
+            t, tzero=self.tzero, amp=self.amp, tconst=self.tconst, offset=self.offset
+        )
 
         self.assertTrue(np.all(np.less_equal(I, self.offset)))
         self.assertTrue(np.all(np.greater_equal(I, self.amp + self.offset)))
@@ -68,7 +74,7 @@ def test_tzero_limits(self):
             amp2=self.amp2,
             tconst1=self.tconst1,
             tconst2=self.tconst2,
-            offset=15
+            offset=15,
         )
 
         # Check that all values before time-zero are the amplitude
@@ -85,7 +91,7 @@ def test_positivity(self):
             amp2=self.amp2,
             tconst1=self.tconst1,
             tconst2=self.tconst2,
-            offset = abs(self.amp1) + abs(self.amp2)
+            offset=abs(self.amp1) + abs(self.amp2),
         )
 
         self.assertTrue(np.all(I >= 0))
@@ -100,7 +106,7 @@ def test_amplitude(self):
             amp2=self.amp2,
             tconst1=self.tconst1,
             tconst2=self.tconst2,
-            offset=10
+            offset=10,
         )
 
         self.assertTrue(np.all(np.less_equal(I, 10)))