From 2fc9a7b0a3933b411fc5f5ff5eac1ed69ab9896e Mon Sep 17 00:00:00 2001
From: Michael Waskom <mwaskom@stanford.edu>
Date: Thu, 16 Jan 2014 15:45:02 -0800
Subject: [PATCH] Allow for precomputed highpass filter matrix

---
 moss/glm.py            | 12 ++++++++++--
 moss/tests/test_glm.py | 14 ++++++++++++++
 2 files changed, 24 insertions(+), 2 deletions(-)

diff --git a/moss/glm.py b/moss/glm.py
index c56a5d1..0426320 100644
--- a/moss/glm.py
+++ b/moss/glm.py
@@ -157,7 +157,7 @@ class DesignMatrix(object):
     """
     def __init__(self, design, hrf_model, ntp, regressors=None, confounds=None,
                  artifacts=None, condition_names=None, confound_pca=False,
-                 tr=2, hpf_cutoff=128, oversampling=16):
+                 tr=2, hpf_cutoff=128, hpf_kernel=None, oversampling=16):
         """Initialize the design matrix object.
 
         Parameters
@@ -197,6 +197,9 @@ def __init__(self, design, hrf_model, ntp, regressors=None, confounds=None,
             sampling interval (in seconds) of the data/design
         hpf_cutoff : float
             filter cutoff (in seconds), or None to skip the filter
+        hpf_kernel : ntp x ntp array
+            precomputed matrix to implement the highpass filter; overrides the
+            hpf_cutoff (i.e. it is not checked if they match).
         oversampling : float
             construction of the condition evs and convolution
             are performed on high-resolution data with this oversampling
@@ -232,6 +235,8 @@ def __init__(self, design, hrf_model, ntp, regressors=None, confounds=None,
         conditions = self._subsample_condition_matrix()
         conditions -= conditions.mean()
         pp_heights = (conditions.max() - conditions.min()).tolist()
+        if hpf_kernel is not None:
+            hpf_cutoff = hpf_kernel
         if hpf_cutoff is not None:
             conditions = self._highpass_filter(conditions, hpf_cutoff)
 
@@ -402,7 +407,10 @@ def _validate_component(self, comp, name_base):
 
     def _highpass_filter(self, mat, cutoff):
         """Highpass-filter each column in mat."""
-        F = fsl_highpass_matrix(self._ntp, cutoff, self.tr)
+        if np.isscalar(cutoff):
+            F = fsl_highpass_matrix(self._ntp, cutoff, self.tr)
+        else:
+            F = cutoff
         for key, vals in mat.iteritems():
             mat[key] = np.dot(F, vals)
         return mat
diff --git a/moss/tests/test_glm.py b/moss/tests/test_glm.py
index 91838dd..aa206ff 100644
--- a/moss/tests/test_glm.py
+++ b/moss/tests/test_glm.py
@@ -308,6 +308,20 @@ def test_design_matrix_confound_pca():
     nt.assert_equal(n_confounds, good_dims)
 
 
+def test_design_matrix_precomputed_kernel():
+    """Test that we can use a precomputed highpass filter matrix."""
+    F = glm.fsl_highpass_matrix(20, 8, 2)
+    hrf = glm.GammaDifferenceHRF()
+
+    design = pd.DataFrame(dict(condition=["one", "two"],
+                               onset=[5, 20]))
+    regressors = rs.randn(20, 2)
+
+    X_1 = glm.DesignMatrix(design, hrf, 20, regressors=regressors, hpf_cutoff=8)
+    X_2 = glm.DesignMatrix(design, hrf, 20, regressors=regressors, hpf_kernel=F)
+
+    npt.assert_array_equal(X_1.design_matrix, X_2.design_matrix)
+
 def test_highpass_matrix_shape():
     """Test the filter matrix is the right shape."""
     for n_tp in 10, 100: