Merge branch 'eg-refactor-regression' into main-master

matthew-brett · matthew-brett · commit 1d4a94f8de04 · 2012-07-19T16:55:55.000-07:00
* eg-refactor-regression:
  BF: fix parallel example to run permutation test
  update parallel fiac example for IPython.parallel
  BF+DOC: fixing permutation test example
diff --git a/examples/fiac/fiac_example.py b/examples/fiac/fiac_example.py
@@ -2,11 +2,19 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Example analyzing the FIAC dataset with NIPY.
 
+* Single run models with per-voxel AR(1).
+* Cross-run, within-subject models with optimal effect estimates.
+* Cross-subject models using fixed / random effects variance ratios.
+* Permutation testing for inference on cross-subject result.
+
+See ``parallel_run.py`` for a rig to run these analysis in parallel using the
+IPython parallel machinery.
+
 This script needs the pre-processed FIAC data.  See ``README.txt`` and
 ``fiac_util.py`` for details.
 
 See ``examples/labs/need_data/first_level_fiac.py`` for an alternative approach
-to this analysis.
+to some of these analyses.
 """
 #-----------------------------------------------------------------------------
 # Imports
@@ -281,17 +289,17 @@ def fixed_effects(subj, design):
     fixdir = pjoin(rootdir, "fixed")
     # Fetch results images from run estimations
     results = futil.results_table(path_dict)
-    # Get our hands on the relevant coordmap to
-    # save our results
+    # Get our hands on the relevant coordmap to save our results
     coordmap = futil.load_image_fiac("fiac_%02d" % subj,
                                      "wanatomical.nii").coordmap
     # Compute the "fixed" effects for each type of contrast
     for con in results:
         fixed_effect = 0
         fixed_var = 0
         for effect, sd in results[con]:
-            effect = load_image(effect); sd = load_image(sd)
-            var = np.array(sd)**2
+            effect = load_image(effect).get_data()
+            sd = load_image(sd).get_data()
+            var = sd ** 2
 
             # The optimal, in terms of minimum variance, combination of the
             # effects has weights 1 / var
@@ -332,9 +340,19 @@ def group_analysis(design, contrast):
 
     # Which subjects have this (contrast, design) pair?
     subj_con_dirs = futil.subj_des_con_dirs(design, contrast)
-
-    sd = array([array(load_image(pjoin(s, "sd.nii"))) for s in subj_con_dirs])
-    Y = array([array(load_image(pjoin(s, "effect.nii"))) for s in subj_con_dirs])
+    if len(subj_con_dirs) == 0:
+        raise ValueError('No subjects for %s, %s' % (design, contrast))
+
+    # Assemble effects and sds into 4D arrays
+    sds = []
+    Ys = []
+    for s in subj_con_dirs:
+        sd_img = load_image(pjoin(s, "sd.nii"))
+        effect_img = load_image(pjoin(s, "effect.nii"))
+        sds.append(sd_img.get_data())
+        Ys.append(effect_img.get_data())
+    sd = array(sds)
+    Y = array(Ys)
 
     # This function estimates the ratio of the fixed effects variance
     # (sum(1/sd**2, 0)) to the estimated random effects variance
@@ -378,7 +396,9 @@ def group_analysis_signs(design, contrast, mask, signs=None):
     ----------
     design: one of 'block', 'event'
     contrast: str
-    mask: array-like
+        name of contrast to estimate
+    mask : ``Image`` instance or array-like
+        image containing mask, or array-like
     signs: ndarray, optional
          Defaults to np.ones. Should have shape (*,nsubj)
          where nsubj is the number of effects combined in the group analysis.
@@ -390,15 +410,25 @@ def group_analysis_signs(design, contrast, mask, signs=None):
     maxT: np.ndarray, maxima of T statistic within mask, one for each
          vector of signs
     """
-    maska = np.asarray(mask).astype(np.bool)
+    if api.is_image(mask):
+        maska = mask.get_data()
+    else:
+        maska = np.asarray(mask)
+    maska = maska.astype(np.bool)
 
     # Which subjects have this (contrast, design) pair?
     subj_con_dirs = futil.subj_des_con_dirs(design, contrast)
 
-    sd = np.array([np.array(load_image(pjoin(s, "sd.nii")))[:,maska]
-                   for s in subj_con_dirs])
-    Y = np.array([np.array(load_image(pjoin(s, "effect.nii")))[:,maska]
-                  for s in subj_con_dirs])
+    # Assemble effects and sds into 4D arrays
+    sds = []
+    Ys = []
+    for s in subj_con_dirs:
+        sd_img = load_image(pjoin(s, "sd.nii"))
+        effect_img = load_image(pjoin(s, "effect.nii"))
+        sds.append(sd_img.get_data()[maska])
+        Ys.append(effect_img.get_data()[maska])
+    sd = np.array(sds)
+    Y = np.array(Ys)
 
     if signs is None:
         signs = np.ones((1, Y.shape[0]))
@@ -428,22 +458,26 @@ def permutation_test(design, contrast, mask=GROUP_MASK, nsample=1000):
 
     Parameters
     ----------
-    design: one of ['block', 'event']
-    contrast: str
-    nsample: int
+    design: str
+        one of ['block', 'event']
+    contrast : str
+        name of contrast to estimate
+    mask : ``Image`` instance or array-like, optional
+        image containing mask, or array-like
+    nsample: int, optional
+        number of permutations
 
     Returns
     -------
     min_vals: np.ndarray
     max_vals: np.ndarray
     """
-    maska = np.asarray(mask).astype(np.bool)
     subj_con_dirs = futil.subj_des_con_dirs(design, contrast)
-    Y = np.array([np.array(load_image(pjoin(s, "effect.nii")))[:,maska]
-                  for s in subj_con_dirs])
-    nsubj = Y.shape[0]
+    nsubj = len(subj_con_dirs)
+    if nsubj == 0:
+        raise ValueError('No subjects have %s, %s' % (design, contrast))
     signs = 2*np.greater(np.random.sample(size=(nsample, nsubj)), 0.5) - 1
-    min_vals, max_vals = group_analysis_signs(design, contrast, maska, signs)
+    min_vals, max_vals = group_analysis_signs(design, contrast, mask, signs)
     return min_vals, max_vals
 
 
diff --git a/examples/fiac/parallel_run.py b/examples/fiac/parallel_run.py
@@ -1,33 +1,38 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Script to run the main analyses in parallel, using the IPython machinery.
+
+See ``fiac_example.py``.
 """
 #-----------------------------------------------------------------------------
 # Imports
 #-----------------------------------------------------------------------------
 
-from IPython.kernel import client
+import os
+
+import numpy as np
 
-# Local imports
-from fiac_example import GROUP_MASK
+from IPython import parallel
 
 #-----------------------------------------------------------------------------
 # Utility functions
 #-----------------------------------------------------------------------------
 
-def setup_mec():
-    """Get a multiengineclient and initialize it.
+_client = None
+def setup_client():
+    """Get a Client and initialize it.
 
     This assumes that all nodes see a shared filesystem.
     """
-    mec = client.MultiEngineClient()
-    # Ensure each engine is in the same directory as we are
-    mydir = os.path.split(os.path.abspath(__file__))[0]
-    mec.execute('''
-import os
-os.chdir(%s)
-    ''' % mydir)
-    return mec
+    global _client
+    if _client is None:
+        _client = parallel.Client()
+        mydir = os.path.split(os.path.abspath(__file__))[0]
+        def cd(path):
+            import os
+            os.chdir(path)
+        _client[:].apply_sync(cd, mydir)
+    return _client
 
 
 def getruns():
@@ -53,66 +58,72 @@ def getvals():
 
 def fitruns():
     """Run the basic model fit."""
-    mec = setup_mec()
-    runs = list(getruns())
-    mec.scatter('runs', runs)
-    mec.execute('''
-import fiac_example
-for subj, run in runs:
-    try:
-        fiac_example.run_model(subj, run)
-    except IOError:
-        pass
-    ''')
+    rc = setup_client()
+    view = rc.load_balanced_view()
+    i_s, j_s = zip(*getruns())
+
+    def _fit(subj, run):
+        import fiac_example
+        try:
+            return fiac_example.run_model(subj, run)
+        except IOError:
+            pass
+
+    return view.map(_fit, i_s, j_s)
 
 
 def fitfixed():
     """Run the fixed effects analysis for all subjects."""
-    mec = setup_mec()
-    mec.scatter('subjects', range(16))
-    mec.execute('''
-import fiac_example
-for s in subjects:
-    try:
-        fiac_example.fixed_effects(s, "block")
-    except IOError:
-        pass
-    try:
-        fiac_example.fixed_effects(s, "event")
-    except IOError:
-        pass
-''')
+    rc = setup_client()
+    view = rc.load_balanced_view()
+    subjects = range(16)
+
+    def _fit(subject):
+        import fiac_example
+        try:
+            fiac_example.fixed_effects(subject, "block")
+        except IOError:
+            pass
+        try:
+            fiac_example.fixed_effects(subject, "event")
+        except IOError:
+            pass
+
+    return view.map(_fit, subjects)
 
 
 def fitgroup():
     """Run the group analysis"""
-    mec = setup_mec()
-    group_vals = list( getvals() )
-    
-    mec.scatter('group_vals', group_vals)
-    mec.execute('''
-import fiac_example
-for d, c in group_vals:
-    fiac_example.group_analysis(d, c)
-''')
+    rc = setup_client()
+    view = rc.load_balanced_view()
+    d_s, c_s = zip(*getvals())
+
+    def _fit(d, c):
+        import fiac_example
+        return fiac_example.group_analysis(d, c)
+
+    return view.map(_fit, d_s, c_s)
 
 
 def run_permute_test(design, contrast, nsample=1000):
-    mec = setup_mec()
-    nnod = len(mec.get_ids())
-    ns_nod = nsample/nnod
-    mec.push(ns_nod=ns_nod, design=design,contrast=contrast)
-    mec.execute('''
-import fiac_example
-from fiac_example import GROUP_MASK
-min_vals, max_vals = fiac_example.permutation_test(design, contrast,
-                               GROUP_MASK, ns_nod)
-    ''')
-    min_vals = mec.gather('min_vals')
-    max_vals = mec.gather('max_vals')
-    return min_vals, max_vals
-
-    
+    rc = setup_client()
+    dview = rc[:]
+    nnod = len(dview)
+    # Samples per node.  Round up
+    ns_nod = np.ceil(nsample / float(nnod))
+
+    def _run_test(n, des, con):
+        import fiac_example
+        from fiac_example import GROUP_MASK
+        min_vals, max_vals = fiac_example.permutation_test(des, con,
+                                       GROUP_MASK, n)
+        return min_vals, max_vals
+
+    ar = dview.apply_async(_run_test, ns_nod, design, contrast)
+    min_vals, max_vals = zip(*[r for r in ar])
+    return np.concatenate(min_vals), np.concatenate(max_vals)
+
+
 #-----------------------------------------------------------------------------
 # Script entry point
 #-----------------------------------------------------------------------------
