Merge pull request #201 from matthew-brett/multi-run-eg-refactor

matthew-brett · matthew-brett · commit b6ce98b60bc7 · 2012-07-19T17:03:04.000-07:00
RF: Multi run example refactor

Large expansion of comments and some refactoring to multi-run example, with JB and I reviewing.
diff --git a/examples/formula/multi_session_contrast.py b/examples/formula/multi_session_contrast.py
@@ -1,85 +1,108 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
-""" This example is now out of date - please remind us to fix it
+""" Example of more than one run in the same model
 """
 import numpy as np
-import sympy
 
 from nipy.algorithms.statistics.api import Term, Formula, Factor
 from nipy.modalities.fmri import utils, hrf
 
-# hrf
+# hrf models we will use for each run.  Just to show it can be done, use a
+# different hrf model for each run
+h1 = hrf.glover
+h2 = hrf.afni
 
-h1 = sympy.Function('hrf1')
-h2 = sympy.Function('hrf2')
+# Symbol for time in general.  The 'events' function below will return models in
+# terms of 't', but we'll want models in terms of 't1' and 't2'.  We need 't'
+# here so we can substitute.
 t = Term('t')
 
-# Session 1
-
-t1 = Term('t1')
-c11 = utils.events([3,7,10], f=h1); c11 = c11.subs(t, t1)
-c21 = utils.events([1,3,9], f=h1); c21 = c21.subs(t, t1)
-c31 = utils.events([2,4,8], f=h1); c31 = c31.subs(t, t1)
-d1 = utils.fourier_basis([0.3,0.5,0.7]); d1 = d1.subs(t, t1)
-tval1 = np.linspace(0,20,101)
-
+# run 1
+t1 = Term('t1') # Time within run 1
+c11 = utils.events([3, 7, 10], f=h1) # Condition 1, run 1
+# The events utility returns a formula in terms of 't' - general time
+c11 = c11.subs(t, t1) # Now make it in terms of time in run 1
+# Same for conditions 2 and 3
+c21 = utils.events([1, 3, 9], f=h1); c21 = c21.subs(t, t1)
+c31 = utils.events([2, 4, 8], f=h1); c31 = c31.subs(t, t1)
+# Add also a fourier basis set for drift with frequencies 0.3, 0.5, 0.7
+d1 = utils.fourier_basis([0.3, 0.5, 0.7]); d1 = d1.subs(t, t1)
+
+# Here's our formula for run 1 signal terms of time in run 1 (t1)
 f1 = Formula([c11,c21,c31]) + d1
 
-# Session 2
-
-t2 = Term('t2')
-c12 = utils.events([3.3,7,10], f=h2); c12 = c12.subs(t, t2)
-c22 = utils.events([1,3.2,9], f=h2); c22 = c22.subs(t, t2)
-c32 = utils.events([2,4.2,8], f=h2); c32 = c32.subs(t, t2)
-d2 = utils.fourier_basis([0.3,0.5,0.7]); d2 = d2.subs(t, t2)
-tval2 = np.linspace(0,10,51)
-
-f2 = Formula([c12,c22,c32]) + d2
-
-sess_factor = Factor('sess', [1,2])
-
-# The multi session formula
-
-f = Formula([sess_factor.terms[0]]) * f1 + Formula([sess_factor.terms[1]]) * f2
-
-# Now, we evaluate the formula
-# the arrays now have 152=101+51 rows...
-
+# run 2
+t2 = Term('t2') # Time within run 2
+# Conditions 1 through 3 in run 2
+c12 = utils.events([3.3, 7, 10], f=h2); c12 = c12.subs(t, t2)
+c22 = utils.events([1, 3.2, 9], f=h2); c22 = c22.subs(t, t2)
+c32 = utils.events([2, 4.2, 8], f=h2); c32 = c32.subs(t, t2)
+d2 = utils.fourier_basis([0.3, 0.5, 0.7]); d2 = d2.subs(t, t2)
+
+# Formula for run 2 signal in terms of time in run 2 (t2)
+f2 = Formula([c12, c22, c32]) + d2
+
+# Factor giving constant for run. The [1, 2] means that there are two levels to
+# this factor, and that when we get to pass in values for this factor,
+# instantiating an actual design matrix (see below), a value of 1 means level
+# 1 and a value of 2 means level 2.
+run_factor = Factor('run', [1, 2])
+run_1_coder = run_factor.get_term(1) # Term coding for level 1
+run_2_coder = run_factor.get_term(2) # Term coding for level 2
+
+# The multi run formula will combine the indicator (dummy value) terms from the
+# run factor with the formulae for the runs (which are functions of (run1, run2)
+# time. The run_factor terms are step functions that are zero when not in the
+# run, 1 when in the run.
+f = Formula([run_1_coder]) * f1 + Formula([run_2_coder]) * f2 + run_factor
+
+# Now, we evaluate the formula.  So far we've been entirely symbolic.  Now we
+# start to think about the values at which we want to evaluate our symbolic
+# formula.
+
+# We'll use these values for time within run 1.  The times are in seconds from
+# the beginning of run 1.  In our case run 1 was 20 seconds long. 101 below
+# gives 101 values from 0 to 20 including the endpoints, giving a dt of 0.2.
+tval1 = np.linspace(0, 20, 101)
+# run 2 lasts 10 seconds.  These are the times in terms of the start of run 2.
+tval2 = np.linspace(0, 10, 51)
+
+# We pad out the tval1 / tval2 time vectors with zeros corresponding to the
+# TRs in run 2 / run 1.
 ttval1 = np.hstack([tval1, np.zeros(tval2.shape)])
 ttval2 = np.hstack([np.zeros(tval1.shape), tval2])
-session = np.array([1]*tval1.shape[0] + [2]*tval2.shape[0])
-
-f = f.subs(h1, hrf.glover)
-f = f.subs(h2, hrf.glover)
+# The arrays above now have 152=101+51 rows...
 
-# Create the recarray that will be used to create the design matrix
+# Vector of run numbers for each time point (with values 1 or 2)
+run_no = np.array([1]*tval1.shape[0] + [2]*tval2.shape[0])
 
-rec = np.array([(t1,t2, s) for t1, t2, s in zip(ttval1, ttval2, session)],
+# Create the recarray that will be used to create the design matrix. The
+# recarray gives the actual values for the symbolic terms in the formulae.  In
+# our case the terms are t1, t2, and the (indicator coding) terms from the run
+# factor.
+rec = np.array([(tv1, tv2, s) for tv1, tv2, s in zip(ttval1, ttval2, run_no)],
                np.dtype([('t1', np.float),
                          ('t2', np.float),
-                         ('sess', np.int)]))
+                         ('run', np.int)]))
 
 # The contrast we care about
-
-# It would be a good idea to be able to create
-# the contrast from the Formula, "f" above,
-# applying all of f's aliases to it....
-
-sess_1 = sess_factor.get_term(1)
-sess_2 = sess_factor.get_term(2)
-contrast = Formula([sess_1*c11-sess_2*c12])
-contrast = contrast.subs(h1, hrf.glover)
-contrast = contrast.subs(h2, hrf.glover)
+contrast = Formula([run_1_coder * c11 - run_2_coder * c12])
 
 # # Create the design matrix
-
 X = f.design(rec, return_float=True)
 
+# Show ourselves the design space covered by the contrast, and the corresponding
+# contrast matrix
 preC = contrast.design(rec, return_float=True)
-
+# C is the matrix such that preC = X.dot(C.T)
 C = np.dot(np.linalg.pinv(X), preC)
 print C
 
-nonzero = np.nonzero(np.fabs(C) >= 0.1)[0]
+# We can also get this by passing the contrast into the design creation.
+X, c = f.design(rec, return_float=True, contrasts=dict(C=contrast))
+assert np.allclose(C, c['C'])
+
+# Show the names of the non-trivial elements of the contrast
+nonzero = np.nonzero(np.fabs(C) >= 1e-5)[0]
 print (f.dtype.names[nonzero[0]], f.dtype.names[nonzero[1]])
-print ((sess_2*c12).subs(h2, hrf.glover), (sess_1*c11).subs(h1, hrf.glover))
+print ((run_1_coder * c11), (run_2_coder * c12))
