# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""The fsl module provides classes for interfacing with the `FSL
<http://www.fmrib.ox.ac.uk/fsl/index.html>`_ command line tools.  This
was written to work with FSL version 4.1.4.

    Change directory to provide relative paths for doctests
    >>> import os
    >>> filepath = os.path.dirname( os.path.realpath( __file__ ) )
    >>> datadir = os.path.realpath(os.path.join(filepath, '../../testing/data'))
    >>> os.chdir(datadir)
"""

import os
from glob import glob
import warnings
from shutil import rmtree

import numpy as np

from nipype.interfaces.fsl.base import (FSLCommand, FSLCommandInputSpec)
from nipype.interfaces.base import (load_template, File, traits, isdefined,
                                    TraitedSpec, BaseInterface, Directory,
                                    InputMultiPath, OutputMultiPath,
                                    BaseInterfaceInputSpec)
from nipype.utils.filemanip import (list_to_filename, filename_to_list)
from nibabel import load

warn = warnings.warn
warnings.filterwarnings('always', category=UserWarning)


class Level1DesignInputSpec(BaseInterfaceInputSpec):
    interscan_interval = traits.Float(mandatory=True,
                desc='Interscan  interval (in secs)')
    session_info = traits.Any(mandatory=True,
                desc='Session specific information generated by ``modelgen.SpecifyModel``')
    bases = traits.Either(traits.Dict(traits.Enum('dgamma'), traits.Dict(traits.Enum('derivs'), traits.Bool)),
                          traits.Dict(traits.Enum('gamma'), traits.Dict(traits.Enum('derivs'), traits.Bool)),
                          traits.Dict(traits.Enum('none'), traits.Enum(None)),
                          mandatory=True,
                          desc="name of basis function and options e.g., {'dgamma': {'derivs': True}}")
    model_serial_correlations = traits.Bool(
        desc="Option to model serial correlations using an \
autoregressive estimator (order 1). Setting this option is only \
useful in the context of the fsf file. If you set this to False, you need to repeat \
this option for FILMGLS by setting autocorr_noestimate to True", mandatory=True)
    contrasts = traits.List(
        traits.Either(traits.Tuple(traits.Str,
                                   traits.Enum('T'),
                                   traits.List(traits.Str),
                                   traits.List(traits.Float)),
                      traits.Tuple(traits.Str,
                                   traits.Enum('T'),
                                   traits.List(traits.Str),
                                   traits.List(traits.Float),
                                   traits.List(traits.Float)),
                      traits.Tuple(traits.Str,
                                   traits.Enum('F'),
                                   traits.List(traits.Either(traits.Tuple(traits.Str,
                                                                          traits.Enum('T'),
                                                                          traits.List(traits.Str),
                                                                          traits.List(traits.Float)),
                                                             traits.Tuple(traits.Str,
                                                                          traits.Enum('T'),
                                                                          traits.List(traits.Str),
                                                                          traits.List(traits.Float),
                                                                          traits.List(traits.Float)))))),
        desc="List of contrasts with each contrast being a list of the form - \
[('name', 'stat', [condition list], [weight list], [session list])]. if \
session list is None or not provided, all sessions are used. For F \
contrasts, the condition list should contain previously defined \
T-contrasts.")


class Level1DesignOutputSpec(TraitedSpec):
    fsf_files = OutputMultiPath(File(exists=True),
                     desc='FSL feat specification files')
    ev_files = OutputMultiPath(traits.List(File(exists=True)),
                     desc='condition information files')


class Level1Design(BaseInterface):
    """Generate FEAT specific files

    Examples
    --------

    >>> level1design = Level1Design()
    >>> level1design.inputs.interscan_interval = 2.5
    >>> level1design.inputs.bases = {'dgamma':{'derivs': False}}
    >>> level1design.inputs.session_info = 'session_info.npz'
    >>> level1design.run() # doctest: +SKIP

    """

    input_spec = Level1DesignInputSpec
    output_spec = Level1DesignOutputSpec

    def _create_ev_file(self, evfname, evinfo):
        f = open(evfname, 'wt')
        for i in evinfo:
            if len(i) == 3:
                f.write('%f %f %f\n' % (i[0], i[1], i[2]))
            else:
                f.write('%f\n' % i[0])
        f.close()

    def _create_ev_files(self, cwd, runinfo, runidx, usetd, contrasts, no_bases,
                         do_tempfilter):
        """Creates EV files from condition and regressor information.

           Parameters:
           -----------

           runinfo : dict
               Generated by `SpecifyModel` and contains information
               about events and other regressors.
           runidx  : int
               Index to run number
           usetd   : int
               Whether or not to use temporal derivatives for
               conditions
           contrasts : list of lists
               Information on contrasts to be evaluated
        """
        conds = {}
        evname = []
        ev_hrf = load_template('feat_ev_hrf.tcl')
        ev_none = load_template('feat_ev_none.tcl')
        ev_ortho = load_template('feat_ev_ortho.tcl')
        ev_txt = ''
        # generate sections for conditions and other nuisance
        # regressors
        num_evs = [0, 0]
        for field in ['cond', 'regress']:
            for i, cond in enumerate(runinfo[field]):
                name = cond['name']
                evname.append(name)
                evfname = os.path.join(cwd, 'ev_%s_%d_%d.txt' % (name, runidx,
                                                                 len(evname)))
                evinfo = []
                num_evs[0] += 1
                num_evs[1] += 1
                if field == 'cond':
                    for j, onset in enumerate(cond['onset']):
                        try:
                            amplitudes = cond['amplitudes']
                            if len(amplitudes) > 1:
                                amp = amplitudes[j]
                            else:
                                amp = amplitudes[0]
                        except KeyError:
                            amp = 1
                        if len(cond['duration']) > 1:
                            evinfo.insert(j, [onset, cond['duration'][j], amp])
                        else:
                            evinfo.insert(j, [onset, cond['duration'][0], amp])
                    if no_bases:
                        ev_txt += ev_none.substitute(ev_num=num_evs[0],
                                                     ev_name=name,
                                                     tempfilt_yn=do_tempfilter,
                                                     cond_file=evfname)
                    else:
                        ev_txt += ev_hrf.substitute(ev_num=num_evs[0],
                                                    ev_name=name,
                                                    tempfilt_yn=do_tempfilter,
                                                    temporalderiv=usetd,
                                                    cond_file=evfname)
                    if usetd:
                        evname.append(name + 'TD')
                        num_evs[1] += 1
                elif field == 'regress':
                    evinfo = [[j] for j in cond['val']]
                    ev_txt += ev_none.substitute(ev_num=num_evs[0],
                                                 ev_name=name,
                                                 tempfilt_yn=do_tempfilter,
                                                 cond_file=evfname)
                ev_txt += "\n"
                conds[name] = evfname
                self._create_ev_file(evfname, evinfo)
        # add ev orthogonalization
        for i in range(1, num_evs[0] + 1):
            for j in range(0, num_evs[0] + 1):
                ev_txt += ev_ortho.substitute(c0=i, c1=j)
                ev_txt += "\n"
        # add contrast info to fsf file
        if isdefined(contrasts):
            contrast_header = load_template('feat_contrast_header.tcl')
            contrast_prolog = load_template('feat_contrast_prolog.tcl')
            contrast_element = load_template('feat_contrast_element.tcl')
            contrast_ftest_element = load_template('feat_contrast_ftest_element.tcl')
            contrastmask_header = load_template('feat_contrastmask_header.tcl')
            contrastmask_footer = load_template('feat_contrastmask_footer.tcl')
            contrastmask_element = load_template('feat_contrastmask_element.tcl')
            # add t/f contrast info
            ev_txt += contrast_header.substitute()
            con_names = []
            for j, con in enumerate(contrasts):
                con_names.append(con[0])
            con_map = {}
            ftest_idx = []
            ttest_idx = []
            for j, con in enumerate(contrasts):
                if con[1] == 'F':
                    ftest_idx.append(j)
                    for c in con[2]:
                        if c[0] not in con_map.keys():
                            con_map[c[0]] = []
                        con_map[c[0]].append(j)
                else:
                    ttest_idx.append(j)

            for ctype in ['real', 'orig']:
                for j, con in enumerate(contrasts):
                    if con[1] == 'F':
                        continue
                    tidx = ttest_idx.index(j) + 1
                    ev_txt += contrast_prolog.substitute(cnum=tidx,
                                                         ctype=ctype,
                                                         cname=con[0])
                    count = 0
                    for c in range(1, len(evname) + 1):
                        if evname[c - 1].endswith('TD') and ctype == 'orig':
                            continue
                        count = count + 1
                        if evname[c - 1] in con[2]:
                            val = con[3][con[2].index(evname[c - 1])]
                        else:
                            val = 0.0
                        ev_txt += contrast_element.substitute(cnum=tidx,
                                                              element=count,
                                                              ctype=ctype, val=val)
                        ev_txt += "\n"
                    if con[0] in con_map.keys():
                        for fconidx in con_map[con[0]]:
                            ev_txt += contrast_ftest_element.substitute(cnum=ftest_idx.index(fconidx) + 1,
                                                                        element=tidx,
                                                                        ctype=ctype,
                                                                        val=1)
                        ev_txt += "\n"

            # add contrast mask info
            ev_txt += contrastmask_header.substitute()
            for j, _ in enumerate(contrasts):
                for k, _ in enumerate(contrasts):
                    if j != k:
                        ev_txt += contrastmask_element.substitute(c1=j + 1,
                                                                  c2=k + 1)
            ev_txt += contrastmask_footer.substitute()
        return num_evs, ev_txt

    def _format_session_info(self, session_info):
        if isinstance(session_info, dict):
            session_info = [session_info]
        return session_info

    def _get_func_files(self, session_info):
        """Returns functional files in the order of runs
        """
        func_files = []
        for i, info in enumerate(session_info):
            func_files.insert(i, info['scans'])
        return func_files

    def _run_interface(self, runtime):
        cwd = os.getcwd()
        fsf_header = load_template('feat_header_l1.tcl')
        fsf_postscript = load_template('feat_nongui.tcl')

        prewhiten = 0
        if isdefined(self.inputs.model_serial_correlations):
            prewhiten = int(self.inputs.model_serial_correlations)
        usetd = 0
        no_bases = False
        basis_key = self.inputs.bases.keys()[0]
        if basis_key in ['dgamma', 'gamma']:
            usetd = int(self.inputs.bases[basis_key]['derivs'])
        if basis_key == 'none':
            no_bases = True
        session_info = self._format_session_info(self.inputs.session_info)
        func_files = self._get_func_files(session_info)
        n_tcon = 0
        n_fcon = 0
        if isdefined(self.inputs.contrasts):
            for i, c in enumerate(self.inputs.contrasts):
                if c[1] == 'T':
                    n_tcon += 1
                elif c[1] == 'F':
                    n_fcon += 1

        for i, info in enumerate(session_info):
            do_tempfilter = 1
            if info['hpf'] == np.inf:
                do_tempfilter = 0
            num_evs, cond_txt = self._create_ev_files(cwd, info, i, usetd,
                                                      self.inputs.contrasts,
                                                      no_bases, do_tempfilter)
            nim = load(func_files[i])
            (_, _, _, timepoints) = nim.get_shape()
            fsf_txt = fsf_header.substitute(run_num=i,
                                            interscan_interval=self.inputs.interscan_interval,
                                            num_vols=timepoints,
                                            prewhiten=prewhiten,
                                            num_evs=num_evs[0],
                                            num_evs_real=num_evs[1],
                                            num_tcon=n_tcon,
                                            num_fcon=n_fcon,
                                            high_pass_filter_cutoff=info['hpf'],
                                            temphp_yn=do_tempfilter,
                                            func_file=func_files[i])
            fsf_txt += cond_txt
            fsf_txt += fsf_postscript.substitute(overwrite=1)

            f = open(os.path.join(cwd, 'run%d.fsf' % i), 'w')
            f.write(fsf_txt)
            f.close()

        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        cwd = os.getcwd()
        outputs['fsf_files'] = []
        outputs['ev_files'] = []
        usetd = 0
        basis_key = self.inputs.bases.keys()[0]
        if basis_key in ['dgamma', 'gamma']:
            usetd = int(self.inputs.bases[basis_key]['derivs'])
        for runno, runinfo in enumerate(self._format_session_info(self.inputs.session_info)):
            outputs['fsf_files'].append(os.path.join(cwd, 'run%d.fsf' % runno))
            outputs['ev_files'].insert(runno, [])
            evname = []
            for field in ['cond', 'regress']:
                for i, cond in enumerate(runinfo[field]):
                    name = cond['name']
                    evname.append(name)
                    evfname = os.path.join(cwd, 'ev_%s_%d_%d.txt' % (name, runno,
                                                                     len(evname)))
                    if field == 'cond':
                        if usetd:
                            evname.append(name + 'TD')
                    outputs['ev_files'][runno].append(os.path.join(cwd, evfname))
        return outputs


class FEATInputSpec(FSLCommandInputSpec):
    fsf_file = File(exist=True, mandatory=True, argstr="%s", position=0,
                    desc="File specifying the feat design spec file")


class FEATOutputSpec(TraitedSpec):
    feat_dir = Directory(exists=True)


class FEAT(FSLCommand):
    """Uses FSL feat to calculate first level stats
    """
    _cmd = 'feat'
    input_spec = FEATInputSpec
    output_spec = FEATOutputSpec

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['feat_dir'] = glob(os.path.join(os.getcwd(), '*feat'))[0]
        return outputs


class FEATModelInputSpec(FSLCommandInputSpec):
    fsf_file = File(exist=True, mandatory=True, argstr="%s", position=0,
                    desc="File specifying the feat design spec file",
                    copyfile=False)
    ev_files = traits.List(File(exists=True),
                              mandatory=True, argstr="%s",
                              desc="Event spec files generated by level1design",
                              position=1, copyfile=False)


class FEATModelOutpuSpec(TraitedSpec):
    design_file = File(exists=True, desc='Mat file containing ascii matrix for design')
    design_image = File(exists=True, desc='Graphical representation of design matrix')
    design_cov = File(exists=True, desc='Graphical representation of design covariance')
    con_file = File(exists=True, desc='Contrast file containing contrast vectors')
    fcon_file = File(desc='Contrast file containing contrast vectors')


class FEATModel(FSLCommand):
    """Uses FSL feat_model to generate design.mat files
    """
    _cmd = 'feat_model'
    input_spec = FEATModelInputSpec
    output_spec = FEATModelOutpuSpec

    def _format_arg(self, name, trait_spec, value):
        if name == 'fsf_file':
            return super(FEATModel, self)._format_arg(name, trait_spec, self._get_design_root(value))
        elif name == 'ev_files':
            return ''
        else:
            return super(FEATModel, self)._format_arg(name, trait_spec, value)

    def _get_design_root(self, infile):
        _, fname = os.path.split(infile)
        return fname.split('.')[0]

    def _list_outputs(self):
        #TODO: figure out file names and get rid off the globs
        outputs = self._outputs().get()
        root = self._get_design_root(list_to_filename(self.inputs.fsf_file))
        design_file = glob(os.path.join(os.getcwd(), '%s*.mat' % root))
        assert len(design_file) == 1, 'No mat file generated by FEAT Model'
        outputs['design_file'] = design_file[0]
        design_image = glob(os.path.join(os.getcwd(), '%s.png' % root))
        assert len(design_image) == 1, 'No design image generated by FEAT Model'
        outputs['design_image'] = design_image[0]
        design_cov = glob(os.path.join(os.getcwd(), '%s_cov.png' % root))
        assert len(design_cov) == 1, 'No covariance image generated by FEAT Model'
        outputs['design_cov'] = design_cov[0]
        con_file = glob(os.path.join(os.getcwd(), '%s*.con' % root))
        assert len(con_file) == 1, 'No con file generated by FEAT Model'
        outputs['con_file'] = con_file[0]
        fcon_file = glob(os.path.join(os.getcwd(), '%s*.fts' % root))
        if fcon_file:
            assert len(fcon_file) == 1, 'No fts file generated by FEAT Model'
            outputs['fcon_file'] = fcon_file[0]
        return outputs


# interface to fsl command line model fit routines
# ohinds: 2009-12-28

class FILMGLSInputSpec(FSLCommandInputSpec):
    in_file = File(exists=True, mandatory=True, position=-3,
                  argstr='%s',
                  desc='input data file')
    design_file = File(exists=True, position=-2,
                       argstr='%s',
                       desc='design matrix file')
    threshold = traits.Float(1000, min=0, argstr='%f',
                             position=-1,
                             desc='threshold')
    smooth_autocorr = traits.Bool(argstr='-sa',
                                  desc='Smooth auto corr estimates')
    mask_size = traits.Int(argstr='-ms %d',
                           desc="susan mask size")
    brightness_threshold = traits.Int(min=0, argstr='-epith %d',
        desc='susan brightness threshold, otherwise it is estimated')
    full_data = traits.Bool(argstr='-v', desc='output full data')
    _estimate_xor = ['autocorr_estimate_only', 'fit_armodel', 'tukey_window',
                     'multitaper_product', 'use_pava', 'autocorr_noestimate']
    autocorr_estimate_only = traits.Bool(argstr='-ac',
                                    xor=_estimate_xor,
                   desc='perform autocorrelation estimatation only')
    fit_armodel = traits.Bool(argstr='-ar', xor=_estimate_xor,
        desc='fits autoregressive model - default is to use tukey with M=sqrt(numvols)')
    tukey_window = traits.Int(argstr='-tukey %d', xor=_estimate_xor,
        desc='tukey window size to estimate autocorr')
    multitaper_product = traits.Int(argstr='-mt %d', xor=_estimate_xor,
               desc='multitapering with slepian tapers and num is the time-bandwidth product')
    use_pava = traits.Bool(argstr='-pava', desc='estimates autocorr using PAVA')
    autocorr_noestimate = traits.Bool(argstr='-noest', xor=_estimate_xor,
                   desc='do not estimate autocorrs')
    output_pwdata = traits.Bool(argstr='-output_pwdata',
                   desc='output prewhitened data and average design matrix')
    results_dir = Directory('results', argstr='-rn %s', usedefault=True,
                            desc='directory to store results in')


class FILMGLSOutputSpec(TraitedSpec):
    param_estimates = OutputMultiPath(File(exists=True),
          desc='Parameter estimates for each column of the design matrix')
    residual4d = File(exists=True,
          desc='Model fit residual mean-squared error for each time point')
    dof_file = File(exists=True, desc='degrees of freedom')
    sigmasquareds = File(exists=True, desc='summary of residuals, See Woolrich, et. al., 2001')
    results_dir = Directory(exists=True,
                         desc='directory storing model estimation output')
    corrections = File(exists=True,
                       desc='statistical corrections used within FILM modelling')
    logfile = File(exists=True,
                   desc='FILM run logfile')


class FILMGLS(FSLCommand):
    """Use FSL film_gls command to fit a design matrix to voxel timeseries

    Examples
    --------

    Initialize with no options, assigning them when calling run:

    >>> from nipype.interfaces import fsl
    >>> fgls = fsl.FILMGLS()
    >>> res = fgls.run('in_file', 'design_file', 'thresh', rn='stats') #doctest: +SKIP

    Assign options through the ``inputs`` attribute:

    >>> fgls = fsl.FILMGLS()
    >>> fgls.inputs.in_file = 'functional.nii'
    >>> fgls.inputs.design_file = 'design.mat'
    >>> fgls.inputs.threshold = 10
    >>> fgls.inputs.results_dir = 'stats'
    >>> res = fgls.run() #doctest: +SKIP

    Specify options when creating an instance:

    >>> fgls = fsl.FILMGLS(in_file='functional.nii', \
design_file='design.mat', \
threshold=10, results_dir='stats')
    >>> res = fgls.run() #doctest: +SKIP

    """

    _cmd = 'film_gls'
    input_spec = FILMGLSInputSpec
    output_spec = FILMGLSOutputSpec

    def _get_pe_files(self, cwd):
        files = None
        if isdefined(self.inputs.design_file):
            fp = open(self.inputs.design_file, 'rt')
            for line in fp.readlines():
                if line.startswith('/NumWaves'):
                    numpes = int(line.split()[-1])
                    files = []
                    for i in range(numpes):
                        files.append(self._gen_fname('pe%d.nii' % (i + 1),
                                                     cwd=cwd))
                    break
            fp.close()
        return files

    def _list_outputs(self):
        outputs = self._outputs().get()
        cwd = os.getcwd()
        results_dir = os.path.join(cwd, self.inputs.results_dir)
        outputs['results_dir'] = results_dir
        pe_files = self._get_pe_files(results_dir)
        if pe_files:
            outputs['param_estimates'] = pe_files
        outputs['residual4d'] = self._gen_fname('res4d.nii', cwd=results_dir)
        outputs['dof_file'] = os.path.join(results_dir, 'dof')
        outputs['sigmasquareds'] = self._gen_fname('sigmasquareds.nii',
                                                   cwd=results_dir)
        outputs['corrections'] = self._gen_fname('corrections.nii',
                                                 cwd=results_dir)
        outputs['logfile'] = self._gen_fname('logfile',
                                             change_ext=False,
                                             cwd=results_dir)
        return outputs


class FEATRegisterInputSpec(BaseInterfaceInputSpec):
    feat_dirs = InputMultiPath(Directory(), exist=True, desc="Lower level feat dirs",
                               mandatory=True)
    reg_image = File(exist=True, desc="image to register to (will be treated as standard)",
                     mandatory=True)
    reg_dof = traits.Int(12, desc="registration degrees of freedom", usedefault=True)


class FEATRegisterOutputSpec(TraitedSpec):
    fsf_file = File(exists=True,
                                desc="FSL feat specification file")


class FEATRegister(BaseInterface):
    """Register feat directories to a specific standard
    """
    input_spec = FEATRegisterInputSpec
    output_spec = FEATRegisterOutputSpec

    def _run_interface(self, runtime):
        fsf_header = load_template('featreg_header.tcl')
        fsf_footer = load_template('feat_nongui.tcl')
        fsf_dirs = load_template('feat_fe_featdirs.tcl')

        num_runs = len(self.inputs.feat_dirs)
        fsf_txt = fsf_header.substitute(num_runs=num_runs,
                                        regimage=self.inputs.reg_image,
                                        regdof=self.inputs.reg_dof)
        for i, rundir in enumerate(filename_to_list(self.inputs.feat_dirs)):
            fsf_txt += fsf_dirs.substitute(runno=i + 1,
                                           rundir=os.path.abspath(rundir))
        fsf_txt += fsf_footer.substitute()
        f = open(os.path.join(os.getcwd(), 'register.fsf'), 'wt')
        f.write(fsf_txt)
        f.close()

        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['fsf_file'] = os.path.abspath(os.path.join(os.getcwd(), 'register.fsf'))
        return outputs


class FLAMEOInputSpec(FSLCommandInputSpec):
    cope_file = File(exists=True, argstr='--copefile=%s', mandatory=True,
                    desc='cope regressor data file')
    var_cope_file = File(exists=True, argstr='--varcopefile=%s',
                       desc='varcope weightings data file')
    dof_var_cope_file = File(exists=True, argstr='--dofvarcopefile=%s',
                          desc='dof data file for varcope data')
    mask_file = File(exists=True, argstr='--maskfile=%s', mandatory=True,
                    desc='mask file')
    design_file = File(exists=True, argstr='--designfile=%s', mandatory=True,
                      desc='design matrix file')
    t_con_file = File(exists=True, argstr='--tcontrastsfile=%s', mandatory=True,
                    desc='ascii matrix specifying t-contrasts')
    f_con_file = File(exists=True, argstr='--fcontrastsfile=%s',
                    desc='ascii matrix specifying f-contrasts')
    cov_split_file = File(exists=True, argstr='--covsplitfile=%s', mandatory=True,
                        desc='ascii matrix specifying the groups the covariance is split into')
    run_mode = traits.Enum('fe', 'ols', 'flame1', 'flame12', argstr='--runmode=%s',
                          mandatory=True, desc='inference to perform')
    n_jumps = traits.Int(argstr='--njumps=%d', desc='number of jumps made by mcmc')
    burnin = traits.Int(argstr='--burnin=%d',
                        desc='number of jumps at start of mcmc to be discarded')
    sample_every = traits.Int(argstr='--sampleevery=%d',
                             desc='number of jumps for each sample')
    fix_mean = traits.Bool(argstr='--fixmean', desc='fix mean for tfit')
    infer_outliers = traits.Bool(argstr='--inferoutliers',
                                desc='infer outliers - not for fe')
    no_pe_outputs = traits.Bool(argstr='--nopeoutput',
                             desc='do not output pe files')
    sigma_dofs = traits.Int(argstr='--sigma_dofs=%d',
                            desc='sigma (in mm) to use for Gaussian smoothing the DOFs in FLAME 2. Default is 1mm, -1 indicates no smoothing')
    outlier_iter = traits.Int(argstr='--ioni=%d',
                              desc='Number of max iterations to use when inferring outliers. Default is 12.')
    log_dir = Directory("stats", argstr='--ld=%s', usedefault=True)  # ohinds
    # no support for ven, vef


class FLAMEOOutputSpec(TraitedSpec):
    pes = OutputMultiPath(exists=True, desc="Parameter estimates for each column of the design matrix" +
                "for each voxel")
    res4d = OutputMultiPath(exists=True, desc="Model fit residual mean-squared error for each time point")
    copes = OutputMultiPath(exists=True, desc="Contrast estimates for each contrast")
    var_copes = OutputMultiPath(exists=True, desc="Variance estimates for each contrast")
    zstats = OutputMultiPath(exists=True, desc="z-stat file for each contrast")
    tstats = OutputMultiPath(exists=True, desc="t-stat file for each contrast")
    mrefvars = OutputMultiPath(exists=True, desc="mean random effect variances for each contrast")
    tdof = OutputMultiPath(exists=True, desc="temporal dof file for each contrast")
    weights = OutputMultiPath(exists=True, desc="weights file for each contrast")
    stats_dir = Directory(exists=True, desc="directory storing model estimation output")


# interface to fsl command line higher level model fit
# satra: 2010-01-09

class FLAMEO(FSLCommand):
    """Use FSL flameo command to perform higher level model fits

    Examples
    --------

    Initialize FLAMEO with no options, assigning them when calling run:

    >>> from nipype.interfaces import fsl
    >>> import os
    >>> flameo = fsl.FLAMEO(cope_file='cope.nii.gz', \
                            var_cope_file='varcope.nii.gz', \
                            cov_split_file='cov_split.mat', \
                            design_file='design.mat', \
                            t_con_file='design.con', \
                            mask_file='mask.nii', \
                            run_mode='fe')
    >>> flameo.cmdline
    'flameo --copefile=cope.nii.gz --covsplitfile=cov_split.mat --designfile=design.mat --ld=stats --maskfile=mask.nii --runmode=fe --tcontrastsfile=design.con --varcopefile=varcope.nii.gz'

    """

    _cmd = 'flameo'
    input_spec = FLAMEOInputSpec
    output_spec = FLAMEOOutputSpec

    # ohinds: 2010-04-06
    def _run_interface(self, runtime):
        log_dir = self.inputs.log_dir
        cwd = os.getcwd()
        if os.access(os.path.join(cwd, log_dir), os.F_OK):
            rmtree(os.path.join(cwd, log_dir))

        return super(FLAMEO, self)._run_interface(runtime)

    # ohinds: 2010-04-06
    # made these compatible with flameo
    def _list_outputs(self):
        outputs = self._outputs().get()
        pth = os.path.join(os.getcwd(), self.inputs.log_dir)

        pes = glob(os.path.join(pth, 'pe[0-9]*.*'))
        assert len(pes) >= 1, 'No pe volumes generated by FSL Estimate'
        outputs['pes'] = pes

        res4d = glob(os.path.join(pth, 'res4d.*'))
        assert len(res4d) == 1, 'No residual volume generated by FSL Estimate'
        outputs['res4d'] = res4d[0]

        copes = glob(os.path.join(pth, 'cope[0-9]*.*'))
        assert len(copes) >= 1, 'No cope volumes generated by FSL CEstimate'
        outputs['copes'] = copes

        var_copes = glob(os.path.join(pth, 'varcope[0-9]*.*'))
        assert len(var_copes) >= 1, 'No varcope volumes generated by FSL CEstimate'
        outputs['var_copes'] = var_copes

        zstats = glob(os.path.join(pth, 'zstat[0-9]*.*'))
        assert len(zstats) >= 1, 'No zstat volumes generated by FSL CEstimate'
        outputs['zstats'] = zstats

        tstats = glob(os.path.join(pth, 'tstat[0-9]*.*'))
        assert len(tstats) >= 1, 'No tstat volumes generated by FSL CEstimate'
        outputs['tstats'] = tstats

        mrefs = glob(os.path.join(pth, 'mean_random_effects_var[0-9]*.*'))
        assert len(mrefs) >= 1, 'No mean random effects volumes generated by FLAMEO'
        outputs['mrefvars'] = mrefs

        tdof = glob(os.path.join(pth, 'tdof_t[0-9]*.*'))
        assert len(tdof) >= 1, 'No T dof volumes generated by FLAMEO'
        outputs['tdof'] = tdof

        weights = glob(os.path.join(pth, 'weights[0-9]*.*'))
        assert len(weights) >= 1, 'No weight volumes generated by FLAMEO'
        outputs['weights'] = weights

        outputs['stats_dir'] = pth

        return outputs


class ContrastMgrInputSpec(FSLCommandInputSpec):
    tcon_file = File(exists=True, mandatory=True,
                     argstr='%s', position=-1,
                     desc='contrast file containing T-contrasts')
    fcon_file = File(exists=True, argstr='-f %s',
                     desc='contrast file containing F-contrasts')
    param_estimates = InputMultiPath(File(exists=True),
                                     argstr='', copyfile=False,
                                     mandatory=True,
          desc='Parameter estimates for each column of the design matrix')
    corrections = File(exists=True, copyfile=False, mandatory=True,
                       desc='statistical corrections used within FILM modelling')
    dof_file = File(exists=True, argstr='', copyfile=False, mandatory=True,
                    desc='degrees of freedom')
    sigmasquareds = File(exists=True, argstr='', position=-2,
                         copyfile=False, mandatory=True,
                         desc='summary of residuals, See Woolrich, et. al., 2001')
    contrast_num = traits.Int(min=1, argstr='-cope',
                desc='contrast number to start labeling copes from')
    suffix = traits.Str(argstr='-suffix %s',
                        desc='suffix to put on the end of the cope filename before the contrast number, default is nothing')


class ContrastMgrOutputSpec(TraitedSpec):
    copes = OutputMultiPath(File(exists=True),
                                 desc='Contrast estimates for each contrast')
    varcopes = OutputMultiPath(File(exists=True),
                                 desc='Variance estimates for each contrast')
    zstats = OutputMultiPath(File(exists=True),
                                 desc='z-stat file for each contrast')
    tstats = OutputMultiPath(File(exists=True),
                                 desc='t-stat file for each contrast')
    fstats = OutputMultiPath(File(exists=True),
                                 desc='f-stat file for each contrast')
    zfstats = OutputMultiPath(File(exists=True),
                              desc='z-stat file for each F contrast')
    neffs = OutputMultiPath(File(exists=True),
                                 desc='neff file ?? for each contrast')


class ContrastMgr(FSLCommand):
    """Use FSL contrast_mgr command to evaluate contrasts

    In interface mode this file assumes that all the required inputs are in the
    same location.
    """

    _cmd = 'contrast_mgr'
    input_spec = ContrastMgrInputSpec
    output_spec = ContrastMgrOutputSpec

    def _run_interface(self, runtime):
        # The returncode is meaningless in ContrastMgr.  So check the output
        # in stderr and if it's set, then update the returncode
        # accordingly.
        runtime = super(ContrastMgr, self)._run_interface(runtime)
        if runtime.stderr:
            self.raise_exception(runtime)
        return runtime

    def _format_arg(self, name, trait_spec, value):
        if name in ['param_estimates', 'corrections', 'dof_file']:
            return ''
        elif name in ['sigmasquareds']:
            path, _ = os.path.split(value)
            return path
        else:
            return super(ContrastMgr, self)._format_arg(name, trait_spec, value)

    def _get_design_root(self, infile):
        _, fname = os.path.split(infile)
        return fname.split('.')[0]

    def _get_numcons(self):
        numtcons = 0
        numfcons = 0
        if isdefined(self.inputs.tcon_file):
            fp = open(self.inputs.tcon_file, 'rt')
            for line in fp.readlines():
                if line.startswith('/NumContrasts'):
                    numtcons = int(line.split()[-1])
                    break
            fp.close()
        if isdefined(self.inputs.fcon_file):
            fp = open(self.inputs.fcon_file, 'rt')
            for line in fp.readlines():
                if line.startswith('/NumContrasts'):
                    numfcons = int(line.split()[-1])
                    break
            fp.close()
        return numtcons, numfcons

    def _list_outputs(self):
        outputs = self._outputs().get()
        pth, _ = os.path.split(self.inputs.sigmasquareds)
        numtcons, numfcons = self._get_numcons()
        base_contrast = 1
        if isdefined(self.inputs.contrast_num):
            base_contrast = self.inputs.contrast_num
        copes = []
        varcopes = []
        zstats = []
        tstats = []
        neffs = []
        for i in range(numtcons):
            copes.append(self._gen_fname('cope%d.nii' % (base_contrast + i),
                                       cwd=pth))
            varcopes.append(self._gen_fname('varcope%d.nii' % (base_contrast + i),
                                          cwd=pth))
            zstats.append(self._gen_fname('zstat%d.nii' % (base_contrast + i),
                                        cwd=pth))
            tstats.append(self._gen_fname('tstat%d.nii' % (base_contrast + i),
                                        cwd=pth))
            neffs.append(self._gen_fname('neff%d.nii' % (base_contrast + i),
                                       cwd=pth))
        if copes:
            outputs['copes'] = copes
            outputs['varcopes'] = varcopes
            outputs['zstats'] = zstats
            outputs['tstats'] = tstats
            outputs['neffs'] = neffs
        fstats = []
        zfstats = []
        for i in range(numfcons):
            fstats.append(self._gen_fname('fstat%d.nii' % (base_contrast + i),
                                        cwd=pth))
            zfstats.append(self._gen_fname('zfstat%d.nii' % (base_contrast + i),
                           cwd=pth))
        if fstats:
            outputs['fstats'] = fstats
            outputs['zfstats'] = zfstats
        return outputs


class L2ModelInputSpec(BaseInterfaceInputSpec):
    num_copes = traits.Int(min=1, mandatory=True,
                             desc='number of copes to be combined')


class L2ModelOutputSpec(TraitedSpec):
    design_mat = File(exists=True, desc='design matrix file')
    design_con = File(exists=True, desc='design contrast file')
    design_grp = File(exists=True, desc='design group file')


class L2Model(BaseInterface):
    """Generate subject specific second level model

    Examples
    --------

    >>> from nipype.interfaces.fsl import L2Model
    >>> model = L2Model(num_copes=3) # 3 sessions

    """

    input_spec = L2ModelInputSpec
    output_spec = L2ModelOutputSpec

    def _run_interface(self, runtime):
        cwd = os.getcwd()
        mat_txt = ['/NumWaves       1',
                   '/NumPoints      %d' % self.inputs.num_copes,
                   '/PPheights      %e' % 1,
                   '',
                   '/Matrix']
        for i in range(self.inputs.num_copes):
            mat_txt += ['%e' % 1]
        mat_txt = '\n'.join(mat_txt)

        con_txt = ['/ContrastName1   group mean',
                   '/NumWaves       1',
                   '/NumContrasts   1',
                   '/PPheights          %e' % 1,
                   '/RequiredEffect     100.0',  # XX where does this
                   #number come from
                   '',
                   '/Matrix',
                   '%e' % 1]
        con_txt = '\n'.join(con_txt)

        grp_txt = ['/NumWaves       1',
                   '/NumPoints      %d' % self.inputs.num_copes,
                   '',
                   '/Matrix']
        for i in range(self.inputs.num_copes):
            grp_txt += ['1']
        grp_txt = '\n'.join(grp_txt)

        txt = {'design.mat': mat_txt,
               'design.con': con_txt,
               'design.grp': grp_txt}

        # write design files
        for i, name in enumerate(['design.mat', 'design.con', 'design.grp']):
            f = open(os.path.join(cwd, name), 'wt')
            f.write(txt[name])
            f.close()

        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        for field in outputs.keys():
            outputs[field] = os.path.join(os.getcwd(),
                                          field.replace('_', '.'))
        return outputs


class MultipleRegressDesignInputSpec(BaseInterfaceInputSpec):
    contrasts = traits.List(
        traits.Either(traits.Tuple(traits.Str,
                                   traits.Enum('T'),
                                   traits.List(traits.Str),
                                   traits.List(traits.Float)),
                      traits.Tuple(traits.Str,
                                   traits.Enum('F'),
                                   traits.List(traits.Tuple(traits.Str,
                                                            traits.Enum('T'),
                                                            traits.List(traits.Str),
                                                            traits.List(traits.Float)),
                                               ))),
        mandatory=True,
        desc="List of contrasts with each contrast being a list of the form - \
[('name', 'stat', [condition list], [weight list])]. if \
session list is None or not provided, all sessions are used. For F \
contrasts, the condition list should contain previously defined \
T-contrasts without any weight list.")
    regressors = traits.Dict(traits.Str, traits.List(traits.Float),
                             mandatory=True,
                             desc='dictionary containing named lists of regressors')
    groups = traits.List(traits.Int,
                  desc='list of group identifiers (defaults to single group)')


class MultipleRegressDesignOutputSpec(TraitedSpec):
    design_mat = File(exists=True, desc='design matrix file')
    design_con = File(exists=True, desc='design t-contrast file')
    design_fts = File(exists=True, desc='design f-contrast file')
    design_grp = File(exists=True, desc='design group file')


class MultipleRegressDesign(BaseInterface):
    """Generate multiple regression design

    .. note::
      FSL does not demean columns for higher level analysis.

    Please see `FSL documentation <http://www.fmrib.ox.ac.uk/fsl/feat5/detail.html#higher>`_
    for more details on model specification for higher level analysis.

    Examples
    --------

    >>> from nipype.interfaces.fsl import MultipleRegressDesign
    >>> model = MultipleRegressDesign()
    >>> model.inputs.contrasts = [['group mean', 'T',['reg1'],[1]]]
    >>> model.inputs.regressors = dict(reg1=[1, 1, 1], reg2=[2.,-4, 3])
    >>> model.run() # doctest: +SKIP

    """

    input_spec = MultipleRegressDesignInputSpec
    output_spec = MultipleRegressDesignOutputSpec

    def _run_interface(self, runtime):
        cwd = os.getcwd()
        regs = sorted(self.inputs.regressors.keys())
        nwaves = len(regs)
        npoints = len(self.inputs.regressors[regs[0]])
        ntcons = sum([1 for con in self.inputs.contrasts if con[1] == 'T'])
        nfcons = sum([1 for con in self.inputs.contrasts if con[1] == 'F'])
        # write mat file
        mat_txt = ['/NumWaves       %d' % nwaves,
                   '/NumPoints      %d' % npoints]
        ppheights = []
        for reg in regs:
            maxreg = np.max(self.inputs.regressors[reg])
            minreg = np.min(self.inputs.regressors[reg])
            if np.sign(maxreg) == np.sign(minreg):
                regheight = max([abs(minreg), abs(maxreg)])
            else:
                regheight = abs(maxreg - minreg)
            ppheights.append('%e' % regheight)
        mat_txt += ['/PPheights      ' + ' '.join(ppheights)]
        mat_txt += ['',
                    '/Matrix']
        for cidx in range(npoints):
            mat_txt.append(' '.join(['%e' % self.inputs.regressors[key][cidx] for key in regs]))
        mat_txt = '\n'.join(mat_txt)
        # write t-con file
        con_txt = []
        counter = 0
        tconmap = {}
        for conidx, con in enumerate(self.inputs.contrasts):
            if con[1] == 'T':
                tconmap[conidx] = counter
                counter += 1
                con_txt += ['/ContrastName%d   %s' % (counter, con[0])]
        con_txt += ['/NumWaves       %d' % nwaves,
                    '/NumContrasts   %d' % ntcons,
                    '/PPheights          %s' % ' '.join(['%e' % 1 for i in range(counter)]),
                    '/RequiredEffect     %s' % ' '.join(['%.3f' % 100 for i in range(counter)]),
                    '',
                    '/Matrix']
        for idx in sorted(tconmap.keys()):
            convals = np.zeros((nwaves, 1))
            for regidx, reg in enumerate(self.inputs.contrasts[idx][2]):
                convals[regs.index(reg)] = self.inputs.contrasts[idx][3][regidx]
            con_txt.append(' '.join(['%e' % val for val in convals]))
        con_txt = '\n'.join(con_txt)
        # write f-con file
        fcon_txt = ''
        if nfcons:
            fcon_txt = ['/NumWaves       %d' % ntcons,
                        '/NumContrasts   %d' % nfcons,
                        '',
                        '/Matrix']
            for conidx, con in enumerate(self.inputs.contrasts):
                if con[1] == 'F':
                    convals = np.zeros((ntcons, 1))
                    for tcon in con[2]:
                        convals[tconmap[self.inputs.contrasts.index(tcon)]] = 1
                    fcon_txt.append(' '.join(['%d' % val for val in convals]))
                    fcon_txt = '\n'.join(fcon_txt)
        # write group file
        grp_txt = ['/NumWaves       1',
                   '/NumPoints      %d' % npoints,
                   '',
                   '/Matrix']
        for i in range(npoints):
            if isdefined(self.inputs.groups):
                grp_txt += ['%d' % self.inputs.groups[i]]
            else:
                grp_txt += ['1']
        grp_txt = '\n'.join(grp_txt)

        txt = {'design.mat': mat_txt,
               'design.con': con_txt,
               'design.fts': fcon_txt,
               'design.grp': grp_txt}

        # write design files
        for key, val in txt.items():
            if ('fts' in key) and (nfcons == 0):
                continue
            filename = key.replace('_', '.')
            f = open(os.path.join(cwd, filename), 'wt')
            f.write(val)
            f.close()

        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        nfcons = sum([1 for con in self.inputs.contrasts if con[1] == 'F'])
        for field in outputs.keys():
            if ('fts' in field) and (nfcons == 0):
                continue
            outputs[field] = os.path.join(os.getcwd(),
                                          field.replace('_', '.'))
        return outputs


class SMMInputSpec(FSLCommandInputSpec):
    spatial_data_file = File(exists=True, position=0, argstr='--sdf="%s"', mandatory=True,
                           desc="statistics spatial map", copyfile=False)
    mask = File(exist=True, position=1, argstr='--mask="%s"', mandatory=True,
                desc="mask file", copyfile=False)
    no_deactivation_class = traits.Bool(position=2, argstr="--zfstatmode",
                             desc="enforces no deactivation class")


class SMMOutputSpec(TraitedSpec):
    null_p_map = File(exists=True)
    activation_p_map = File(exists=True)
    deactivation_p_map = File(exists=True)


class SMM(FSLCommand):
    '''
    Spatial Mixture Modelling. For more detail on the spatial mixture modelling see
    Mixture Models with Adaptive Spatial Regularisation for Segmentation with an Application to FMRI Data;
    Woolrich, M., Behrens, T., Beckmann, C., and Smith, S.; IEEE Trans. Medical Imaging, 24(1):1-11, 2005.
    '''
    _cmd = 'mm --ld=logdir'
    input_spec = SMMInputSpec
    output_spec = SMMOutputSpec

    def _list_outputs(self):
        outputs = self._outputs().get()
        #TODO get the true logdir from the stdout
        outputs['null_p_map'] = self._gen_fname(basename="w1_mean", cwd="logdir")
        outputs['activation_p_map'] = self._gen_fname(basename="w2_mean", cwd="logdir")
        if not isdefined(self.inputs.no_deactivation_class) or not self.inputs.no_deactivation_class:
            outputs['deactivation_p_map'] = self._gen_fname(basename="w3_mean", cwd="logdir")
        return outputs


class MELODICInputSpec(FSLCommandInputSpec):
    in_files = InputMultiPath(File(exists=True), argstr="-i %s", mandatory=True, position=0,
                              desc="input file names (either single file name or a list)")
    out_dir = Directory(argstr="-o %s", desc="output directory name", genfile=True)
    mask = File(exists=True, argstr="-m %s", desc="file name of mask for thresholding")
    no_mask = traits.Bool(argstr="--nomask", desc="switch off masking")
    update_mask = traits.Bool(argstr="--update_mask", desc="switch off mask updating")
    no_bet = traits.Bool(argstr="--nobet", desc="switch off BET")
    bg_threshold = traits.Float(argstr="--bgthreshold=%f", desc="brain/non-brain threshold used to mask non-brain voxels, as a percentage (only if --nobet selected)")
    dim = traits.Int(argstr="-d %d", desc="dimensionality reduction into #num dimensions"\
                     "(default: automatic estimation)")
    dim_est = traits.Str(argstr="--dimest=%s", desc="use specific dim. estimation technique:"\
                         " lap, bic, mdl, aic, mean (default: lap)")
    sep_whiten = traits.Bool(argstr="--sep_whiten", desc="switch on separate whitening")
    sep_vn = traits.Bool(argstr="--sep_vn", desc="switch off joined variance normalization")
    num_ICs = traits.Int(argstr="-n %d", desc="number of IC's to extract (for deflation approach)")
    approach = traits.Str(argstr="-a %s", desc="approach for decomposition, 2D: defl, symm (default), "\
                          " 3D: tica (default), concat")
    non_linearity = traits.Str(argstr="--nl=%s", desc="nonlinearity: gauss, tanh, pow3, pow4")
    var_norm = traits.Bool(argstr="--vn", desc="switch off variance normalization")
    pbsc = traits.Bool(argstr="--pbsc", desc="switch off conversion to percent BOLD signal change")
    cov_weight = traits.Float(argstr="--covarweight=%f", desc="voxel-wise weights for the covariance "\
                              "matrix (e.g. segmentation information)")
    epsilon = traits.Float(argstr="--eps=%f", desc="minimum error change")
    epsilonS = traits.Float(argstr="--epsS=%f", desc="minimum error change for rank-1 approximation in TICA")
    maxit = traits.Int(argstr="--maxit=%d", desc="maximum number of iterations before restart")
    max_restart = traits.Int(argstr="--maxrestart=%d", desc="maximum number of restarts")
    mm_thresh = traits.Float(argstr="--mmthresh=%f", desc="threshold for Mixture Model based inference")
    no_mm = traits.Bool(argstr="--no_mm", desc="switch off mixture modelling on IC maps")
    ICs = File(exists=True, argstr="--ICs=%s", desc="filename of the IC components file for mixture modelling")
    mix = File(exists=True, argstr="--mix=%s", desc="mixing matrix for mixture modelling / filtering")
    smode = File(exists=True, argstr="--smode=%s", desc="matrix of session modes for report generation")
    rem_cmp = traits.List(traits.Int, argstr="-f %d", desc="component numbers to remove")
    report = traits.Bool(argstr="--report", desc="generate Melodic web report")
    bg_image = File(exists=True, argstr="--bgimage=%s", desc="specify background image for report"\
                    " (default: mean image)")
    tr_sec = traits.Float(argstr="--tr=%f", desc="TR in seconds")
    log_power = traits.Bool(argstr="--logPower", desc="calculate log of power for frequency spectrum")
    t_des = File(exists=True, argstr="--Tdes=%s", desc="design matrix across time-domain")
    t_con = File(exists=True, argstr="--Tcon=%s", desc="t-contrast matrix across time-domain")
    s_des = File(exists=True, argstr="--Sdes=%s", desc="design matrix across subject-domain")
    s_con = File(exists=True, argstr="--Scon=%s", desc="t-contrast matrix across subject-domain")
    out_all = traits.Bool(argstr="--Oall", desc="output everything")
    out_unmix = traits.Bool(argstr="--Ounmix", desc="output unmixing matrix")
    out_stats = traits.Bool(argstr="--Ostats", desc="output thresholded maps and probability maps")
    out_pca = traits.Bool(argstr="--Opca", desc="output PCA results")
    out_white = traits.Bool(argstr="--Owhite", desc="output whitening/dewhitening matrices")
    out_orig = traits.Bool(argstr="--Oorig", desc="output the original ICs")
    out_mean = traits.Bool(argstr="--Omean", desc="output mean volume")
    report_maps = traits.Str(argstr="--report_maps=%s", desc="control string for spatial map images (see slicer)")
    remove_deriv = traits.Bool(argstr="--remove_deriv", desc="removes every second entry in paradigm"\
                               " file (EV derivatives)")


class MELODICOutputSpec(TraitedSpec):
    out_dir = Directory(exists=True)
    report_dir = Directory(exists=True)


class MELODIC(FSLCommand):
    """Multivariate Exploratory Linear Optimised Decomposition into Independent Components

    Examples
    --------

    >>> melodic_setup = MELODIC()
    >>> melodic_setup.inputs.approach = 'tica'
    >>> melodic_setup.inputs.in_files = ['functional.nii', 'functional2.nii', 'functional3.nii']
    >>> melodic_setup.inputs.no_bet = True
    >>> melodic_setup.inputs.bg_threshold = 10
    >>> melodic_setup.inputs.tr_sec = 1.5
    >>> melodic_setup.inputs.mm_thresh = 0.5
    >>> melodic_setup.inputs.out_stats = True
    >>> melodic_setup.inputs.t_des = 'timeDesign.mat'
    >>> melodic_setup.inputs.t_con = 'timeDesign.con'
    >>> melodic_setup.inputs.s_des = 'subjectDesign.mat'
    >>> melodic_setup.inputs.s_con = 'subjectDesign.con'
    >>> melodic_setup.inputs.out_dir = 'groupICA.out'
    >>> melodic_setup.run() # doctest: +SKIP


    """
    input_spec = MELODICInputSpec
    output_spec = MELODICOutputSpec
    _cmd = 'melodic'

    def _list_outputs(self):
        outputs = self.output_spec().get()
        outputs['out_dir'] = self.inputs.out_dir
        if not isdefined(outputs['out_dir']):
            outputs['out_dir'] = self._gen_filename("out_dir")
        if isdefined(self.inputs.report) and self.inputs.report:
            outputs['report_dir'] = os.path.join(self._gen_filename("out_dir"), "report")
        return outputs

    def _gen_filename(self, name):
        if name == "out_dir":
            return os.getcwd()


class SmoothEstimateInputSpec(FSLCommandInputSpec):
    dof = traits.Int(argstr='--dof=%d', mandatory=True,
                     xor=['zstat_file'],
                     desc='number of degrees of freedom')
    mask_file = File(argstr='--mask=%s',
                     exists=True, mandatory=True,
                     desc='brain mask volume')
    residual_fit_file = File(argstr='--res=%s',
                             exists=True, requires=['dof'],
                             desc='residual-fit image file')
    zstat_file = File(argstr='--zstat=%s',
                      exists=True, xor=['dof'],
                      desc='zstat image file')


class SmoothEstimateOutputSpec(TraitedSpec):
    dlh = traits.Float(desc='smoothness estimate sqrt(det(Lambda))')
    volume = traits.Int(desc='number of voxels in mask')
    resels = traits.Float(desc='number of resels')


class SmoothEstimate(FSLCommand):
    """ Estimates the smoothness of an image

    Examples
    --------

    >>> est = SmoothEstimate()
    >>> est.inputs.zstat_file = 'zstat1.nii.gz'
    >>> est.inputs.mask_file = 'mask.nii'
    >>> est.cmdline
    'smoothest --mask=mask.nii --zstat=zstat1.nii.gz'

    """

    input_spec = SmoothEstimateInputSpec
    output_spec = SmoothEstimateOutputSpec
    _cmd = 'smoothest'

    def aggregate_outputs(self, runtime=None, needed_outputs=None):
        outputs = self._outputs()
        stdout = runtime.stdout.split('\n')
        outputs.dlh = float(stdout[0].split()[1])
        outputs.volume = int(stdout[1].split()[1])
        outputs.resels = float(stdout[2].split()[1])
        return outputs


class ClusterInputSpec(FSLCommandInputSpec):
    in_file = File(argstr='--in=%s', mandatory=True,
                   exists=True, desc='input volume')
    threshold = traits.Float(argstr='--thresh=%.10f',
                             mandatory=True,
                             desc='threshold for input volume')
    out_index_file = traits.Either(traits.Bool, File,
                                   argstr='--oindex=%s',
                                   desc='output of cluster index (in size order)', hash_files=False)
    out_threshold_file = traits.Either(traits.Bool, File,
                                       argstr='--othresh=%s',
                                       desc='thresholded image', hash_files=False)
    out_localmax_txt_file = traits.Either(traits.Bool, File,
                                          argstr='--olmax=%s',
                                          desc='local maxima text file', hash_files=False)
    out_localmax_vol_file = traits.Either(traits.Bool, File,
                                          argstr='--olmaxim=%s',
                                          desc='output of local maxima volume', hash_files=False)
    out_size_file = traits.Either(traits.Bool, File,
                                  argstr='--osize=%s',
                                  desc='filename for output of size image', hash_files=False)
    out_max_file = traits.Either(traits.Bool, File,
                                 argstr='--omax=%s',
                                 desc='filename for output of max image', hash_files=False)
    out_mean_file = traits.Either(traits.Bool, File,
                                  argstr='--omean=%s',
                                  desc='filename for output of mean image', hash_files=False)
    out_pval_file = traits.Either(traits.Bool, File,
                                  argstr='--opvals=%s',
                                  desc='filename for image output of log pvals', hash_files=False)
    pthreshold = traits.Float(argstr='--pthresh=%.10f',
                              requires=['dlh', 'volume'],
                              desc='p-threshold for clusters')
    peak_distance = traits.Float(argstr='--peakdist=%.10f',
                                 desc='minimum distance between local maxima/minima, in mm (default 0)')
    cope_file = traits.File(argstr='--cope=%s',
                            desc='cope volume')
    volume = traits.Int(argstr='--volume=%d',
                        desc='number of voxels in the mask')
    dlh = traits.Float(argstr='--dlh=%.10f',
                       desc='smoothness estimate = sqrt(det(Lambda))')
    fractional = traits.Bool('--fractional',
                             desc='interprets the threshold as a fraction of the robust range')
    connectivity = traits.Int(argstr='--connectivity=%d',
                              desc='the connectivity of voxels (default 26)')
    use_mm = traits.Bool('--mm', desc='use mm, not voxel, coordinates')
    find_min = traits.Bool('--min', desc='find minima instead of maxima')
    no_table = traits.Bool('--no_table', desc='suppresses printing of the table info')
    minclustersize = traits.Bool(argstr='--minclustersize',
                                 desc='prints out minimum significant cluster size')
    xfm_file = File(argstr='--xfm=%s',
                    desc='filename for Linear: input->standard-space transform. Non-linear: input->highres transform')
    std_space_file = File(argstr='--stdvol=%s',
                          desc='filename for standard-space volume')
    num_maxima = traits.Int(argstr='--num=%d',
                            desc='no of local maxima to report')
    warpfield_file = File(argstr='--warpvol=%s',
                          desc='file contining warpfield')


class ClusterOutputSpec(TraitedSpec):
    index_file = File(desc='output of cluster index (in size order)')
    threshold_file = File(desc='thresholded image')
    localmax_txt_file = File(desc='local maxima text file')
    localmax_vol_file = File(desc='output of local maxima volume')
    size_file = File(desc='filename for output of size image')
    max_file = File(desc='filename for output of max image')
    mean_file = File(desc='filename for output of mean image')
    pval_file = File(desc='filename for image output of log pvals')


class Cluster(FSLCommand):
    """ Uses FSL cluster to perform clustering on statistical output

    Examples
    --------

    >>> cl = Cluster()
    >>> cl.inputs.threshold = 2.3
    >>> cl.inputs.in_file = 'zstat1.nii.gz'
    >>> cl.inputs.out_localmax_txt_file = 'stats.txt'
    >>> cl.cmdline
    'cluster --in=zstat1.nii.gz --olmax=stats.txt --thresh=2.3000000000'

    """
    input_spec = ClusterInputSpec
    output_spec = ClusterOutputSpec
    _cmd = 'cluster'

    filemap = {'out_index_file': 'index', 'out_threshold_file':'threshold',
               'out_localmax_txt_file': 'localmax.txt',
               'out_localmax_vol_file': 'localmax',
               'out_size_file': 'size', 'out_max_file': 'max',
               'out_mean_file': 'mean', 'out_pval_file': 'pval'}

    def _list_outputs(self):
        outputs = self.output_spec().get()
        for key, suffix in self.filemap.items():
            outkey = key[4:]
            inval = getattr(self.inputs, key)
            if isdefined(inval):
                if isinstance(inval, bool):
                    if inval:
                        change_ext = True
                        if suffix.endswith('.txt'):
                            change_ext=False
                        outputs[outkey] = self._gen_fname(self.inputs.in_file,
                                                          suffix='_' + suffix,
                                                          change_ext=change_ext)
                else:
                    outputs[outkey] = os.pardir.abspath(inval)
        return outputs

    def _format_arg(self, name, spec, value):
        if name in self.filemap.keys():
            if isinstance(value, bool):
                fname = self._list_outputs()[name[4:]]
            else:
                fname = value
            return spec.argstr % fname
        return super(Cluster, self)._format_arg(name, spec, value)


class RandomiseInputSpec(FSLCommandInputSpec):
    in_file = File(exists=True, desc='4D input file', argstr='-i %s', position=0, mandatory=True)
    base_name = traits.Str('tbss_', desc='the rootname that all generated files will have',
                          argstr='-o "%s"', position=1, usedefault=True)
    design_mat = File(exists=True, desc='design matrix file', argstr='-d %s', position=2)
    tcon = File(exists=True, desc='t contrasts file', argstr='-t %s', position=3)
    fcon = File(exists=True, desc='f contrasts file', argstr='-f %s')
    mask = File(exists=True, desc='mask image', argstr='-m %s')
    x_block_labels = File(exists=True, desc='exchangeability block labels file', argstr='-e %s')
    demean = traits.Bool(desc='demean data temporally before model fitting', argstr='-D')
    one_sample_group_mean = traits.Bool(desc='perform 1-sample group-mean test instead of generic permutation test',
                                  argstr='-1')
    show_total_perms = traits.Bool(desc='print out how many unique permutations would be generated and exit',
                                 argstr='-q')
    show_info_parallel_mode = traits.Bool(desc='print out information required for parallel mode and exit',
                                  argstr='-Q')
    vox_p_values = traits.Bool(desc='output voxelwise (corrected and uncorrected) p-value images',
                            argstr='-x')
    tfce = traits.Bool(desc='carry out Threshold-Free Cluster Enhancement', argstr='-T')
    tfce2D = traits.Bool(desc='carry out Threshold-Free Cluster Enhancement with 2D optimisation',
                         argstr='--T2')
    f_only = traits.Bool(desc='calculate f-statistics only', argstr='--f_only')
    raw_stats_imgs = traits.Bool(desc='output raw ( unpermuted ) statistic images', argstr='-R')
    p_vec_n_dist_files = traits.Bool(desc='output permutation vector and null distribution text files',
                                 argstr='-P')
    num_perm = traits.Int(argstr='-n %d', desc='number of permutations (default 5000, set to 0 for exhaustive)')
    seed = traits.Int(argstr='--seed %d', desc='specific integer seed for random number generator')
    var_smooth = traits.Int(argstr='-v %d', desc='use variance smoothing (std is in mm)')
    c_thresh = traits.Float(argstr='-c %.2f', desc='carry out cluster-based thresholding')
    cm_thresh = traits.Float(argstr='-C %.2f', desc='carry out cluster-mass-based thresholding')
    f_c_thresh = traits.Float(argstr='-F %.2f', desc='carry out f cluster thresholding')
    f_cm_thresh = traits.Float(argstr='-S %.2f', desc='carry out f cluster-mass thresholding')
    tfce_H = traits.Float(argstr='--tfce_H %.2f', desc='TFCE height parameter (default=2)')
    tfce_E = traits.Float(argstr='--tfce_E %.2f', desc='TFCE extent parameter (default=0.5)')
    tfce_C = traits.Float(argstr='--tfce_C %.2f', desc='TFCE connectivity (6 or 26; default=6)')
    vxl = traits.List(traits.Int, argstr='--vxl %d', desc='list of numbers indicating voxelwise EVs' +
                      'position in the design matrix (list order corresponds to files in vxf option)')
    vxf = traits.List(traits.Int, argstr='--vxf %d', desc='list of 4D images containing voxelwise EVs' +
                      '(list order corresponds to numbers in vxl option)')


class RandomiseOutputSpec(TraitedSpec):
    tstat_files = traits.List(
        File(exists=True),
        desc='t contrast raw statistic')
    fstat_files = traits.List(
        File(exists=True),
        desc='f contrast raw statistic')
    t_p_files = traits.List(
        File(exists=True),
        desc='f contrast uncorrected p values files')
    f_p_files = traits.List(
        File(exists=True),
        desc='f contrast uncorrected p values files')
    t_corrected_p_files = traits.List(
        File(exists=True),
        desc='t contrast FWE (Family-wise error) corrected p values files')
    f_corrected_p_files = traits.List(
        File(exists=True),
        desc='f contrast FWE (Family-wise error) corrected p values files')


class Randomise(FSLCommand):
    """XXX UNSTABLE DO NOT USE

    FSL Randomise: feeds the 4D projected FA data into GLM
    modelling and thresholding
    in order to find voxels which correlate with your model

    Example
    -------
    >>> import nipype.interfaces.fsl as fsl
    >>> rand = fsl.Randomise(in_file='allFA.nii', \
    mask = 'mask.nii', \
    tcon='design.con', \
    design_mat='design.mat')
    >>> rand.cmdline
    'randomise -i allFA.nii -o tbss_ -d design.mat -t design.con -m mask.nii'

    """

    _cmd = 'randomise'
    input_spec = RandomiseInputSpec
    output_spec = RandomiseOutputSpec

    def _list_outputs(self):
        outputs = self.output_spec().get()
        outputs['tstat_files'] = glob(self._gen_fname(\
                '%s_tstat*.nii' % self.inputs.base_name))
        outputs['fstat_files'] = glob(self._gen_fname(\
                '%s_fstat*.nii' % self.inputs.base_name))
        prefix = False
        if self.inputs.tfce or self.inputs.tfce2D:
            prefix = 'tfce'
        elif self.inputs.vox_p_values:
            prefix = 'vox'
        elif self.inputs.c_thresh or self.inputs.f_c_thresh:
            prefix = 'clustere'
        elif self.inputs.cm_thresh or self.inputs.f_cm_thresh:
            prefix = 'clusterm'
        if prefix:
            outputs['t_p_files'] = glob(self._gen_fname(\
                     '%s_%s_p_tstat*' % (self.inputs.base_name, prefix)))
            outputs['t_corrected_p_files'] = glob(self._gen_fname(\
                '%s_%s_corrp_tstat*.nii' % (self.inputs.base_name, prefix)))

            outputs['f_p_files'] = glob(self._gen_fname(\
                '%s_%s_p_fstat*.nii' % (self.inputs.base_name, prefix)))
            outputs['f_corrected_p_files'] = glob(self._gen_fname(\
                '%s_%s_corrp_fstat*.nii' % (self.inputs.base_name, prefix)))
        return outputs
