def run(args):
    verbose(1, "Loading %d result files" % len(args.data))
    # TODO: support hdf5 datasets
    nis = [nib.load(f) for f in args.data]
    data = np.asarray([ni.get_data() for ni in nis])

    if args.mask:
        mask = nib.load(args.mask).get_data()
        out_of_mask = mask == 0
    else:
        # just take where no voxel had a value
        out_of_mask = np.sum(data != 0, axis=0)==0

    t, p = ttest_1samp(data, popmean=args.chance_level, axis=0,
                       alternative=args.alternative)

    if args.stat == 'z':
        if args.alternative == 'two-sided':
            s = stats.norm.isf(p/2)
        else:
            s = stats.norm.isf(p)
        # take the sign of the original t
        s = np.abs(s) * np.sign(t)
    elif args.stat == 'p':
        s = p
    elif args.stat == 't':
        s = t
    else:
        raise ValueError('WTF you gave me? have no clue about %r' % (args.stat,))

    s[out_of_mask] = 0

    verbose(1, "Saving to %s" % args.output)
    nib.Nifti1Image(s, None, header=nis[0].header).to_filename(args.output)
    return s

 def test_verbose_below(self):
     """Test if outputs at lower levels and indents
     by default with spaces
     """
     verbose(2, self.msg)
     self.failUnlessEqual(self.sout.getvalue(),
                          "  %s\n" % self.msg)

 def test_verbose_indent(self):
     """Test indent symbol
     """
     verbose.indent = "."
     verbose(2, self.msg)
     self.failUnlessEqual(self.sout.getvalue(), "..%s\n" % self.msg)
     verbose.indent = " "            # restore

def run(args):
    print args.data
    dss = [arg2ds(d)[:,:100] for d in args.data]
    verbose(1, "Loaded %i input datasets" % len(dss))
    if __debug__:
        for i, ds in enumerate(dss):
            debug('CMDLINE', "dataset %i: %s" % (i, str(ds)))
    # TODO at this point more check could be done, e.g. ref_ds > len(dss)
    # assemble parameters
    params = dict([(param, getattr(args, param)) for param in _supported_parameters])
    if __debug__:
        debug('CMDLINE', "configured parameters: '%s'" % params)
    # assemble CAs
    enabled_ca = [ca for ca in _supported_cas if getattr(args, ca)]
    if __debug__:
        debug('CMDLINE', "enabled conditional attributes: '%s'" % enabled_ca)
    hyper = Hyperalignment(enable_ca=enabled_ca,
                           alignment=ProcrusteanMapper(svd='dgesvd',
                                                       space='commonspace'),
                           **params)
    verbose(1, "Running hyperalignment")
    promappers = hyper(dss)
    verbose(2, "Alignment reference is dataset %i" % hyper.ca.chosen_ref_ds)
    verbose(1, "Writing output")
    # save on memory and remove the training data
    del dss
    if args.commonspace:
        if __debug__:
            debug('CMDLINE', "write commonspace as hdf5")
        h5save('%s%s.hdf5' % (args.output_prefix,
                              _output_specs['commonspace']['output_suffix']),
               hyper.commonspace,
               compression=args.hdf5_compression)
    for ca in _supported_cas:
        if __debug__:
            debug('CMDLINE', "check conditional attribute: '%s'" % ca)
        if getattr(args, ca):
            if __debug__:
                debug('CMDLINE', "store conditional attribute: '%s'" % ca)
            np.savetxt('%s%s' % (args.output_prefix,
                                 _supported_cas[ca]['output_suffix']),
                       hyper.ca[ca].value.samples)
    if args.store_transformation:
        for i, pm in enumerate(promappers):
            if __debug__:
                debug('CMDLINE', "store mapper %i: %s" % (i, str(pm)))
            h5save('%s%s.hdf5' % (args.output_prefix, '_map%.3i' % i),
                   pm, compression=args.hdf5_compression)
    if args.transform:
        tdss, dss = _transform_dss(args.transform, promappers, args)
        del dss
        verbose(1, "Store transformed datasets")
        for i, td in enumerate(tdss):
            if __debug__:
                debug('CMDLINE', "store transformed data %i: %s" % (i, str(td)))
            h5save('%s%s.hdf5' % (args.output_prefix, '_transformed%.3i' % i),
                   td, compression=args.hdf5_compression)

def _transform_dss(srcs, mappers, args):
    if __debug__:
        debug('CMDLINE', "loading to-be-transformed data from %s" % srcs)
    dss = [arg2ds(d) for d in srcs]
    verbose(1, "Loaded %i to-be-transformed datasets" % len(dss))
    if __debug__:
        debug('CMDLINE', "transform datasets")
    tdss = [ mappers[i].forward(td) for i, td in enumerate(dss)]
    return tdss, dss

def run(args):
    """Run it"""
    verbose(1, "Loading %d result files" % len(args.data))

    filetype_in = guess_backend(args.data[0])

    if filetype_in == 'nifti':
        dss = [fmri_dataset(f) for f in args.data]
    elif filetype_in == 'hdf5':
        dss = [h5load(f) for f in args.data]
    data = np.asarray([d.samples[args.isample] for d in dss])

    if args.mask:
        filetype_mask = guess_backend(args.mask)
        if filetype_mask == 'nifti':
            mask = nib.load(args.mask).get_data()
        elif filetype_mask == 'hdf5':
            mask = h5load(args.mask).samples
        out_of_mask = mask == 0
    else:
        # just take where no voxel had a value
        out_of_mask = np.sum(data != 0, axis=0)==0

    t, p = ttest_1samp(data, popmean=args.chance_level, axis=0,
                       alternative=args.alternative)

    if args.stat == 'z':
        if args.alternative == 'two-sided':
            s = stats.norm.isf(p/2)
        else:
            s = stats.norm.isf(p)
        # take the sign of the original t
        s = np.abs(s) * np.sign(t)
    elif args.stat == 'p':
        s = p
    elif args.stat == 't':
        s = t
    else:
        raise ValueError('WTF you gave me? have no clue about %r' % (args.stat,))

    if s.shape != out_of_mask.shape:
        try:
            out_of_mask = out_of_mask.reshape(s.shape)
        except ValueError:
            raise ValueError('Cannot use mask of shape {0} with '
                             'data of shape {1}'.format(out_of_mask.shape, s.shape))
    s[out_of_mask] = 0

    verbose(1, "Saving to %s" % args.output)
    filetype_out = guess_backend(args.output)
    if filetype_out == 'nifti':
        map2nifti(dss[0], data=s).to_filename(args.output)
    else:  # filetype_out is hdf5
        s = Dataset(np.atleast_2d(s), fa=dss[0].fa, a=dss[0].a)
        h5save(args.output, s)
    return s

 def test_cr(self):
     """Test if works fine with carriage return (cr) symbol"""
     verbose(2, self.msg, cr=True)
     verbose(2, "rewrite", cr=True)
     verbose(1, "rewrite 2", cr=True)
     verbose(1, " add", cr=False, lf=False)
     verbose(1, " finish")
     target = '\r  %s\r              \rrewrite' % self.msg + \
              '\r       \rrewrite 2 add finish\n'
     self.failUnlessEqual(self.sout.getvalue(), target)

def run(args):
    ds = arg2ds(args.data)
    verbose(3, "Concatenation yielded %i samples with %i features" % ds.shape)
    if args.numpy_xfm is not None:
        from mvpa2.mappers.fx import FxMapper

        fx, axis = args.numpy_xfm
        mapper = FxMapper(axis, fx)
        ds = ds.get_mapped(mapper)
    info_fx[args.report](ds, args)

    def _set_active(self, active):
        """Set active logging set
        """
        # just unique entries... we could have simply stored Set I guess,
        # but then smth like debug.active += ["BLAH"] would not work
        from mvpa2.base import verbose
        self.__active = []
        registered_keys = self.__registered.keys()
        for item in list(set(active)):
            if item == '':
                continue
            if isinstance(item, basestring):
                if item in ['?', 'list', 'help']:
                    self.print_registered(detailed=(item != '?'))
                    raise SystemExit(0)
                if item.upper() == "ALL":
                    verbose(2, "Enabling all registered debug handlers")
                    self.__active = registered_keys
                    break
                # try to match item as it is regexp
                regexp_str = "^%s$" % item
                try:
                    regexp = re.compile(regexp_str)
                except:
                    raise ValueError, \
                          "Unable to create regular expression out of  %s" % item
                matching_keys = filter(regexp.match, registered_keys)
                toactivate = matching_keys
                if len(toactivate) == 0:
                    ids = self.registered.keys()
                    ids.sort()
                    raise ValueError, \
                          "Unknown debug ID '%s' was asked to become active," \
                          " or regular expression '%s' did not get any match" \
                          " among known ids: %s" \
                          % (item, regexp_str, ids)
            else:
                toactivate = [item]

            # Lets check if asked items are known
            for item_ in toactivate:
                if not (item_ in registered_keys):
                    raise ValueError, \
                          "Unknown debug ID %s was asked to become active" \
                          % item_
            self.__active += toactivate

        self.__active = list(set(self.__active)) # select just unique ones
        self.__maxstrlength = max([len(str(x)) for x in self.__active] + [0])
        if len(self.__active):
            verbose(2, "Enabling debug handlers: %s" % `self.__active`)

def run(args):
    import pylab as pl
    from mvpa2.base import verbose

    # segments x [subjects x timepoints x properties]
    data = [np.array(s) for s in args.segment]

    # put in standard property order: first translation, then rotation
    if args.estimate_order == 'rottrans':
        data = [d[:, :, (3, 4, 5, 0, 1, 2)] for d in data]

    # convert rotations, now known to be last
    if args.rad2deg:
        for d in data:
            v = d[:, :, 3:]
            np.rad2deg(v, v)

    # and plot
    # figure setup
    fig = pl.figure(figsize=(12, 5))
    # translation
    ax = pl.subplot(211)
    outlier = motionqc_plot(
        [d[..., :3] for d in data],
        args.outlier_minthresh,
        args.outlier_stdthresh,
        "translation\nestimate L2-norm")
    if outlier:
        verbose(
            0,
            "Detected per-segment translation outlier input samples {0} (zero-based)".format(
                outlier))
    # rotation
    ax = pl.subplot(212)
    outlier = motionqc_plot(
        [d[..., 3:] for d in data],
        args.outlier_minthresh,
        args.outlier_stdthresh,
        "rotation\nestimate L2-norm")
    if outlier:
        verbose(
            0,
            "Detected per-segment rotation outlier input samples {0} (zero-based)".format(
                outlier))

    if args.savefig is None:
        pl.show()
    else:
        pl.savefig(args.savefig[0])

 def test_no_lf(self):
     """Test if it works fine with no newline (LF) symbol"""
     verbose(2, self.msg, lf=False)
     verbose(2, " continue ", lf=False)
     verbose(2, "end")
     verbose(0, "new %s" % self.msg)
     self.failUnlessEqual(self.sout.getvalue(),
                          "  %s continue end\nnew %s\n" % \
                          (self.msg, self.msg))

 def bxplot(stats, label):
     stats = concat_ts_boxplot_stats(stats)
     verbose(0, "List of outlier time series follows (if any)")
     for i, run in enumerate([np.where(np.sum(np.logical_not(o.mask), axis=0))
                           for o in stats[1]]):
         sids = run[0]
         if len(sids):
             verbose(0, "%s r%.3i: %s"
                        % (label, i + 1, [s + 1 for s in sids]))
     timeseries_boxplot(stats[0]['median'],
             mean=stats[0]['mean'], std=stats[0]['std'], n=stats[0]['n'],
             min=stats[0]['min'], max=stats[0]['max'],
             p25=stats[0]['p25'], p75=stats[0]['p75'],
             outlierd=stats[1], segment_sizes=segment_sizes)
     pl.title(label)
     xp, xl = pl.xticks()
     pl.xticks(xp, ['' for i in xl])
     pl.xlim((0, len(stats[0]['n'])))
     pl.ylabel(plt_props[label])

def ds2hdf5(ds, fname, compression=None):
    """Save one or more datasets into an HDF5 file.

    Parameters
    ----------
    ds : Dataset or list(Dataset)
      One or more datasets to store
    fname : str
      Filename of the output file. If it doesn't end with '.hdf5', such an
      extension will be appended.
    compression : {'gzip','lzf','szip'} or 1-9
      compression type for HDF5 storage. Available values depend on the specific
      HDF5 installation.
    """
    # this one doesn't actually check what it stores
    from mvpa2.base.hdf5 import h5save
    if not fname.endswith('.hdf5'):
        fname = '%s.hdf5' % fname
    verbose(1, "Save dataset to '%s'" % fname)
    h5save(fname, ds, mkdir=True, compression=compression)

def run(args):
    dss = hdf2ds(args.data)
    verbose(3, 'Loaded %i dataset(s)' % len(dss))
    ds = vstack(dss)
    verbose(3, 'Concatenation yielded %i samples with %i features' % ds.shape)
    # get CV instance
    cv = get_crossvalidation_instance(
            args.learner, args.partitioner, args.errorfx, args.sampling_repetitions,
            args.learner_space, args.balance_training, args.permutations,
            args.avg_datafold_results, args.prob_tail)
    res = cv(ds)
    # some meaningful output
    # XXX make condition on classification analysis only?
    print cv.ca.stats
    print 'Results\n-------'
    if args.permutations > 0:
        nprob =  cv.ca.null_prob.samples
    if res.shape[1] == 1:
        # simple result structure
        if args.permutations > 0:
            p=', p-value (%s tail)' % args.prob_tail
        else:
            p=''
        print 'Fold, Result%s' % p
        for i in xrange(len(res)):
            if args.permutations > 0:
                p = ', %f' % nprob[i, 0]
            else:
                p = ''
            print '%s, %f%s' % (res.sa.cvfolds[i], res.samples[i, 0], p)
    # and store
    ds2hdf5(res, args.output, compression=args.hdf5_compression)
    if args.permutations > 0:
        if args.output.endswith('.hdf5'):
            args.output = args.output[:-5]
        ds2hdf5(cv.ca.null_prob, '%s_nullprob' % args.output,
                compression=args.hdf5_compression)
    return res

def run(args):
    dss = hdf2ds(args.data)
    verbose(3, 'Loaded %i dataset(s)' % len(dss))
    ds = vstack(dss)
    verbose(3, 'Concatenation yielded %i samples with %i features' % ds.shape)
    # slicing
    sliceme = {'samples': slice(None), 'features': slice(None)}
    # indices
    for opt, col, which in ((args.samples_by_index, ds.sa, 'samples'),
                     (args.features_by_index, ds.fa, 'features')):
        if opt is None:
            continue
        if len(opt) == 1 and opt[0].count(':'):
            # slice spec
            arg = opt[0].split(':')
            spec = []
            for a in arg:
                if not len(a):
                    spec.append(None)
                else:
                    spec.append(int(a))
            sliceme[which] = slice(*spec)
        else:
            # actual indices
            sliceme[which] = [int(o) for o in opt]
    # attribute evaluation
    for opt, col, which in ((args.samples_by_attr, ds.sa, 'samples'),
                     (args.features_by_attr, ds.fa, 'features')):
        if opt is None:
            continue
        sliceme[which] = _eval_attr_expr(opt, col)

    # apply selection
    ds = ds.__getitem__((sliceme['samples'], sliceme['features']))
    verbose(1, 'Selected %i samples with %i features' % ds.shape)

    # strip attributes
    for attrarg, col, descr in ((args.strip_sa, ds.sa, 'sample '),
                                (args.strip_fa, ds.fa, 'feature '),
                                (args.strip_da, ds.a, '')):
        if not attrarg is None:
            for attr in attrarg:
                try:
                    del col[attr]
                except KeyError:
                    warning("dataset has no %sattribute '%s' to remove"
                            % (descr, attr))
    # and store
    ds2hdf5(ds, args.output, compression=args.hdf5_compression)
    return ds

def process_common_dsattr_opts(ds, args):
    """Goes through an argument namespace and processes attribute options"""
    # legacy support
    if not args.add_sa_attr is None:
        from mvpa2.misc.io.base import SampleAttributes
        smpl_attrs = SampleAttributes(args.add_sa_attr)
        for a in ('targets', 'chunks'):
            verbose(2, "Add sample attribute '%s' from sample attributes file"
                       % a)
            ds.sa[a] = getattr(smpl_attrs, a)
    # loop over all attribute configurations that we know
    attr_cfgs = (# var, dst_collection, loader
            ('--add-sa', args.add_sa, ds.sa, _load_from_cmdline),
            ('--add-fa', args.add_fa, ds.fa, _load_from_cmdline),
            ('--add-sa-txt', args.add_sa_txt, ds.sa, _load_from_txt),
            ('--add-fa-txt', args.add_fa_txt, ds.fa, _load_from_txt),
            ('--add-sa-npy', args.add_sa_npy, ds.sa, _load_from_npy),
            ('--add-fa-npy', args.add_fa_npy, ds.fa, _load_from_npy),
        )
    for varid, srcvar, dst_collection, loader in attr_cfgs:
        if not srcvar is None:
            for spec in srcvar:
                attr_name = spec[0]
                if not len(spec) > 1:
                    raise argparse.ArgumentTypeError(
                        "%s option need at least two values " % varid +
                        "(attribute name and source filename (got: %s)" % spec)
                if dst_collection is ds.sa:
                    verbose(2, "Add sample attribute '%s' from '%s'"
                               % (attr_name, spec[1]))
                else:
                    verbose(2, "Add feature attribute '%s' from '%s'"
                               % (attr_name, spec[1]))
                attr = loader(spec[1:])
                try:
                    dst_collection[attr_name] = attr
                except ValueError, e:
                    # try making the exception more readable
                    e_str = str(e)
                    if e_str.startswith('Collectable'):
                        raise ValueError('attribute %s' % e_str[12:])
                    else:
                        raise e

def extract_lsall(data, TR, time_res,
                  hrf_gen, F,
                  good_ons,
                  good_evs,
                  desmat,
                  extract_evs=None):
    ntp, nvox = data.shape

    hrf = hrf_gen(time_res)
    # Set up the high time-resolution design matrix
    time_up = N.arange(0, TR*ntp+time_res, time_res)

    all_onsets = []
    all_durations = []
    all_conds = []  # condition marker

    if extract_evs is None:
        extract_evs = range(len(good_evs))

    nuisance_evs = sorted(list(set(range(desmat.mat.shape[1])).difference(
        [good_evs[e] for e in extract_evs])))

    for e in extract_evs:
        ev = good_evs[e]
        all_onsets    = N.hstack((all_onsets,    good_ons[e].onsets))
        all_durations = N.hstack((all_durations, good_ons[e].durations))
        # yoh: ad-hoc warning -- it is marking with (ev/2)+1 (I guess)
        # assuming presence of derivatives EVs
        all_conds     = N.hstack((all_conds,     N.ones(len(good_ons[e].onsets))*((ev/2)+1)))

    #all_onsets=N.round(all_onsets/TR)  # round to nearest TR number
    ntrials = len(all_onsets)
    glm_res_full = N.zeros((nvox, ntrials))
    dm_trials = N.zeros((ntp, ntrials))
    dm_full = []
    for t in range(ntrials):
        verbose(2, "Estimating for trial %d" % t)

        ## yoh: TODO -- filter outside
        ## if all_onsets[t] > max_evtime:
        ##     continue
        # build model for each trial
        dm_trial = N.zeros(len(time_up))
        window_ons = [N.where(time_up==x)[0][0]
                      for x in time_up
                      if all_onsets[t] <= x < all_onsets[t] + all_durations[t]]
        dm_trial[window_ons] = 1
        dm_trial_up = N.convolve(dm_trial, hrf)
        dm_trial_down = dm_trial_up[0:ntp/time_res*TR:(TR/time_res)]
        dm_trials[:, t] = dm_trial_down

    # filter the desmtx, except for the nuisance part (which is already filtered)
    # since it is taken from a loaded FSL
    dm_full = N.dot(F, dm_trials)

    # mean center trials models
    dm_trials -= dm_trials.mean(0)

    if len(nuisance_evs) > 0:
        # and stick nuisance evs if any to the back
        dm_full = N.hstack((dm_full, desmat.mat[:, nuisance_evs]))

    dm_full = N.hstack((dm_full, N.ones((ntp, 1))))
    glm_res_full = N.dot(N.linalg.pinv(dm_full), data.samples)
    glm_res_full = glm_res_full[:ntrials]

    return all_conds, glm_res_full

def plot_scatter(dataXd, mask=None, masked_opacity=0.,
                 labels=None, colors=True,
                 dimcolor=1, title=None, limits='auto',
                 thresholds=None, hint_opacity=0.9,
                 x_jitter=None, y_jitter=None,
                 fig=None,
                 ax_scatter=None, ax_hist_x=None, ax_hist_y=None,
                 bp_location='scatter',
                 xlim=None, ylim=None,
                 rasterized=None,
                 uniq=False,
                 include_stats=False,
                 ):
    """
    Parameters
    ----------
    dataXd: array
      The volumetric (or not) data to plot where first dimension
      should only have 2 items
    mask: array, optional
      Additional mask to specify which values do not consider to plot.
      By default values with 0s in both dimensions are not plotted.
    masked_opacity: float, optional
      By default masked out values are not plotted at all.  Value in
      (0,1] will make them visible with this specified opacity
    labels: list of str, optional
      Labels to place for x and y axes
    colors: bool or string or colormap, optional
      Either to use colors to associate with physical location and
      what colormap to use (jet by default if colors=True)
    dimcolor: int
      If `colors`, then which dimension (within given 3D volume) to
      "track"
    limits: 'auto', 'same', 'per-axis' or (min, max)
      Limits for axes: when 'auto' if data ranges overlap is more than
      50% of the union range, 'same' is considered.  When 'same' --
      the same limits on both axes as determined by data.  If
      two-element tuple or list is provided, then that range is
      applied to both axes.
    hint_opacity: float, optional
      If `colors` is True, to then a "slice" of the volumetric data
      is plotted in the specified opacity to hint about the location
      of points in the original Xd data in `dimcolor` dimension
    x_jitter: float, optional
      Half-width of uniform noise added to x values.  Might be useful if data
      is quantized so it is valuable to jitter points a bit.
    y_jitter: float, optional
      Half-width of uniform noise added to y values.  Might be useful if data
      is quantized so it is valuable to jitter points a bit
    fig : Figure, optional
      Figure to plot on, otherwise new one created
    ax_*: axes, optional
      Axes for the scatter plot and histograms. If none of them is specified
      (which is the default) then 'classical' plot is rendered with histograms
      above and to the right
    bp_location: ('scatter', 'hist', None), optional
      Where to place boxplots depicting data range
    xlim: tuple, optional
    ylim: tuple, optional
      To fix plotted range
    rasterized: bool, optional
      Passed to scatter call, to allow rasterization of heavy scatter plots
    uniq: bool, optional
      Plot uniq values (those present in one but not in the other) along
      each axis with crosses
    include_stats: bool, optional
      Whether to report additional statistics on the data. Stats are also
      reported via verbose at level 2
    """
    if len(dataXd) != 2:
        raise ValueError("First axis of dataXd can only have two dimensions, "
                         "got {0}".format(len(dataXd)))
    dataXd = np.asanyarray(dataXd)      # TODO: allow to operate on list of arrays to not waste RAM/cycles
    data = dataXd.reshape((2, -1))
    if dataXd.ndim < 5:
        ntimepoints = 1
    elif dataXd.ndim == 5:
        ntimepoints = dataXd.shape[-1]
    else:
        raise ValueError("Do not know how to handle data with %d dimensions" % (dataXd.ndim - 1))
    if x_jitter or y_jitter:
        data = data.copy()              # lazy and wasteful
        def jitter_me(x, w):
            x += np.random.uniform(-w, w, size=data.shape[-1])
        if x_jitter:
            jitter_me(data[0, :], x_jitter)
        if y_jitter:
            jitter_me(data[1, :], y_jitter)

    finites = np.isfinite(data)
    nz = np.logical_and(data != 0, finites)
    # TODO : avoid doing data !=0 and just use provided utter mask
    #nz[:, 80000:] = False # for quick testing

    nzsum = np.sum(nz, axis=0)

    intersection = nzsum == 2
    # for coloring we would need to know all the indices
    union = nzsum > 0
    x, y = datainter = data[:, intersection]
#.........这里部分代码省略.........

    if args.atlasFile is None:
        if args.atlasPath is None:
            args.atlasPath = KNOWN_ATLASES[args.atlasName]
        args.atlasFile = args.atlasPath % ( {'name': args.atlasName} )

    akwargs_common = {}
    if args.atlasImageFile:
        akwargs_common['image_file'] = args.atlasImageFile

    if not args.forbidDirectMapping \
           and niftiInput is not None and not args.transformationFile:
        akwargs = {'resolution': niftiInput.get_header().get_zooms()[0]}
        query_voxel = True   # if we can query directly by voxel, do so

        akwargs.update(akwargs_common)
        verbose(1, "Will attempt direct mapping from input voxels into atlas "
                   "voxels at resolution %.2f" % akwargs['resolution'])

        atlas = Atlas(args.atlasFile, **akwargs)

        # verify that we got the same qforms in atlas and in the data file
        if atlas.space != args.inputSpace:
            verbose(0,
                "Cannot do direct mapping between input image in %s space and"
                " atlas in %s space. Use -I switch to override input space if"
                " it misspecified, or use -T to provide transformation. Trying"
                " to proceed" %(args.inputSpace, atlas.space))
            query_voxel = False
        elif not (niftiInput.get_header().get_qform() == atlas._image.get_header().get_qform()).all():
            if args.atlasImageFile is None:
                warning(
                    "Cannot do direct mapping between files with different qforms."

def run(args):
    if os.path.isfile(args.payload) and args.payload.endswith('.py'):
        measure = script2obj(args.payload)
    elif args.payload == 'cv':
        if args.cv_learner is None or args.cv_partitioner is None:
            raise ValueError('cross-validation payload requires --learner and --partitioner')
        # get CV instance
        measure = get_crossvalidation_instance(
                    args.cv_learner, args.cv_partitioner, args.cv_errorfx,
                    args.cv_sampling_repetitions, args.cv_learner_space,
                    args.cv_balance_training, args.cv_permutations,
                    args.cv_avg_datafold_results, args.cv_prob_tail)
    else:
        raise RuntimeError("this should not happen")
    ds = arg2ds(args.data)
    if args.ds_preproc_fx is not None:
        ds = args.ds_preproc_fx(ds)
    # setup neighborhood
    # XXX add big switch to allow for setting up surface-based neighborhoods
    from mvpa2.misc.neighborhood import IndexQueryEngine
    qe = IndexQueryEngine(**dict(args.neighbors))
    # determine ROIs
    rids = None     # all by default
    aggregate_fx = args.aggregate_fx
    if args.roi_attr is not None:
        # first figure out which roi features should be processed
        if len(args.roi_attr) == 1 and args.roi_attr[0] in ds.fa.keys():
            # name of an attribute -> pull non-zeroes
            rids = ds.fa[args.roi_attr[0]].value.nonzero()[0]
        else:
            # an expression?
            from .cmd_select import _eval_attr_expr
            rids = _eval_attr_expr(args.roi_attr, ds.fa).nonzero()[0]

    seed_ids = None
    if args.scatter_rois is not None:
        # scatter_neighborhoods among available ids if was requested
        from mvpa2.misc.neighborhood import scatter_neighborhoods
        attr, nb = args.scatter_rois
        coords = ds.fa[attr].value
        if rids is not None:
            # select only those which were chosen by ROI
            coords = coords[rids]
        _, seed_ids = scatter_neighborhoods(nb, coords)
        if aggregate_fx is None:
            # no custom one given -> use default "fill in" function
            aggregate_fx = _fill_in_scattered_results
            if args.enable_ca is None:
                args.enable_ca = ['roi_feature_ids']
            elif 'roi_feature_ids' not in args.enable_ca:
                args.enable_ca += ['roi_feature_ids']

    if seed_ids is None:
        roi_ids = rids
    else:
        if rids is not None:
            # we had to sub-select by scatterring among available rids
            # so we would need to get original ids
            roi_ids = rids[seed_ids]
        else:
            # scattering happened on entire feature-set
            roi_ids = seed_ids

    verbose(3, 'Attempting %i ROI analyses'
               % ((roi_ids is None) and ds.nfeatures or len(roi_ids)))

    from mvpa2.measures.searchlight import Searchlight

    sl = Searchlight(measure,
                     queryengine=qe,
                     roi_ids=roi_ids,
                     nproc=args.nproc,
                     results_backend=args.multiproc_backend,
                     results_fx=aggregate_fx,
                     enable_ca=args.enable_ca,
                     disable_ca=args.disable_ca)
    # XXX support me too!
    #                 add_center_fa
    #                 tmp_prefix
    #                 nblocks
    #                 null_dist
    # run
    res = sl(ds)
    if (seed_ids is not None) and ('mapper' in res.a):
        # strip the last mapper link in the chain, which would be the seed ID selection
        res.a['mapper'] = res.a.mapper[:-1]
    # XXX create more output
    # and store
    ds2hdf5(res, args.output, compression=args.hdf5_compression)
    return res