# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
'''
Miscellaneous algorithms

    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
import os.path as op

import nibabel as nb
import numpy as np
from math import floor, ceil
from scipy.ndimage.morphology import grey_dilation
from scipy.ndimage.morphology import binary_erosion
from scipy.spatial.distance import cdist, euclidean, dice, jaccard
from scipy.ndimage.measurements import center_of_mass, label
from scipy.special import legendre
import scipy.io as sio
import itertools
import scipy.stats as stats

from nipype import logging

import warnings

import metrics as nam
from ..interfaces.base import (BaseInterface, traits, TraitedSpec, File,
                               InputMultiPath, OutputMultiPath,
                               BaseInterfaceInputSpec, isdefined,
                               DynamicTraitedSpec)
from nipype.utils.filemanip import fname_presuffix, split_filename
iflogger = logging.getLogger('interface')


class PickAtlasInputSpec(BaseInterfaceInputSpec):
    atlas = File(exists=True, desc="Location of the atlas that will be used.",
                 mandatory=True)
    labels = traits.Either(
        traits.Int, traits.List(traits.Int),
        desc=("Labels of regions that will be included in the mask. Must be\
        compatible with the atlas used."),
        mandatory=True
    )
    hemi = traits.Enum(
        'both', 'left', 'right',
        desc="Restrict the mask to only one hemisphere: left or right",
        usedefault=True
    )
    dilation_size = traits.Int(
        usedefault=True,
        desc="Defines how much the mask will be dilated (expanded in 3D)."
    )
    output_file = File(desc="Where to store the output mask.")


class PickAtlasOutputSpec(TraitedSpec):
    mask_file = File(exists=True, desc="output mask file")


class PickAtlas(BaseInterface):
    """Returns ROI masks given an atlas and a list of labels. Supports dilation
    and left right masking (assuming the atlas is properly aligned).
    """

    input_spec = PickAtlasInputSpec
    output_spec = PickAtlasOutputSpec

    def _run_interface(self, runtime):
        nim = self._get_brodmann_area()
        nb.save(nim, self._gen_output_filename())

        return runtime

    def _gen_output_filename(self):
        if not isdefined(self.inputs.output_file):
            output = fname_presuffix(fname=self.inputs.atlas, suffix="_mask",
                                     newpath=os.getcwd(), use_ext=True)
        else:
            output = os.path.realpath(self.inputs.output_file)
        return output

    def _get_brodmann_area(self):
        nii = nb.load(self.inputs.atlas)
        origdata = nii.get_data()
        newdata = np.zeros(origdata.shape)

        if not isinstance(self.inputs.labels, list):
            labels = [self.inputs.labels]
        else:
            labels = self.inputs.labels
        for lab in labels:
            newdata[origdata == lab] = 1
        if self.inputs.hemi == 'right':
            newdata[floor(float(origdata.shape[0]) / 2):, :, :] = 0
        elif self.inputs.hemi == 'left':
            newdata[:ceil(float(origdata.shape[0]) / 2), :, :] = 0

        if self.inputs.dilation_size != 0:
            newdata = grey_dilation(
                newdata, (2 * self.inputs.dilation_size + 1,
                          2 * self.inputs.dilation_size +
                          1,
                          2 * self.inputs.dilation_size + 1))

        return nb.Nifti1Image(newdata, nii.get_affine(), nii.get_header())

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['mask_file'] = self._gen_output_filename()
        return outputs


class SimpleThresholdInputSpec(BaseInterfaceInputSpec):
    volumes = InputMultiPath(
        File(exists=True), desc='volumes to be thresholded', mandatory=True)
    threshold = traits.Float(
        desc='volumes to be thresholdedeverything below this value will be set\
        to zero',
        mandatory=True
    )


class SimpleThresholdOutputSpec(TraitedSpec):
    thresholded_volumes = OutputMultiPath(
        File(exists=True), desc="thresholded volumes")


class SimpleThreshold(BaseInterface):
    """Applies a threshold to input volumes
    """
    input_spec = SimpleThresholdInputSpec
    output_spec = SimpleThresholdOutputSpec

    def _run_interface(self, runtime):
        for fname in self.inputs.volumes:
            img = nb.load(fname)
            data = np.array(img.get_data())

            active_map = data > self.inputs.threshold

            thresholded_map = np.zeros(data.shape)
            thresholded_map[active_map] = data[active_map]

            new_img = nb.Nifti1Image(
                thresholded_map, img.get_affine(), img.get_header())
            _, base, _ = split_filename(fname)
            nb.save(new_img, base + '_thresholded.nii')

        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs["thresholded_volumes"] = []
        for fname in self.inputs.volumes:
            _, base, _ = split_filename(fname)
            outputs["thresholded_volumes"].append(
                os.path.abspath(base + '_thresholded.nii'))
        return outputs


class ModifyAffineInputSpec(BaseInterfaceInputSpec):
    volumes = InputMultiPath(
        File(exists=True),
        desc='volumes which affine matrices will be modified',
        mandatory=True
    )
    transformation_matrix = traits.Array(
        value=np.eye(4),
        shape=(4, 4),
        desc="transformation matrix that will be left multiplied by the\
        affine matrix",
        usedefault=True
    )


class ModifyAffineOutputSpec(TraitedSpec):
    transformed_volumes = OutputMultiPath(File(exist=True))


class ModifyAffine(BaseInterface):
    """Left multiplies the affine matrix with a specified values. Saves the volume
    as a nifti file.
    """
    input_spec = ModifyAffineInputSpec
    output_spec = ModifyAffineOutputSpec

    def _gen_output_filename(self, name):
        _, base, _ = split_filename(name)
        return os.path.abspath(base + "_transformed.nii")

    def _run_interface(self, runtime):
        for fname in self.inputs.volumes:
            img = nb.load(fname)

            affine = img.get_affine()
            affine = np.dot(self.inputs.transformation_matrix, affine)

            nb.save(nb.Nifti1Image(img.get_data(), affine,
                    img.get_header()), self._gen_output_filename(fname))

        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['transformed_volumes'] = []
        for fname in self.inputs.volumes:
            outputs['transformed_volumes'].append(
                self._gen_output_filename(fname))
        return outputs


class CreateNiftiInputSpec(BaseInterfaceInputSpec):
    data_file = File(exists=True, mandatory=True, desc="ANALYZE img file")
    header_file = File(
        exists=True, mandatory=True, desc="corresponding ANALYZE hdr file")
    affine = traits.Array(desc="affine transformation array")


class CreateNiftiOutputSpec(TraitedSpec):
    nifti_file = File(exists=True)


class CreateNifti(BaseInterface):
    """Creates a nifti volume
    """
    input_spec = CreateNiftiInputSpec
    output_spec = CreateNiftiOutputSpec

    def _gen_output_file_name(self):
        _, base, _ = split_filename(self.inputs.data_file)
        return os.path.abspath(base + ".nii")

    def _run_interface(self, runtime):
        hdr = nb.AnalyzeHeader.from_fileobj(
            open(self.inputs.header_file, 'rb'))

        if isdefined(self.inputs.affine):
            affine = self.inputs.affine
        else:
            affine = None

        data = hdr.data_from_fileobj(open(self.inputs.data_file, 'rb'))
        img = nb.Nifti1Image(data, affine, hdr)
        nb.save(img, self._gen_output_file_name())

        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['nifti_file'] = self._gen_output_file_name()
        return outputs


class TSNRInputSpec(BaseInterfaceInputSpec):
    in_file = InputMultiPath(File(exists=True), mandatory=True,
                             desc='realigned 4D file or a list of 3D files')
    regress_poly = traits.Range(low=1, desc='Remove polynomials')


class TSNROutputSpec(TraitedSpec):
    tsnr_file = File(exists=True, desc='tsnr image file')
    mean_file = File(exists=True, desc='mean image file')
    stddev_file = File(exists=True, desc='std dev image file')
    detrended_file = File(desc='detrended input file')


class TSNR(BaseInterface):
    """Computes the time-course SNR for a time series

    Typically you want to run this on a realigned time-series.

    Example
    -------

    >>> tsnr = TSNR()
    >>> tsnr.inputs.in_file = 'functional.nii'
    >>> res = tsnr.run() # doctest: +SKIP

    """
    input_spec = TSNRInputSpec
    output_spec = TSNROutputSpec

    def _gen_output_file_name(self, suffix=None):
        _, base, ext = split_filename(self.inputs.in_file[0])
        if suffix in ['mean', 'stddev']:
            return os.path.abspath(base + "_tsnr_" + suffix + ext)
        elif suffix in ['detrended']:
            return os.path.abspath(base + "_" + suffix + ext)
        else:
            return os.path.abspath(base + "_tsnr" + ext)

    def _run_interface(self, runtime):
        img = nb.load(self.inputs.in_file[0])
        header = img.get_header().copy()
        vollist = [nb.load(filename) for filename in self.inputs.in_file]
        data = np.concatenate([vol.get_data().reshape(
            vol.get_shape()[:3] + (-1,)) for vol in vollist], axis=3)
        if data.dtype.kind == 'i':
            header.set_data_dtype(np.float32)
            data = data.astype(np.float32)
        if isdefined(self.inputs.regress_poly):
            timepoints = img.get_shape()[-1]
            X = np.ones((timepoints, 1))
            for i in range(self.inputs.regress_poly):
                X = np.hstack((X, legendre(
                    i + 1)(np.linspace(-1, 1, timepoints))[:, None]))
            betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2))
            datahat = np.rollaxis(np.dot(X[:, 1:],
                                         np.rollaxis(
                                             betas[1:, :, :, :], 0, 3)),
                                  0, 4)
            data = data - datahat
            img = nb.Nifti1Image(data, img.get_affine(), header)
            nb.save(img, self._gen_output_file_name('detrended'))
        meanimg = np.mean(data, axis=3)
        stddevimg = np.std(data, axis=3)
        tsnr = meanimg / stddevimg
        img = nb.Nifti1Image(tsnr, img.get_affine(), header)
        nb.save(img, self._gen_output_file_name())
        img = nb.Nifti1Image(meanimg, img.get_affine(), header)
        nb.save(img, self._gen_output_file_name('mean'))
        img = nb.Nifti1Image(stddevimg, img.get_affine(), header)
        nb.save(img, self._gen_output_file_name('stddev'))
        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['tsnr_file'] = self._gen_output_file_name()
        outputs['mean_file'] = self._gen_output_file_name('mean')
        outputs['stddev_file'] = self._gen_output_file_name('stddev')
        if isdefined(self.inputs.regress_poly):
            outputs['detrended_file'] = self._gen_output_file_name('detrended')
        return outputs


class GunzipInputSpec(BaseInterfaceInputSpec):
    in_file = File(exists=True, mandatory=True)


class GunzipOutputSpec(TraitedSpec):
    out_file = File(exists=True)


class Gunzip(BaseInterface):
    """Gunzip wrapper
    """
    input_spec = GunzipInputSpec
    output_spec = GunzipOutputSpec

    def _gen_output_file_name(self):
        _, base, ext = split_filename(self.inputs.in_file)
        if ext[-2:].lower() == ".gz":
            ext = ext[:-3]
        return os.path.abspath(base + ext[:-3])

    def _run_interface(self, runtime):
        import gzip
        in_file = gzip.open(self.inputs.in_file, 'rb')
        out_file = open(self._gen_output_file_name(), 'wb')
        out_file.write(in_file.read())
        out_file.close()
        in_file.close()
        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['out_file'] = self._gen_output_file_name()
        return outputs


def replaceext(in_list, ext):
    out_list = list()
    for filename in in_list:
        path, name, _ = split_filename(op.abspath(filename))
        out_name = op.join(path, name) + ext
        out_list.append(out_name)
    return out_list


def matlab2csv(in_array, name, reshape):
    output_array = np.asarray(in_array)
    if reshape:
        if len(np.shape(output_array)) > 1:
            output_array = np.reshape(output_array, (
                np.shape(output_array)[0]*np.shape(output_array)[1], 1))
            iflogger.info(np.shape(output_array))
    output_name = op.abspath(name + '.csv')
    np.savetxt(output_name, output_array, delimiter=',')
    return output_name


class Matlab2CSVInputSpec(TraitedSpec):
    in_file = File(exists=True, mandatory=True, desc='Input MATLAB .mat file')
    reshape_matrix = traits.Bool(
        True, usedefault=True,
        desc='The output of this interface is meant for R, so matrices will be\
        reshaped to vectors by default.'
    )


class Matlab2CSVOutputSpec(TraitedSpec):
    csv_files = OutputMultiPath(
        File(desc='Output CSV files for each variable saved in the input .mat\
        file')
    )


class Matlab2CSV(BaseInterface):
    """Simple interface to save the components of a MATLAB .mat file as a text
    file with comma-separated values (CSVs).

    CSV files are easily loaded in R, for use in statistical processing.
    For further information, see cran.r-project.org/doc/manuals/R-data.pdf

    Example
    -------

    >>> from nipype.algorithms import misc
    >>> mat2csv = misc.Matlab2CSV()
    >>> mat2csv.inputs.in_file = 'cmatrix.mat'
    >>> mat2csv.run() # doctest: +SKIP
    """
    input_spec = Matlab2CSVInputSpec
    output_spec = Matlab2CSVOutputSpec

    def _run_interface(self, runtime):
        in_dict = sio.loadmat(op.abspath(self.inputs.in_file))

        # Check if the file has multiple variables in it. If it does, loop
        # through them and save them as individual CSV files.
        # If not, save the variable as a single CSV file using the input file
        # name and a .csv extension.

        saved_variables = list()
        for key in in_dict.keys():
            if not key.startswith('__'):
                if isinstance(in_dict[key][0], np.ndarray):
                    saved_variables.append(key)
                else:
                    iflogger.info('One of the keys in the input file, {k}, is not a Numpy array'.format(k=key))

        if len(saved_variables) > 1:
            iflogger.info(
                '{N} variables found:'.format(N=len(saved_variables)))
            iflogger.info(saved_variables)
            for variable in saved_variables:
                iflogger.info(
                    '...Converting {var} - type {ty} - to\
                    CSV'.format(var=variable, ty=type(in_dict[variable]))
                )
                matlab2csv(
                    in_dict[variable], variable, self.inputs.reshape_matrix)
        elif len(saved_variables) == 1:
            _, name, _ = split_filename(self.inputs.in_file)
            variable = saved_variables[0]
            iflogger.info('Single variable found {var}, type {ty}:'.format(
                var=variable, ty=type(in_dict[variable])))
            iflogger.info('...Converting {var} to CSV from {f}'.format(
                var=variable, f=self.inputs.in_file))
            matlab2csv(in_dict[variable], name, self.inputs.reshape_matrix)
        else:
            iflogger.error('No values in the MATLAB file?!')
        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        in_dict = sio.loadmat(op.abspath(self.inputs.in_file))
        saved_variables = list()
        for key in in_dict.keys():
            if not key.startswith('__'):
                if isinstance(in_dict[key][0], np.ndarray):
                    saved_variables.append(key)
                else:
                    iflogger.error('One of the keys in the input file, {k}, is\
                                   not a Numpy array'.format(k=key))

        if len(saved_variables) > 1:
            outputs['csv_files'] = replaceext(saved_variables, '.csv')
        elif len(saved_variables) == 1:
            _, name, ext = split_filename(self.inputs.in_file)
            outputs['csv_files'] = op.abspath(name + '.csv')
        else:
            iflogger.error('No values in the MATLAB file?!')
        return outputs


def merge_csvs(in_list):
    for idx, in_file in enumerate(in_list):
        try:
            in_array = np.loadtxt(in_file, delimiter=',')
        except ValueError, ex:
            try:
                in_array = np.loadtxt(in_file, delimiter=',', skiprows=1)
            except ValueError, ex:
                first = open(in_file, 'r')
                header_line = first.readline()
                header_list = header_line.split(',')
                n_cols = len(header_list)
                try:
                    in_array = np.loadtxt(
                        in_file, delimiter=',', skiprows=1,
                        usecols=range(1, n_cols)
                    )
                except ValueError, ex:
                    in_array = np.loadtxt(
                        in_file, delimiter=',', skiprows=1, usecols=range(1, n_cols-1))
        if idx == 0:
            out_array = in_array
        else:
            out_array = np.dstack((out_array, in_array))
    out_array = np.squeeze(out_array)
    iflogger.info('Final output array shape:')
    iflogger.info(np.shape(out_array))
    return out_array


def remove_identical_paths(in_files):
    import os.path as op
    from nipype.utils.filemanip import split_filename
    if len(in_files) > 1:
        out_names = list()
        commonprefix = op.commonprefix(in_files)
        lastslash = commonprefix.rfind('/')
        commonpath = commonprefix[0:(lastslash+1)]
        for fileidx, in_file in enumerate(in_files):
            path, name, ext = split_filename(in_file)
            in_file = op.join(path, name)
            name = in_file.replace(commonpath, '')
            name = name.replace('_subject_id_', '')
            out_names.append(name)
    else:
        path, name, ext = split_filename(in_files[0])
        out_names = [name]
    return out_names


def maketypelist(rowheadings, shape, extraheadingBool, extraheading):
    typelist = []
    if rowheadings:
        typelist.append(('heading', 'a40'))
    if len(shape) > 1:
        for idx in range(1, (min(shape)+1)):
            typelist.append((str(idx), float))
    else:
        for idx in range(1, (shape[0]+1)):
            typelist.append((str(idx), float))
    if extraheadingBool:
        typelist.append((extraheading, 'a40'))
    iflogger.info(typelist)
    return typelist


def makefmtlist(output_array, typelist, rowheadingsBool,
                shape, extraheadingBool):
    fmtlist = []
    if rowheadingsBool:
        fmtlist.append('%s')
    if len(shape) > 1:
        output = np.zeros(max(shape), typelist)
        for idx in range(1, min(shape)+1):
            output[str(idx)] = output_array[:, idx-1]
            fmtlist.append('%f')
    else:
        output = np.zeros(1, typelist)
        for idx in range(1, len(output_array)+1):
            output[str(idx)] = output_array[idx-1]
            fmtlist.append('%f')
    if extraheadingBool:
        fmtlist.append('%s')
    fmt = ','.join(fmtlist)
    return fmt, output


class MergeCSVFilesInputSpec(TraitedSpec):
    in_files = InputMultiPath(File(exists=True), mandatory=True,
                              desc='Input comma-separated value (CSV) files')
    out_file = File('merged.csv', usedefault=True,
                    desc='Output filename for merged CSV file')
    column_headings = traits.List(
        traits.Str, desc='List of column headings to save in merged CSV file\
        (must be equal to number of input files). If left undefined, these\
        will be pulled from the input filenames.')
    row_headings = traits.List(
        traits.Str, desc='List of row headings to save in merged CSV file\
        (must be equal to number of rows in the input files).')
    row_heading_title = traits.Str(
        'label', usedefault=True, desc='Column heading for the row headings\
         added')
    extra_column_heading = traits.Str(
        desc='New heading to add for the added field.')
    extra_field = traits.Str(
        desc='New field to add to each row. This is useful for saving the\
        group or subject ID in the file.')


class MergeCSVFilesOutputSpec(TraitedSpec):
    csv_file = File(desc='Output CSV file containing columns ')


class MergeCSVFiles(BaseInterface):
    """This interface is designed to facilitate data loading in the R environment.
    It takes input CSV files and merges them into a single CSV file.
    If provided, it will also incorporate column heading names into the
    resulting CSV file.

    CSV files are easily loaded in R, for use in statistical processing.
    For further information, see cran.r-project.org/doc/manuals/R-data.pdf

    Example
    -------

    >>> from nipype.algorithms import misc
    >>> mat2csv = misc.MergeCSVFiles()
    >>> mat2csv.inputs.in_files = ['degree.mat','clustering.mat']
    >>> mat2csv.inputs.column_headings = ['degree','clustering']
    >>> mat2csv.run() # doctest: +SKIP
    """
    input_spec = MergeCSVFilesInputSpec
    output_spec = MergeCSVFilesOutputSpec

    def _run_interface(self, runtime):
        extraheadingBool = False
        extraheading = ''
        rowheadingsBool = False
        """
        This block defines the column headings.
        """
        if isdefined(self.inputs.column_headings):
            iflogger.info('Column headings have been provided:')
            headings = self.inputs.column_headings
        else:
            iflogger.info(
                'Column headings not provided! Pulled from input filenames:')
            headings = remove_identical_paths(self.inputs.in_files)

        if isdefined(self.inputs.extra_field):
            if isdefined(self.inputs.extra_column_heading):
                extraheading = self.inputs.extra_column_heading
                iflogger.info('Extra column heading provided: {col}'.format(
                    col=extraheading))
            else:
                extraheading = 'type'
                iflogger.info(
                    'Extra column heading was not defined. Using "type"')
            headings.append(extraheading)
            extraheadingBool = True

        if len(self.inputs.in_files) == 1:
            iflogger.warn('Only one file input!')

        if isdefined(self.inputs.row_headings):
            iflogger.info('Row headings have been provided. Adding "labels"\
                          column header.')
            prefix = '"{p}","'.format(p=self.inputs.row_heading_title)
            csv_headings = prefix + '","'.join(itertools.chain(
                headings)) + '"\n'
            rowheadingsBool = True
        else:
            iflogger.info('Row headings have not been provided.')
            csv_headings = '"' + '","'.join(itertools.chain(headings)) + '"\n'

        iflogger.info('Final Headings:')
        iflogger.info(csv_headings)

        """
        Next we merge the arrays and define the output text file
        """

        output_array = merge_csvs(self.inputs.in_files)
        _, name, ext = split_filename(self.inputs.out_file)
        if not ext == '.csv':
            ext = '.csv'

        out_file = op.abspath(name + ext)
        file_handle = open(out_file, 'w')
        file_handle.write(csv_headings)

        shape = np.shape(output_array)
        typelist = maketypelist(
            rowheadingsBool, shape, extraheadingBool, extraheading)
        fmt, output = makefmtlist(
            output_array, typelist, rowheadingsBool, shape, extraheadingBool)

        if rowheadingsBool:
            row_heading_list = self.inputs.row_headings
            row_heading_list_with_quotes = []
            for row_heading in row_heading_list:
                row_heading_with_quotes = '"' + row_heading + '"'
                row_heading_list_with_quotes.append(row_heading_with_quotes)
            row_headings = np.array(row_heading_list_with_quotes, dtype='|S40')
            output['heading'] = row_headings

        if isdefined(self.inputs.extra_field):
            extrafieldlist = []
            if len(shape) > 1:
                mx = shape[0]
            else:
                mx = 1
            for idx in range(0, mx):
                extrafieldlist.append(self.inputs.extra_field)
            iflogger.info(len(extrafieldlist))
            output[extraheading] = extrafieldlist
        iflogger.info(output)
        iflogger.info(fmt)
        np.savetxt(file_handle, output, fmt, delimiter=',')
        file_handle.close()
        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        _, name, ext = split_filename(self.inputs.out_file)
        if not ext == '.csv':
            ext = '.csv'
        out_file = op.abspath(name + ext)
        outputs['csv_file'] = out_file
        return outputs


class AddCSVColumnInputSpec(TraitedSpec):
    in_file = File(exists=True, mandatory=True,
                   desc='Input comma-separated value (CSV) files')
    out_file = File('extra_heading.csv', usedefault=True,
                    desc='Output filename for merged CSV file')
    extra_column_heading = traits.Str(
        desc='New heading to add for the added field.')
    extra_field = traits.Str(
        desc='New field to add to each row. This is useful for saving the\
        group or subject ID in the file.')


class AddCSVColumnOutputSpec(TraitedSpec):
    csv_file = File(desc='Output CSV file containing columns ')


class AddCSVColumn(BaseInterface):
    """Short interface to add an extra column and field to a text file

    Example
    -------

    >>> from nipype.algorithms import misc
    >>> addcol = misc.AddCSVColumn()
    >>> addcol.inputs.in_file = 'degree.csv'
    >>> addcol.inputs.extra_column_heading = 'group'
    >>> addcol.inputs.extra_field = 'male'
    >>> addcol.run() # doctest: +SKIP
    """
    input_spec = AddCSVColumnInputSpec
    output_spec = AddCSVColumnOutputSpec

    def _run_interface(self, runtime):
        in_file = open(self.inputs.in_file, 'r')
        _, name, ext = split_filename(self.inputs.out_file)
        if not ext == '.csv':
            ext = '.csv'
        out_file = op.abspath(name + ext)

        out_file = open(out_file, 'w')
        firstline = in_file.readline()
        firstline = firstline.replace('\n', '')
        new_firstline = firstline + ',"' + \
            self.inputs.extra_column_heading + '"\n'
        out_file.write(new_firstline)
        for line in in_file:
            new_line = line.replace('\n', '')
            new_line = new_line + ',' + self.inputs.extra_field + '\n'
            out_file.write(new_line)
        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        _, name, ext = split_filename(self.inputs.out_file)
        if not ext == '.csv':
            ext = '.csv'
        out_file = op.abspath(name + ext)
        outputs['csv_file'] = out_file
        return outputs


class AddCSVRowInputSpec(DynamicTraitedSpec, BaseInterfaceInputSpec):
    in_file = traits.File(mandatory=True, desc='Input comma-separated value (CSV) files')
    _outputs = traits.Dict(traits.Any, value={}, usedefault=True)

    def __setattr__(self, key, value):
        if key not in self.copyable_trait_names():
            if not isdefined(value):
                super(AddCSVRowInputSpec, self).__setattr__(key, value)
            self._outputs[key] = value
        else:
            if key in self._outputs:
                self._outputs[key] = value
            super(AddCSVRowInputSpec, self).__setattr__(key, value)


class AddCSVRowOutputSpec(TraitedSpec):
    csv_file = File(desc='Output CSV file containing rows ')


class AddCSVRow(BaseInterface):
    """Simple interface to add an extra row to a csv file

    .. note:: Requires `pandas <http://pandas.pydata.org/>`_

    .. warning:: Multi-platform thread-safe execution is possible with
        `lockfile <https://pythonhosted.org/lockfile/lockfile.html>`_. Please recall that (1)
        this module is alpha software; and (2) it should be installed for thread-safe writing.
        If lockfile is not installed, then the interface is not thread-safe.


    Example
    -------

    >>> from nipype.algorithms import misc
    >>> addrow = misc.AddCSVRow()
    >>> addrow.inputs.in_file = 'scores.csv'
    >>> addrow.inputs.si = 0.74
    >>> addrow.inputs.di = 0.93
    >>> addrow.subject_id = 'S400'
    >>> addrow.inputs.list_of_values = [ 0.4, 0.7, 0.3 ]
    >>> addrow.run() # doctest: +SKIP
    """
    input_spec = AddCSVRowInputSpec
    output_spec = AddCSVRowOutputSpec

    def __init__(self, infields=None, force_run=True, **kwargs):
        super(AddCSVRow, self).__init__(**kwargs)
        undefined_traits = {}
        self._infields = infields
        self._have_lock = False
        self._lock = None

        if infields:
            for key in infields:
                self.inputs.add_trait(key, traits.Any)
                self.inputs._outputs[key] = Undefined
                undefined_traits[key] = Undefined
        self.inputs.trait_set(trait_change_notify=False, **undefined_traits)

        if force_run:
            self._always_run = True

    def _run_interface(self, runtime):
        try:
            import pandas as pd
        except ImportError:
            raise ImportError('This interface requires pandas (http://pandas.pydata.org/) to run.')

        try:
            import lockfile as pl
            self._have_lock = True
        except ImportError:
            import warnings
            warnings.warn(('Python module lockfile was not found: AddCSVRow will not be thread-safe '
                          'in multi-processor execution'))

        input_dict = {}
        for key, val in self.inputs._outputs.items():
            # expand lists to several columns
            if isinstance(val, list):
                for i,v in enumerate(val):
                    input_dict['%s_%d' % (key,i)]=v
            else:
                input_dict[key] = val

        df = pd.DataFrame([input_dict])

        if self._have_lock:
            self._lock = pl.FileLock(self.inputs.in_file)

            # Acquire lock
            self._lock.acquire()

        if op.exists(self.inputs.in_file):
            formerdf = pd.read_csv(self.inputs.in_file, index_col=0)
            df = pd.concat([formerdf, df], ignore_index=True)

        with open(self.inputs.in_file, 'w') as f:
            df.to_csv(f)

        if self._have_lock:
            self._lock.release()

        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        outputs['csv_file'] = self.inputs.in_file
        return outputs

    def _outputs(self):
        return self._add_output_traits(super(AddCSVRow, self)._outputs())

    def _add_output_traits(self, base):
        return base


class CalculateNormalizedMomentsInputSpec(TraitedSpec):
    timeseries_file = File(
        exists=True, mandatory=True,
        desc='Text file with timeseries in columns and timepoints in rows,\
        whitespace separated')
    moment = traits.Int(
        mandatory=True,
        desc="Define which moment should be calculated, 3 for skewness, 4 for\
        kurtosis.")


class CalculateNormalizedMomentsOutputSpec(TraitedSpec):
    moments = traits.List(traits.Float(), desc='Moments')


class CalculateNormalizedMoments(BaseInterface):
    """Calculates moments of timeseries.

    Example
    -------

    >>> from nipype.algorithms import misc
    >>> skew = misc.CalculateNormalizedMoments()
    >>> skew.inputs.moment = 3
    >>> skew.inputs.timeseries_file = 'timeseries.txt'
    >>> skew.run() # doctest: +SKIP
    """
    input_spec = CalculateNormalizedMomentsInputSpec
    output_spec = CalculateNormalizedMomentsOutputSpec

    def _run_interface(self, runtime):

        self._moments = calc_moments(
            self.inputs.timeseries_file, self.inputs.moment)
        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        outputs['skewness'] = self._moments
        return outputs


def calc_moments(timeseries_file, moment):
    """Returns nth moment (3 for skewness, 4 for kurtosis) of timeseries
    (list of values; one per timeseries).

    Keyword arguments:
    timeseries_file -- text file with white space separated timepoints in rows

    """
    timeseries = np.genfromtxt(timeseries_file)

    m2 = stats.moment(timeseries, 2, axis=0)
    m3 = stats.moment(timeseries, moment, axis=0)
    zero = (m2 == 0)
    return np.where(zero, 0, m3 / m2**(moment/2.0))


class AddNoiseInputSpec(TraitedSpec):
    in_file = File(exists=True, mandatory=True,
                   desc='input image that will be corrupted with noise')
    in_mask = File(exists=True, desc=('input mask, voxels outside this mask '
                   'will be considered background'))
    snr = traits.Float(10.0, desc='desired output SNR in dB', usedefault=True)
    dist = traits.Enum('normal', 'rician', usedefault=True, mandatory=True,
                       desc=('desired noise distribution'))
    bg_dist = traits.Enum('normal', 'rayleigh', usedefault=True, mandatory=True,
                          desc=('desired noise distribution, currently '
                          'only normal is implemented'))
    out_file = File(desc='desired output filename')


class AddNoiseOutputSpec(TraitedSpec):
    out_file = File(exists=True, desc='corrupted image')


class AddNoise(BaseInterface):
    """
    Corrupts with noise the input image


    Example
    -------
    >>> from nipype.algorithms.misc import AddNoise
    >>> noise = AddNoise()
    >>> noise.inputs.in_file = 'T1.nii'
    >>> noise.inputs.in_mask = 'mask.nii'
    >>> noise.snr = 30.0
    >>> noise.run() # doctest: +SKIP

    """
    input_spec = AddNoiseInputSpec
    output_spec = AddNoiseOutputSpec

    def _run_interface(self, runtime):
        in_image = nb.load(self.inputs.in_file)
        in_data = in_image.get_data()
        snr = self.inputs.snr

        if isdefined(self.inputs.in_mask):
            in_mask = nb.load(self.inputs.in_mask).get_data()
        else:
            in_mask = np.ones_like(in_data)

        result = self.gen_noise(in_data, mask=in_mask, snr_db=snr,
                                dist=self.inputs.dist, bg_dist=self.inputs.bg_dist)
        res_im = nb.Nifti1Image(result, in_image.get_affine(), in_image.get_header())
        res_im.to_filename(self._gen_output_filename())
        return runtime

    def _gen_output_filename(self):
        if not isdefined(self.inputs.out_file):
            _, base, ext = split_filename(self.inputs.in_file)
            out_file = os.path.abspath('%s_SNR%03.2f%s' % (base, self.inputs.snr, ext))
        else:
            out_file = self.inputs.out_file

        return out_file

    def _list_outputs(self):
        outputs = self.output_spec().get()
        outputs['out_file'] = self._gen_output_filename()
        return outputs

    def gen_noise(self, image, mask=None, snr_db=10.0, dist='normal', bg_dist='normal'):
        """
        Generates a copy of an image with a certain amount of
        added gaussian noise (rayleigh for background in mask)
        """
        from math import sqrt
        snr = sqrt(np.power(10.0, snr_db/10.0))

        if mask is None:
            mask = np.ones_like(image)
        else:
            mask[mask > 0] = 1
            mask[mask < 1] = 0

            if mask.ndim < image.ndim:
                mask = np.rollaxis(np.array([mask]*image.shape[3]), 0, 4)

        signal = image[mask > 0].reshape(-1)

        if dist == 'normal':
            signal = signal - signal.mean()
            sigma_n = sqrt(signal.var()/snr)
            noise = np.random.normal(size=image.shape, scale=sigma_n)

            if (np.any(mask == 0)) and (bg_dist == 'rayleigh'):
                bg_noise = np.random.rayleigh(size=image.shape, scale=sigma_n)
                noise[mask == 0] = bg_noise[mask == 0]

            im_noise = image + noise

        elif dist == 'rician':
            sigma_n = signal.mean()/snr
            n_1 = np.random.normal(size=image.shape, scale=sigma_n)
            n_2 = np.random.normal(size=image.shape, scale=sigma_n)
            stde_1 = n_1/sqrt(2.0)
            stde_2 = n_2/sqrt(2.0)
            im_noise = np.sqrt((image + stde_1)**2 + (stde_2)**2)
        else:
            raise NotImplementedError(('Only normal and rician distributions '
                                      'are supported'))

        return im_noise


class NormalizeProbabilityMapSetInputSpec(TraitedSpec):
    in_files = InputMultiPath(File(exists=True, mandatory=True,
                              desc='The tpms to be normalized'))
    in_mask = File(exists=True,
                   desc='Masked voxels must sum up 1.0, 0.0 otherwise.')


class NormalizeProbabilityMapSetOutputSpec(TraitedSpec):
    out_files = OutputMultiPath(File(exists=True),
                                desc="normalized maps")


class NormalizeProbabilityMapSet(BaseInterface):
    """ Returns the input tissue probability maps (tpms, aka volume fractions)
    normalized to sum up 1.0 at each voxel within the mask.

    .. note:: Please recall this is not a spatial normalization algorithm


    Example
    -------

    >>> from nipype.algorithms import misc
    >>> normalize = misc.NormalizeProbabilityMapSet()
    >>> normalize.inputs.in_files = [ 'tpm_00.nii.gz', 'tpm_01.nii.gz', \
'tpm_02.nii.gz' ]
    >>> normalize.inputs.in_mask = 'tpms_msk.nii.gz'
    >>> normalize.run() # doctest: +SKIP
    """
    input_spec = NormalizeProbabilityMapSetInputSpec
    output_spec = NormalizeProbabilityMapSetOutputSpec

    def _run_interface(self, runtime):
        mask = None

        if isdefined(self.inputs.in_mask):
            mask = self.inputs.in_mask

        self._out_filenames = normalize_tpms(self.inputs.in_files, mask)
        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        outputs['out_files'] = self._out_filenames
        return outputs


class SplitROIsInputSpec(TraitedSpec):
    in_file = File(exists=True, mandatory=True,
                   desc='file to be splitted')
    in_mask = File(exists=True, desc='only process files inside mask')
    roi_size = traits.Tuple(traits.Int, traits.Int, traits.Int,
                            desc='desired ROI size')


class SplitROIsOutputSpec(TraitedSpec):
    out_files = OutputMultiPath(File(exists=True),
                                desc='the resulting ROIs')
    out_masks = OutputMultiPath(File(exists=True),
                                desc='a mask indicating valid values')
    out_index = OutputMultiPath(File(exists=True),
                                desc='arrays keeping original locations')


class SplitROIs(BaseInterface):
    """
    Splits a 3D image in small chunks to enable parallel processing.
    ROIs keep time series structure in 4D images.
    >>> from nipype.algorithms import misc
    >>> rois = misc.SplitROIs()
    >>> rois.inputs.in_file = 'diffusion.nii'
    >>> rois.inputs.in_mask = 'mask.nii'
    >>> rois.run() # doctest: +SKIP

    """
    input_spec = SplitROIsInputSpec
    output_spec = SplitROIsOutputSpec

    def _run_interface(self, runtime):
        mask = None
        roisize = None
        self._outnames = {}

        if isdefined(self.inputs.in_mask):
            mask = self.inputs.in_mask
        if isdefined(self.inputs.roi_size):
            roisize = self.inputs.roi_size

        res = split_rois(self.inputs.in_file,
                         mask, roisize)
        self._outnames['out_files'] = res[0]
        self._outnames['out_masks'] = res[1]
        self._outnames['out_index'] = res[2]
        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        for k, v in self._outnames.iteritems():
            outputs[k] = v
        return outputs


class MergeROIsInputSpec(TraitedSpec):
    in_files = InputMultiPath(File(exists=True, mandatory=True,
                                   desc='files to be re-merged'))
    in_index = InputMultiPath(File(exists=True, mandatory=True),
                              desc='array keeping original locations')
    in_reference = File(exists=True, desc='reference file')


class MergeROIsOutputSpec(TraitedSpec):
    merged_file = File(exists=True, desc='the recomposed file')


class MergeROIs(BaseInterface):
    """
    Splits a 3D image in small chunks to enable parallel processing.
    ROIs keep time series structure in 4D images.

    Example
    -------

    >>> from nipype.algorithms import misc
    >>> rois = misc.MergeROIs()
    >>> rois.inputs.in_files = ['roi%02d.nii' % i for i in xrange(1, 6)]
    >>> rois.inputs.in_reference = 'mask.nii'
    >>> rois.inputs.in_index = ['roi%02d_idx.npz' % i for i in xrange(1, 6)]
    >>> rois.run() # doctest: +SKIP

    """
    input_spec = MergeROIsInputSpec
    output_spec = MergeROIsOutputSpec

    def _run_interface(self, runtime):
        res = merge_rois(self.inputs.in_files,
                         self.inputs.in_index,
                         self.inputs.in_reference)
        self._merged = res
        return runtime

    def _list_outputs(self):
        outputs = self.output_spec().get()
        outputs['merged_file'] = self._merged
        return outputs


def normalize_tpms(in_files, in_mask=None, out_files=[]):
    """
    Returns the input tissue probability maps (tpms, aka volume fractions)
    normalized to sum up 1.0 at each voxel within the mask.
    """
    import nibabel as nib
    import numpy as np
    import os.path as op

    in_files = np.atleast_1d(in_files).tolist()

    if len(out_files) != len(in_files):
        for i,finname in enumerate(in_files):
            fname,fext = op.splitext(op.basename(finname))
            if fext == '.gz':
                fname,fext2 = op.splitext(fname)
                fext = fext2 + fext

            out_file = op.abspath('%s_norm_%02d%s' % (fname,i,fext))
            out_files+= [out_file]

    imgs = [nib.load(fim) for fim in in_files]

    if len(in_files)==1:
        img_data = imgs[0].get_data()
        img_data[img_data>0.0] = 1.0
        hdr = imgs[0].get_header().copy()
        hdr['data_type']= 16
        hdr.set_data_dtype(np.float32)
        nib.save(nib.Nifti1Image(img_data.astype(np.float32), imgs[0].get_affine(), hdr), out_files[0])
        return out_files[0]

    img_data = np.array([im.get_data() for im in imgs]).astype(np.float32)
    #img_data[img_data>1.0] = 1.0
    img_data[img_data<0.0] = 0.0
    weights = np.sum(img_data, axis=0)

    msk = np.ones_like(imgs[0].get_data())
    msk[ weights<= 0 ] = 0

    if not in_mask is None:
        msk = nib.load(in_mask).get_data()
        msk[ msk<=0 ] = 0
        msk[ msk>0 ] = 1

    msk = np.ma.masked_equal(msk, 0)


    for i,out_file in enumerate(out_files):
        data = np.ma.masked_equal(img_data[i], 0)
        probmap = data / weights
        hdr = imgs[i].get_header().copy()
        hdr['data_type']= 16
        hdr.set_data_dtype('float32')
        nib.save(nib.Nifti1Image(probmap.astype(np.float32), imgs[i].get_affine(), hdr), out_file)

    return out_files


def split_rois(in_file, mask=None, roishape=None):
    """
    Splits an image in ROIs for parallel processing
    """
    import nibabel as nb
    import numpy as np
    from math import sqrt, ceil
    import os.path as op

    if roishape is None:
        roishape = (10, 10, 1)

    im = nb.load(in_file)
    imshape = im.get_shape()
    dshape = imshape[:3]
    nvols = imshape[-1]
    roisize = roishape[0] * roishape[1] * roishape[2]
    droishape = (roishape[0], roishape[1], roishape[2], nvols)

    if mask is not None:
        mask = nb.load(mask).get_data()
        mask[mask > 0] = 1
        mask[mask < 1] = 0
    else:
        mask = np.ones(dshape)

    mask = mask.reshape(-1).astype(np.uint8)
    nzels = np.nonzero(mask)
    els = np.sum(mask)
    nrois = int(ceil(els/roisize))

    data = im.get_data().reshape((mask.size, -1))
    data = np.squeeze(data.take(nzels, axis=0))
    nvols = data.shape[-1]

    roidefname = op.abspath('onesmask.nii.gz')
    nb.Nifti1Image(np.ones(roishape, dtype=np.uint8), None,
                   None).to_filename(roidefname)

    out_files = []
    out_mask = []
    out_idxs = []

    for i in xrange(nrois):
        first = i * roisize
        last = (i+1) * roisize
        fill = 0

        if last > els:
            fill = last - els
            last = els

        droi = data[first:last, ...]
        iname = op.abspath('roi%010d_idx' % i)
        out_idxs.append(iname+'.npz')
        np.savez(iname, (nzels[0][first:last],))

        if fill > 0:
            droi = np.vstack((droi, np.zeros((fill, nvols), dtype=np.float32)))
            partialmsk = np.ones((roisize,), dtype=np.uint8)
            partialmsk[-fill:] = 0
            partname = op.abspath('partialmask.nii.gz')
            nb.Nifti1Image(partialmsk.reshape(roishape), None,
                           None).to_filename(partname)
            out_mask.append(partname)
        else:
            out_mask.append(roidefname)

        fname = op.abspath('roi%010d.nii.gz' % i)
        nb.Nifti1Image(droi.reshape(droishape),
                       None, None).to_filename(fname)
        out_files.append(fname)
    return out_files, out_mask, out_idxs


def merge_rois(in_files, in_idxs, in_ref,
               dtype=None, out_file=None):
    """
    Re-builds an image resulting from a parallelized processing
    """
    import nibabel as nb
    import numpy as np
    import os.path as op
    import subprocess as sp

    if out_file is None:
        out_file = op.abspath('merged.nii.gz')

    if dtype is None:
        dtype = np.float32

    # if file is compressed, uncompress using os
    # to avoid memory errors
    if op.splitext(in_ref)[1] == '.gz':
        try:
            iflogger.info('uncompress %i' % in_ref)
            sp.check_call(['gunzip', in_ref], stdout=sp.PIPE, shell=True)
            in_ref = op.splitext(in_ref)[0]
        except:
            pass

    ref = nb.load(in_ref)
    aff = ref.get_affine()
    hdr = ref.get_header().copy()
    rsh = ref.get_shape()
    del ref
    npix = rsh[0] * rsh[1] * rsh[2]
    fcdata = nb.load(in_files[0]).get_data()

    if fcdata.ndim == 4:
        ndirs = fcdata.shape[-1]
    else:
        ndirs = 1
    newshape = (rsh[0], rsh[1], rsh[2], ndirs)
    hdr.set_data_dtype(dtype)
    hdr.set_xyzt_units('mm', 'sec')

    if ndirs < 300:
        data = np.zeros((npix, ndirs))
        for cname, iname in zip(in_files, in_idxs):
            with np.load(iname) as f:
                idxs = np.squeeze(f['arr_0'])
            cdata = nb.load(cname).get_data().reshape(-1, ndirs)
            nels = len(idxs)
            idata = (idxs, )
            try:
                data[idata, ...] = cdata[0:nels, ...]
            except:
                print(data.shape, cdata.shape)
                raise

        hdr.set_data_shape(newshape)

        nb.Nifti1Image(data.reshape(newshape).astype(dtype),
                       aff, hdr).to_filename(out_file)


    else:
        hdr.set_data_shape(rsh[:3])
        nii = []
        for d in xrange(ndirs):
            fname = op.abspath('vol%06d.nii' % d)
            nb.Nifti1Image(np.zeros(rsh[:3]), aff, hdr).to_filename(fname)
            nii.append(fname)

        for cname, iname in zip(in_files, in_idxs):
            with np.load(iname) as f:
                idxs = np.squeeze(f['arr_0'])

            for d, fname in enumerate(nii):
                data = nb.load(fname).get_data().reshape(-1)
                cdata = nb.load(cname).get_data().reshape(-1, ndirs)[:, d]
                nels = len(idxs)
                idata = (idxs, )
                data[idata] = cdata[0:nels]
                nb.Nifti1Image(data.reshape(rsh[:3]), aff, hdr).to_filename(fname)

        imgs = [nb.load(im) for im in nii]
        allim = nb.concat_images(imgs)
        allim.to_filename(out_file)

    return out_file


# Deprecated interfaces ------------------------------------------------------

class Distance(nam.Distance):
    """Calculates distance between two volumes.

    .. deprecated:: 0.10.0
       Use :py:class:`nipype.algorithms.metrics.Distance` instead.
    """
    def __init__(self, **inputs):
        super(nam.Distance, self).__init__(**inputs)
        warnings.warn(("This interface has been deprecated since 0.10.0,"
                      " please use nipype.algorithms.metrics.Distance"),
                      DeprecationWarning)


class Overlap(nam.Overlap):
    """Calculates various overlap measures between two maps.

    .. deprecated:: 0.10.0
       Use :py:class:`nipype.algorithms.metrics.Overlap` instead.
    """
    def __init__(self, **inputs):
        super(nam.Overlap, self).__init__(**inputs)
        warnings.warn(("This interface has been deprecated since 0.10.0,"
                      " please use nipype.algorithms.metrics.Overlap"),
                      DeprecationWarning)


class FuzzyOverlap(nam.FuzzyOverlap):
    """Calculates various overlap measures between two maps, using a fuzzy
    definition.

    .. deprecated:: 0.10.0
       Use :py:class:`nipype.algorithms.metrics.FuzzyOverlap` instead.
    """
    def __init__(self, **inputs):
        super(nam.FuzzyOverlap, self).__init__(**inputs)
        warnings.warn(("This interface has been deprecated since 0.10.0,"
                      " please use nipype.algorithms.metrics.FuzzyOverlap"),
                      DeprecationWarning)
