Skip to content
Browse files

ENH: Some preliminary work on an error function.

  • Loading branch information...
1 parent d364a6c commit 6b40ad7de51e97adcf4f8b15aef8c950dae86410 @arokem committed with Sep 26, 2011
Showing with 75 additions and 9 deletions.
  1. +35 −6 nitime/fmri/prf.py
  2. +40 −3 nitime/fmri/tests/test_prf.py
View
41 nitime/fmri/prf.py
@@ -232,13 +232,17 @@ def stat_non_lin(x, kwargs): return x
if non_lin_params is None:
# Just make up a dummy dict:
- non_lin_params = dict(dummy=None)
+ non_lin_params = dict(kwargs=None)
+ # Reshape these, before performing the cross-multiplication:
+ # For stimulus, time is the last dimension and remains:
stim_re = np.matrix(stim.reshape(np.prod(stim.shape[:-1]),stim.shape[-1]))
+ # For the PRF, both dims are space, so we ravel:
p_re = np.matrix(p.ravel())
+ # Cross-multiply the matrices and apply the static non-linearity:
neural_response = stat_non_lin(np.array(stim_re.T * p_re.T).squeeze(),
- non_lin_params)
+ **non_lin_params)
# If no h function/vector is provided, the output is just the neural response:
if h is None:
@@ -247,7 +251,7 @@ def stat_non_lin(x, kwargs): return x
# If it is a function, we need to generate the vector of values for the
# convolution:
if callable(h):
- h = h(h_dur, h_params)
+ h = h(h_dur, **h_params)
# Now do the convolution with the result of that operation:
bold_response = np.convolve(neural_response, h)
@@ -383,13 +387,38 @@ def projection_matrix(A, norm=None):
return np.eye(A.shape[0]) - A * ols_matrix(A, norm)
-def exponent(x, kwargs):
+def exponent(x, n=None):
"""
A simple static non-linearity for :func:`response`
"""
- if kwargs is None:
+ if n is None:
return x
else:
- return x ** kwargs['n']
+ return x ** n
+
+
+def err_func(bold, stim, prf, hrf, prf_params=None, hrf_params=None,
+ extra_params=None, non_lin=None, non_lin_params=None):
+ """
+ This is the error-function to optimize on.
+
+ Given BOLD data from a voxel, prf starting params and hrf starting params,
+ find the best fit params for these factors.
+
+ bold, stim are given, but the rest should be fit.
+
+ """
+
+ p = prf(n_x=stim.shape[0], n_y=stim.shape[1])
+
+ if callable(hrf):
+ h = hrf(33)
+
+ r = response(stim, p, h=h, h_dur=33,
+ h_params=hrf_params, stat_non_lin=non_lin,
+ non_lin_params=non_lin_params)
+
+ # This is the error:
+ return r - bold
View
43 nitime/fmri/tests/test_prf.py
@@ -2,6 +2,7 @@
import numpy.testing as npt
import nitime.fmri.prf as prf
+import nitime.fmri.hrf as hrf
def test_norms():
"""
@@ -74,8 +75,44 @@ def test_response():
stim = np.random.randn(n,m,t)
p = prf.gabor(n,m)
- # The simplest call, no static non-linearity, no
- prf.response(stim,p)
+ # The simplest call, no static non-linearity, no HRF:
+ r1 = prf.response(stim,p)
# Static non-linearity, but no hrf:
- prf.response(stim,p,stat_non_lin=prf.exponent,non_lin_params=dict(n=0.5))
+ r2 = prf.response(stim,p,stat_non_lin=prf.exponent,non_lin_params=dict(n=0.5))
+
+ # HRF function, but no static non-linearity:
+ r3 = prf.response(stim,p,h=hrf.gamma, h_dur=33, h_params=dict(tau=1.5, Fs=0.5))
+
+ # HRF vector and no static non-linearity:
+ r4 = prf.response(stim,p,h=hrf.gamma(33,tau=1.5,Fs=0.5))
+
+ # These two last ones should be equal:
+ npt.assert_equal(r3,r4)
+
+ # The full monty:
+ r5 = prf.response(stim,p,
+ h=hrf.gamma, h_dur=33, h_params=dict(tau=1.5, Fs=0.5),
+ stat_non_lin=prf.exponent, non_lin_params=dict(n=0.5))
+
+def test_errfunc():
+ """
+ Testing that the error-function does what it's supposed to.
+
+ """
+ stim_dur = 180
+ stim = np.random.randn(100,100,stim_dur)
+ real_p = prf.gaussian(100,100)
+ real_h = hrf.gamma(33)
+ bold = prf.response(stim,real_p,real_h)
+
+ e = prf.err_func(bold, stim, prf.gaussian, hrf.gamma)
+
+ # In this case, the prf and hrf are right on, so there should be no error:
+ npt.assert_equal(np.all(e==0),True)
+
+ # Now move the PRF slightly away:
+ real_p = prf.gaussian(100,100,x0=60)
+ bold = prf.response(stim,real_p,real_h)
+
+ e = prf.err_func(bold, stim, prf.gaussian, hrf.gamma)

0 comments on commit 6b40ad7

Please sign in to comment.
Something went wrong with that request. Please try again.