<reponame>mwaskom/nipype<filename>nipype/algorithms/confounds.py<gh_stars>0
# -*- coding: utf-8 -*-
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
'''
Algorithms to compute confounds in :abbr:`fMRI (functional MRI)`

    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)

'''
from __future__ import print_function, division, unicode_literals, absolute_import
from builtins import range

import os
import os.path as op

import nibabel as nb
import numpy as np
from numpy.polynomial import Legendre
from scipy import linalg

from .. import logging
from ..external.due import BibTeX
from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
                               BaseInterfaceInputSpec, File, isdefined,
                               InputMultiPath, OutputMultiPath)
from ..utils import NUMPY_MMAP
from ..utils.misc import normalize_mc_params

IFLOG = logging.getLogger('interface')


class ComputeDVARSInputSpec(BaseInterfaceInputSpec):
    in_file = File(exists=True, mandatory=True, desc='functional data, after HMC')
    in_mask = File(exists=True, mandatory=True, desc='a brain mask')
    remove_zerovariance = traits.Bool(True, usedefault=True,
                                      desc='remove voxels with zero variance')
    save_std = traits.Bool(True, usedefault=True,
                           desc='save standardized DVARS')
    save_nstd = traits.Bool(False, usedefault=True,
                            desc='save non-standardized DVARS')
    save_vxstd = traits.Bool(False, usedefault=True,
                             desc='save voxel-wise standardized DVARS')
    save_all = traits.Bool(False, usedefault=True, desc='output all DVARS')

    series_tr = traits.Float(desc='repetition time in sec.')
    save_plot = traits.Bool(False, usedefault=True, desc='write DVARS plot')
    figdpi = traits.Int(100, usedefault=True, desc='output dpi for the plot')
    figsize = traits.Tuple(traits.Float(11.7), traits.Float(2.3), usedefault=True,
                           desc='output figure size')
    figformat = traits.Enum('png', 'pdf', 'svg', usedefault=True,
                            desc='output format for figures')
    intensity_normalization = traits.Float(1000.0, usedefault=True,
                              desc='Divide value in each voxel at each timepoint '
                                   'by the median calculated across all voxels'
                                   'and timepoints within the mask (if specified)'
                                   'and then multiply by the value specified by'
                                   'this parameter. By using the default (1000)' \
                                   'output DVARS will be expressed in ' \
                                   'x10 % BOLD units compatible with Power et al.' \
                                   '2012. Set this to 0 to disable intensity' \
                                   'normalization altogether.')



class ComputeDVARSOutputSpec(TraitedSpec):
    out_std = File(exists=True, desc='output text file')
    out_nstd = File(exists=True, desc='output text file')
    out_vxstd = File(exists=True, desc='output text file')
    out_all = File(exists=True, desc='output text file')
    avg_std = traits.Float()
    avg_nstd = traits.Float()
    avg_vxstd = traits.Float()
    fig_std = File(exists=True, desc='output DVARS plot')
    fig_nstd = File(exists=True, desc='output DVARS plot')
    fig_vxstd = File(exists=True, desc='output DVARS plot')


class ComputeDVARS(BaseInterface):
    """
    Computes the DVARS.
    """
    input_spec = ComputeDVARSInputSpec
    output_spec = ComputeDVARSOutputSpec
    references_ = [{
        'entry': BibTeX("""\
@techreport{nichols_notes_2013,
    address = {Coventry, UK},
    title = {Notes on {Creating} a {Standardized} {Version} of {DVARS}},
    url = {http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-\
research/nichols/scripts/fsl/standardizeddvars.pdf},
    urldate = {2016-08-16},
    institution = {University of Warwick},
    author = {<NAME>},
    year = {2013}
}"""),
        'tags': ['method']
    }, {
        'entry': BibTeX("""\
@article{power_spurious_2012,
    title = {Spurious but systematic correlations in functional connectivity {MRI} networks \
arise from subject motion},
    volume = {59},
    doi = {10.1016/j.neuroimage.2011.10.018},
    number = {3},
    urldate = {2016-08-16},
    journal = {NeuroImage},
    author = {Power, <NAME>. and Barnes, <NAME>. and Snyder, <NAME>. and Schlaggar, \
<NAME>. and Petersen, <NAME>.},
    year = {2012},
    pages = {2142--2154},
}
"""),
        'tags': ['method']
    }]

    def __init__(self, **inputs):
        self._results = {}
        super(ComputeDVARS, self).__init__(**inputs)

    def _gen_fname(self, suffix, ext=None):
        fname, in_ext = op.splitext(op.basename(
            self.inputs.in_file))

        if in_ext == '.gz':
            fname, in_ext2 = op.splitext(fname)
            in_ext = in_ext2 + in_ext

        if ext is None:
            ext = in_ext

        if ext.startswith('.'):
            ext = ext[1:]

        return op.abspath('{}_{}.{}'.format(fname, suffix, ext))

    def _run_interface(self, runtime):
        dvars = compute_dvars(self.inputs.in_file, self.inputs.in_mask,
                              remove_zerovariance=self.inputs.remove_zerovariance,
                              intensity_normalization=self.inputs.intensity_normalization)

        (self._results['avg_std'],
         self._results['avg_nstd'],
         self._results['avg_vxstd']) = np.mean(dvars, axis=1).astype(float)

        tr = None
        if isdefined(self.inputs.series_tr):
            tr = self.inputs.series_tr

        if self.inputs.save_std:
            out_file = self._gen_fname('dvars_std', ext='tsv')
            np.savetxt(out_file, dvars[0], fmt=b'%0.6f')
            self._results['out_std'] = out_file

            if self.inputs.save_plot:
                self._results['fig_std'] = self._gen_fname(
                    'dvars_std', ext=self.inputs.figformat)
                fig = plot_confound(dvars[0], self.inputs.figsize, 'Standardized DVARS',
                                    series_tr=tr)
                fig.savefig(self._results['fig_std'], dpi=float(self.inputs.figdpi),
                        format=self.inputs.figformat,
                        bbox_inches='tight')
                fig.clf()

        if self.inputs.save_nstd:
            out_file = self._gen_fname('dvars_nstd', ext='tsv')
            np.savetxt(out_file, dvars[1], fmt=b'%0.6f')
            self._results['out_nstd'] = out_file

            if self.inputs.save_plot:
                self._results['fig_nstd'] = self._gen_fname(
                    'dvars_nstd', ext=self.inputs.figformat)
                fig = plot_confound(dvars[1], self.inputs.figsize, 'DVARS', series_tr=tr)
                fig.savefig(self._results['fig_nstd'], dpi=float(self.inputs.figdpi),
                            format=self.inputs.figformat,
                            bbox_inches='tight')
                fig.clf()

        if self.inputs.save_vxstd:
            out_file = self._gen_fname('dvars_vxstd', ext='tsv')
            np.savetxt(out_file, dvars[2], fmt=b'%0.6f')
            self._results['out_vxstd'] = out_file

            if self.inputs.save_plot:
                self._results['fig_vxstd'] = self._gen_fname(
                    'dvars_vxstd', ext=self.inputs.figformat)
                fig = plot_confound(dvars[2], self.inputs.figsize, 'Voxelwise std DVARS',
                                    series_tr=tr)
                fig.savefig(self._results['fig_vxstd'], dpi=float(self.inputs.figdpi),
                            format=self.inputs.figformat,
                            bbox_inches='tight')
                fig.clf()

        if self.inputs.save_all:
            out_file = self._gen_fname('dvars', ext='tsv')
            np.savetxt(out_file, np.vstack(dvars).T, fmt=b'%0.8f', delimiter=b'\t',
                       header='std DVARS\tnon-std DVARS\tvx-wise std DVARS', comments='')
            self._results['out_all'] = out_file

        return runtime

    def _list_outputs(self):
        return self._results


class FramewiseDisplacementInputSpec(BaseInterfaceInputSpec):
    in_file = File(exists=True, mandatory=True, desc='motion parameters')
    parameter_source = traits.Enum("FSL", "AFNI", "SPM", "FSFAST", "NIPY",
                                   desc="Source of movement parameters",
                                   mandatory=True)
    radius = traits.Float(50, usedefault=True,
                          desc='radius in mm to calculate angular FDs, 50mm is the '
                               'default since it is used in Power et al. 2012')
    out_file = File('fd_power_2012.txt', usedefault=True, desc='output file name')
    out_figure = File('fd_power_2012.pdf', usedefault=True, desc='output figure name')
    series_tr = traits.Float(desc='repetition time in sec.')
    save_plot = traits.Bool(False, usedefault=True, desc='write FD plot')
    normalize = traits.Bool(False, usedefault=True, desc='calculate FD in mm/s')
    figdpi = traits.Int(100, usedefault=True, desc='output dpi for the FD plot')
    figsize = traits.Tuple(traits.Float(11.7), traits.Float(2.3), usedefault=True,
                           desc='output figure size')


class FramewiseDisplacementOutputSpec(TraitedSpec):
    out_file = File(desc='calculated FD per timestep')
    out_figure = File(desc='output image file')
    fd_average = traits.Float(desc='average FD')


class FramewiseDisplacement(BaseInterface):
    """
    Calculate the :abbr:`FD (framewise displacement)` as in [Power2012]_.
    This implementation reproduces the calculation in fsl_motion_outliers

    .. [Power2012] Power et al., Spurious but systematic correlations in functional
         connectivity MRI networks arise from subject motion, NeuroImage 59(3),
         2012. doi:`10.1016/j.neuroimage.2011.10.018
         <http://dx.doi.org/10.1016/j.neuroimage.2011.10.018>`_.


    """

    input_spec = FramewiseDisplacementInputSpec
    output_spec = FramewiseDisplacementOutputSpec

    references_ = [{
        'entry': BibTeX("""\
@article{power_spurious_2012,
    title = {Spurious but systematic correlations in functional connectivity {MRI} networks \
arise from subject motion},
    volume = {59},
    doi = {10.1016/j.neuroimage.2011.10.018},
    number = {3},
    urldate = {2016-08-16},
    journal = {NeuroImage},
    author = {Power, <NAME>. and Barnes, <NAME>. and Snyder, <NAME>. and Schlaggar, \
<NAME>. and Petersen, <NAME>.},
    year = {2012},
    pages = {2142--2154},
}
"""),
        'tags': ['method']
    }]

    def _run_interface(self, runtime):
        mpars = np.loadtxt(self.inputs.in_file)  # mpars is N_t x 6
        mpars = np.apply_along_axis(func1d=normalize_mc_params,
                                    axis=1, arr=mpars,
                                    source=self.inputs.parameter_source)
        diff = mpars[:-1, :6] - mpars[1:, :6]
        diff[:, 3:6] *= self.inputs.radius
        fd_res = np.abs(diff).sum(axis=1)

        self._results = {
            'out_file': op.abspath(self.inputs.out_file),
            'fd_average': float(fd_res.mean())
        }
        np.savetxt(self.inputs.out_file, fd_res, header='FramewiseDisplacement', comments='')

        if self.inputs.save_plot:
            tr = None
            if isdefined(self.inputs.series_tr):
                tr = self.inputs.series_tr

            if self.inputs.normalize and tr is None:
                IFLOG.warn('FD plot cannot be normalized if TR is not set')

            self._results['out_figure'] = op.abspath(self.inputs.out_figure)
            fig = plot_confound(fd_res, self.inputs.figsize, 'FD', units='mm',
                                series_tr=tr, normalize=self.inputs.normalize)
            fig.savefig(self._results['out_figure'], dpi=float(self.inputs.figdpi),
                        format=self.inputs.out_figure[-3:],
                        bbox_inches='tight')
            fig.clf()

        return runtime

    def _list_outputs(self):
        return self._results


class CompCorInputSpec(BaseInterfaceInputSpec):
    realigned_file = File(exists=True, mandatory=True,
                          desc='already realigned brain image (4D)')
    mask_files = InputMultiPath(File(exists=True),
                                desc=('One or more mask files that determines '
                                      'ROI (3D). When more that one file is '
                                      'provided `merge_method` or '
                                      '`merge_index` must be provided'))
    merge_method = traits.Enum('union', 'intersect', 'none', xor=['mask_index'],
                               requires=['mask_files'],
                               desc=('Merge method if multiple masks are '
                                     'present - `union` uses voxels included in'
                                     ' at least one input mask, `intersect` '
                                     'uses only voxels present in all input '
                                     'masks, `none` performs CompCor on '
                                     'each mask individually'))
    mask_index = traits.Range(low=0, xor=['merge_method'],
                              requires=['mask_files'],
                              desc=('Position of mask in `mask_files` to use - '
                                    'first is the default.'))
    components_file = traits.Str('components_file.txt', usedefault=True,
                           desc='Filename to store physiological components')
    num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL
    use_regress_poly = traits.Bool(True, usedefault=True,
                                   desc=('use polynomial regression '
                                         'pre-component extraction'))
    regress_poly_degree = traits.Range(low=1, default=1, usedefault=True,
                                       desc='the degree polynomial to use')
    header_prefix = traits.Str(desc=('the desired header for the output tsv '
                                     'file (one column). If undefined, will '
                                     'default to "CompCor"'))


class CompCorOutputSpec(TraitedSpec):
    components_file = File(exists=True,
                           desc='text file containing the noise components')


class CompCor(BaseInterface):
    """
    Interface with core CompCor computation, used in aCompCor and tCompCor

    Example
    -------

    >>> ccinterface = CompCor()
    >>> ccinterface.inputs.realigned_file = 'functional.nii'
    >>> ccinterface.inputs.mask_files = 'mask.nii'
    >>> ccinterface.inputs.num_components = 1
    >>> ccinterface.inputs.use_regress_poly = True
    >>> ccinterface.inputs.regress_poly_degree = 2

    """
    input_spec = CompCorInputSpec
    output_spec = CompCorOutputSpec
    references_ = [{'entry': BibTeX("@article{compcor_2007,"
                                    "title = {A component based noise correction method (CompCor) for BOLD and perfusion based},"
                                    "volume = {37},"
                                    "number = {1},"
                                    "doi = {10.1016/j.neuroimage.2007.04.042},"
                                    "urldate = {2016-08-13},"
                                    "journal = {NeuroImage},"
                                    "author = {<NAME> and <NAME> <NAME>.},"
                                    "year = {2007},"
                                    "pages = {90-101},}"
                                   ),
                    'tags': ['method', 'implementation']
                   }]

    def __init__(self, *args, **kwargs):
        ''' exactly the same as compcor except the header '''
        super(CompCor, self).__init__(*args, **kwargs)
        self._header = 'CompCor'

    def _run_interface(self, runtime):
        mask_images = []
        if isdefined(self.inputs.mask_files):
            mask_images = combine_mask_files(self.inputs.mask_files,
                                             self.inputs.merge_method,
                                             self.inputs.mask_index)

        degree = (self.inputs.regress_poly_degree if
                  self.inputs.use_regress_poly else 0)

        imgseries = nb.load(self.inputs.realigned_file,
                            mmap=NUMPY_MMAP)

        if len(imgseries.shape) != 4:
            raise ValueError('tCompCor expected a 4-D nifti file. Input {} has '
                            '{} dimensions (shape {})'.format(
                            self.inputs.realigned_file, len(imgseries.shape),
                            imgseries.shape))

        if len(mask_images) == 0:
            img = nb.Nifti1Image(np.ones(imgseries.shape[:3], dtype=np.bool),
                                 affine=imgseries.affine,
                                 header=imgseries.header)
            mask_images = [img]

        mask_images = self._process_masks(mask_images, imgseries.get_data())

        components = compute_noise_components(imgseries.get_data(),
                                              mask_images, degree,
                                              self.inputs.num_components)

        components_file = os.path.join(os.getcwd(), self.inputs.components_file)
        np.savetxt(components_file, components, fmt=b"%.10f", delimiter='\t',
                   header=self._make_headers(components.shape[1]), comments='')
        return runtime

    def _process_masks(self, mask_images, timeseries=None):
        return mask_images

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['components_file'] = os.path.abspath(self.inputs.components_file)
        return outputs

    def _make_headers(self, num_col):
        headers = []
        header = self.inputs.header_prefix if \
            isdefined(self.inputs.header_prefix) else self._header
        for i in range(num_col):
            headers.append(header + '{:02d}'.format(i))
        return '\t'.join(headers)


class ACompCor(CompCor):
    """
    Anatomical compcor: for inputs and outputs, see CompCor.
    When the mask provided is an anatomical mask, then CompCor
    is equivalent to ACompCor.
    """

    def __init__(self, *args, **kwargs):
        ''' exactly the same as compcor except the header '''
        super(ACompCor, self).__init__(*args, **kwargs)
        self._header = 'aCompCor'


class TCompCorInputSpec(CompCorInputSpec):
    # and all the fields in CompCorInputSpec
    percentile_threshold = traits.Range(low=0., high=1., value=.02,
                                        exclude_low=True, exclude_high=True,
                                        usedefault=True, desc='the percentile '
                                        'used to select highest-variance '
                                        'voxels, represented by a number '
                                        'between 0 and 1, exclusive. By '
                                        'default, this value is set to .02. '
                                        'That is, the 2% of voxels '
                                        'with the highest variance are used.')


class TCompCorOutputSpec(CompCorOutputSpec):
    # and all the fields in CompCorOutputSpec
    high_variance_masks = OutputMultiPath(File(exists=True),
                                          desc=(("voxels exceeding the variance"
                                                 " threshold")))


class TCompCor(CompCor):
    """
    Interface for tCompCor. Computes a ROI mask based on variance of voxels.

    Example
    -------

    >>> ccinterface = TCompCor()
    >>> ccinterface.inputs.realigned_file = 'functional.nii'
    >>> ccinterface.inputs.mask_files = 'mask.nii'
    >>> ccinterface.inputs.num_components = 1
    >>> ccinterface.inputs.use_regress_poly = True
    >>> ccinterface.inputs.regress_poly_degree = 2
    >>> ccinterface.inputs.percentile_threshold = .03

    """

    input_spec = TCompCorInputSpec
    output_spec = TCompCorOutputSpec

    def __init__(self, *args, **kwargs):
        ''' exactly the same as compcor except the header '''
        super(TCompCor, self).__init__(*args, **kwargs)
        self._header = 'tCompCor'
        self._mask_files = []

    def _process_masks(self, mask_images, timeseries=None):
        out_images = []
        self._mask_files = []
        for i, img in enumerate(mask_images):
            mask = img.get_data().astype(np.bool)
            imgseries = timeseries[mask, :]
            imgseries = regress_poly(2, imgseries)
            tSTD = _compute_tSTD(imgseries, 0, axis=-1)
            threshold_std = np.percentile(tSTD, np.round(100. *
                           (1. - self.inputs.percentile_threshold)).astype(int))
            mask_data = np.zeros_like(mask)
            mask_data[mask != 0] = tSTD >= threshold_std
            out_image = nb.Nifti1Image(mask_data, affine=img.affine,
                                       header=img.header)

            # save mask
            mask_file = os.path.abspath('mask_{:03d}.nii.gz'.format(i))
            out_image.to_filename(mask_file)
            IFLOG.debug('tCompcor computed and saved mask of shape {} to '
                        'mask_file {}'.format(mask.shape, mask_file))
            self._mask_files.append(mask_file)
            out_images.append(out_image)
        return out_images

    def _list_outputs(self):
        outputs = super(TCompCor, self)._list_outputs()
        outputs['high_variance_masks'] = self._mask_files
        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')
    tsnr_file = File('tsnr.nii.gz', usedefault=True, hash_files=False,
                     desc='output tSNR file')
    mean_file = File('mean.nii.gz', usedefault=True, hash_files=False,
                     desc='output mean file')
    stddev_file = File('stdev.nii.gz', usedefault=True, hash_files=False,
                       desc='output tSNR file')
    detrended_file = File('detrend.nii.gz', usedefault=True, hash_files=False,
                          desc='input file after detrending')


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 _run_interface(self, runtime):
        img = nb.load(self.inputs.in_file[0], mmap=NUMPY_MMAP)
        header = img.header.copy()
        vollist = [nb.load(filename, mmap=NUMPY_MMAP) for filename in self.inputs.in_file]
        data = np.concatenate([vol.get_data().reshape(
            vol.shape[:3] + (-1,)) for vol in vollist], axis=3)
        data = np.nan_to_num(data)

        if data.dtype.kind == 'i':
            header.set_data_dtype(np.float32)
            data = data.astype(np.float32)

        if isdefined(self.inputs.regress_poly):
            data = regress_poly(self.inputs.regress_poly, data, remove_mean=False)
            img = nb.Nifti1Image(data, img.affine, header)
            nb.save(img, op.abspath(self.inputs.detrended_file))

        meanimg = np.mean(data, axis=3)
        stddevimg = np.std(data, axis=3)
        tsnr = np.zeros_like(meanimg)
        tsnr[stddevimg > 1.e-3] = meanimg[stddevimg > 1.e-3] / stddevimg[stddevimg > 1.e-3]
        img = nb.Nifti1Image(tsnr, img.affine, header)
        nb.save(img, op.abspath(self.inputs.tsnr_file))
        img = nb.Nifti1Image(meanimg, img.affine, header)
        nb.save(img, op.abspath(self.inputs.mean_file))
        img = nb.Nifti1Image(stddevimg, img.affine, header)
        nb.save(img, op.abspath(self.inputs.stddev_file))
        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        for k in ['tsnr_file', 'mean_file', 'stddev_file']:
            outputs[k] = op.abspath(getattr(self.inputs, k))

        if isdefined(self.inputs.regress_poly):
            outputs['detrended_file'] = op.abspath(self.inputs.detrended_file)
        return outputs


class NonSteadyStateDetectorInputSpec(BaseInterfaceInputSpec):
    in_file = File(exists=True, mandatory=True, desc='4D NIFTI EPI file')


class NonSteadyStateDetectorOutputSpec(TraitedSpec):
    n_volumes_to_discard = traits.Int(desc='Number of non-steady state volumes'
                                           'detected in the beginning of the scan.')


class NonSteadyStateDetector(BaseInterface):
    """
    Returns the number of non-steady state volumes detected at the beginning
    of the scan.
    """

    input_spec = NonSteadyStateDetectorInputSpec
    output_spec = NonSteadyStateDetectorOutputSpec

    def _run_interface(self, runtime):
        in_nii = nb.load(self.inputs.in_file)
        global_signal = in_nii.get_data()[:,:,:,:50].mean(axis=0).mean(axis=0).mean(axis=0)

        self._results = {
            'n_volumes_to_discard': _is_outlier(global_signal)
        }

        return runtime

    def _list_outputs(self):
        return self._results


def compute_dvars(in_file, in_mask, remove_zerovariance=False,
                  intensity_normalization=1000):
    """
    Compute the :abbr:`DVARS (D referring to temporal
    derivative of timecourses, VARS referring to RMS variance over voxels)`
    [Power2012]_.

    Particularly, the *standardized* :abbr:`DVARS (D referring to temporal
    derivative of timecourses, VARS referring to RMS variance over voxels)`
    [Nichols2013]_ are computed.

    .. [Nichols2013] <NAME>, `Notes on creating a standardized version of
         DVARS <http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-\
research/nichols/scripts/fsl/standardizeddvars.pdf>`_, 2013.

    .. note:: Implementation details

      Uses the implementation of the `Yule-Walker equations
      from nitime
      <http://nipy.org/nitime/api/generated/nitime.algorithms.autoregressive.html\
#nitime.algorithms.autoregressive.AR_est_YW>`_
      for the :abbr:`AR (auto-regressive)` filtering of the fMRI signal.

    :param numpy.ndarray func: functional data, after head-motion-correction.
    :param numpy.ndarray mask: a 3D mask of the brain
    :param bool output_all: write out all dvars
    :param str out_file: a path to which the standardized dvars should be saved.
    :return: the standardized DVARS

    """
    import numpy as np
    import nibabel as nb
    from nitime.algorithms import AR_est_YW
    import warnings

    func = nb.load(in_file, mmap=NUMPY_MMAP).get_data().astype(np.float32)
    mask = nb.load(in_mask, mmap=NUMPY_MMAP).get_data().astype(np.uint8)

    if len(func.shape) != 4:
        raise RuntimeError(
            "Input fMRI dataset should be 4-dimensional")

    idx = np.where(mask > 0)
    mfunc = func[idx[0], idx[1], idx[2], :]

    if intensity_normalization != 0:
        mfunc = (mfunc / np.median(mfunc)) * intensity_normalization

    # Robust standard deviation (we are using "lower" interpolation
    # because this is what FSL is doing
    func_sd = (np.percentile(mfunc, 75, axis=1, interpolation="lower") -
               np.percentile(mfunc, 25, axis=1, interpolation="lower")) / 1.349

    if remove_zerovariance:
        mfunc = mfunc[func_sd != 0, :]
        func_sd = func_sd[func_sd != 0]

    # Compute (non-robust) estimate of lag-1 autocorrelation
    ar1 = np.apply_along_axis(AR_est_YW, 1,
                              regress_poly(0, mfunc, remove_mean=True).astype(
                                  np.float32), 1)[:, 0]

    # Compute (predicted) standard deviation of temporal difference time series
    diff_sdhat = np.squeeze(np.sqrt(((1 - ar1) * 2).tolist())) * func_sd
    diff_sd_mean = diff_sdhat.mean()

    # Compute temporal difference time series
    func_diff = np.diff(mfunc, axis=1)

    # DVARS (no standardization)
    dvars_nstd = np.sqrt(np.square(func_diff).mean(axis=0))

    # standardization
    dvars_stdz = dvars_nstd / diff_sd_mean

    with warnings.catch_warnings():  # catch, e.g., divide by zero errors
        warnings.filterwarnings('error')

        # voxelwise standardization
        diff_vx_stdz = np.square(
            func_diff / np.array([diff_sdhat] * func_diff.shape[-1]).T)
        dvars_vx_stdz = np.sqrt(diff_vx_stdz.mean(axis=0))

    return (dvars_stdz, dvars_nstd, dvars_vx_stdz)


def plot_confound(tseries, figsize, name, units=None,
                  series_tr=None, normalize=False):
    """
    A helper function to plot :abbr:`fMRI (functional MRI)` confounds.

    """
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.gridspec import GridSpec
    from matplotlib.backends.backend_pdf import FigureCanvasPdf as FigureCanvas
    import seaborn as sns

    fig = plt.Figure(figsize=figsize)
    FigureCanvas(fig)
    grid = GridSpec(1, 2, width_ratios=[3, 1], wspace=0.025)
    grid.update(hspace=1.0, right=0.95, left=0.1, bottom=0.2)

    ax = fig.add_subplot(grid[0, :-1])
    if normalize and series_tr is not None:
        tseries /= series_tr

    ax.plot(tseries)
    ax.set_xlim((0, len(tseries)))
    ylabel = name
    if units is not None:
        ylabel += (' speed [{}/s]' if normalize else ' [{}]').format(units)
    ax.set_ylabel(ylabel)

    xlabel = 'Frame #'
    if series_tr is not None:
        xlabel = 'Frame # ({} sec TR)'.format(series_tr)
    ax.set_xlabel(xlabel)
    ylim = ax.get_ylim()

    ax = fig.add_subplot(grid[0, -1])
    sns.distplot(tseries, vertical=True, ax=ax)
    ax.set_xlabel('Frames')
    ax.set_ylim(ylim)
    ax.set_yticklabels([])
    return fig


def is_outlier(points, thresh=3.5):
    """
    Returns a boolean array with True if points are outliers and False
    otherwise.

    :param nparray points: an numobservations by numdimensions numpy array of observations
    :param float thresh: the modified z-score to use as a threshold. Observations with
        a modified z-score (based on the median absolute deviation) greater
        than this value will be classified as outliers.

    :return: A bolean mask, of size numobservations-length array.

    .. note:: References

        <NAME> and <NAME> (1993), "Volume 16: How to Detect and
        Handle Outliers", The ASQC Basic References in Quality Control:
        Statistical Techniques, <NAME>, Ph.D., Editor.

    """
    if len(points.shape) == 1:
        points = points[:, None]
    median = np.median(points, axis=0)
    diff = np.sum((points - median) ** 2, axis=-1)
    diff = np.sqrt(diff)
    med_abs_deviation = np.median(diff)

    modified_z_score = 0.6745 * diff / med_abs_deviation

    timepoints_to_discard = 0
    for i in range(len(modified_z_score)):
        if modified_z_score[i] <= thresh:
            break
        else:
            timepoints_to_discard += 1

    return timepoints_to_discard


def regress_poly(degree, data, remove_mean=True, axis=-1):
    """
    Returns data with degree polynomial regressed out.

    :param bool remove_mean: whether or not demean data (i.e. degree 0),
    :param int axis: numpy array axes along which regression is performed

    """
    IFLOG.debug('Performing polynomial regression on data of shape ' + str(data.shape))

    datashape = data.shape
    timepoints = datashape[axis]

    # Rearrange all voxel-wise time-series in rows
    data = data.reshape((-1, timepoints))

    # Generate design matrix
    X = np.ones((timepoints, 1)) # quick way to calc degree 0
    for i in range(degree):
        polynomial_func = Legendre.basis(i + 1)
        value_array = np.linspace(-1, 1, timepoints)
        X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis]))

    # Calculate coefficients
    betas = np.linalg.pinv(X).dot(data.T)

    # Estimation
    if remove_mean:
        datahat = X.dot(betas).T
    else: # disregard the first layer of X, which is degree 0
        datahat = X[:, 1:].dot(betas[1:, ...]).T
    regressed_data = data - datahat

    # Back to original shape
    return regressed_data.reshape(datashape)


def combine_mask_files(mask_files, mask_method=None, mask_index=None):
    """Combines input mask files into a single nibabel image

    A helper function for CompCor

    mask_files: a list
        one or more binary mask files
    mask_method: enum ('union', 'intersect', 'none')
        determines how to combine masks
    mask_index: an integer
        determines which file to return (mutually exclusive with mask_method)

    returns: a list of nibabel images
    """

    if isdefined(mask_index) or not isdefined(mask_method):
        if not isdefined(mask_index):
            if len(mask_files) == 1:
                mask_index = 0
            else:
                raise ValueError(('When more than one mask file is provided, '
                                  'one of merge_method or mask_index must be '
                                  'set'))
        if mask_index < len(mask_files):
            mask = nb.load(mask_files[mask_index], mmap=NUMPY_MMAP)
            return [mask]
        raise ValueError(('mask_index {0} must be less than number of mask '
                          'files {1}').format(mask_index, len(mask_files)))
    masks = []
    if mask_method == 'none':
        for filename in mask_files:
            masks.append(nb.load(filename, mmap=NUMPY_MMAP))
        return masks

    if mask_method == 'union':
        mask = None
        for filename in mask_files:
            img = nb.load(filename, mmap=NUMPY_MMAP)
            if mask is None:
                mask = img.get_data() > 0
            np.logical_or(mask, img.get_data() > 0, mask)
        img = nb.Nifti1Image(mask, img.affine, header=img.header)
        return [img]

    if mask_method == 'intersect':
        mask = None
        for filename in mask_files:
            img = nb.load(filename, mmap=NUMPY_MMAP)
            if mask is None:
                mask = img.get_data() > 0
            np.logical_and(mask, img.get_data() > 0, mask)
        img = nb.Nifti1Image(mask, img.affine, header=img.header)
        return [img]


def compute_noise_components(imgseries, mask_images, degree, num_components):
    """Compute the noise components from the imgseries for each mask

    imgseries: a nibabel img
    mask_images: a list of nibabel images
    degree: order of polynomial used to remove trends from the timeseries
    num_components: number of noise components to return

    returns:

    components: a numpy array

    """
    components = None
    for img in mask_images:
        mask = img.get_data().astype(np.bool)
        if imgseries.shape[:3] != mask.shape:
            raise ValueError('Inputs for CompCor, timeseries and mask, '
                             'do not have matching spatial dimensions '
                             '({} and {}, respectively)'.format(
                imgseries.shape[:3], mask.shape))

        voxel_timecourses = imgseries[mask, :]

        # Zero-out any bad values
        voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0

        # from paper:
        # "The constant and linear trends of the columns in the matrix M were
        # removed [prior to ...]"
        voxel_timecourses = regress_poly(degree, voxel_timecourses)

        # "Voxel time series from the noise ROI (either anatomical or tSTD) were
        # placed in a matrix M of size Nxm, with time along the row dimension
        # and voxels along the column dimension."
        M = voxel_timecourses.T

        # "[... were removed] prior to column-wise variance normalization."
        M = M / _compute_tSTD(M, 1.)

        # "The covariance matrix C = MMT was constructed and decomposed into its
        # principal components using a singular value decomposition."
        u, _, _ = linalg.svd(M, full_matrices=False)
        if components is None:
            components = u[:, :num_components]
        else:
            components = np.hstack((components,
                                    u[:, :num_components]))
    if components is None and num_components > 0:
        raise ValueError('No components found')
    return components


def _compute_tSTD(M, x, axis=0):
    stdM = np.std(M, axis=axis)
    # set bad values to x
    stdM[stdM == 0] = x
    stdM[np.isnan(stdM)] = x
    return stdM

# -*- coding: utf-8 -*-
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
'''
Algorithms to compute confounds in :abbr:`fMRI (functional MRI)`

    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)

'''
from __future__ import print_function, division, unicode_literals, absolute_import
from builtins import range

import os
import os.path as op

import nibabel as nb
import numpy as np
from numpy.polynomial import Legendre
from scipy import linalg

from .. import logging
from ..external.due import BibTeX
from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
                               BaseInterfaceInputSpec, File, isdefined,
                               InputMultiPath, OutputMultiPath)
from ..utils import NUMPY_MMAP
from ..utils.misc import normalize_mc_params

IFLOG = logging.getLogger('interface')


class ComputeDVARSInputSpec(BaseInterfaceInputSpec):
    in_file = File(exists=True, mandatory=True, desc='functional data, after HMC')
    in_mask = File(exists=True, mandatory=True, desc='a brain mask')
    remove_zerovariance = traits.Bool(True, usedefault=True,
                                      desc='remove voxels with zero variance')
    save_std = traits.Bool(True, usedefault=True,
                           desc='save standardized DVARS')
    save_nstd = traits.Bool(False, usedefault=True,
                            desc='save non-standardized DVARS')
    save_vxstd = traits.Bool(False, usedefault=True,
                             desc='save voxel-wise standardized DVARS')
    save_all = traits.Bool(False, usedefault=True, desc='output all DVARS')

    series_tr = traits.Float(desc='repetition time in sec.')
    save_plot = traits.Bool(False, usedefault=True, desc='write DVARS plot')
    figdpi = traits.Int(100, usedefault=True, desc='output dpi for the plot')
    figsize = traits.Tuple(traits.Float(11.7), traits.Float(2.3), usedefault=True,
                           desc='output figure size')
    figformat = traits.Enum('png', 'pdf', 'svg', usedefault=True,
                            desc='output format for figures')
    intensity_normalization = traits.Float(1000.0, usedefault=True,
                              desc='Divide value in each voxel at each timepoint '
                                   'by the median calculated across all voxels'
                                   'and timepoints within the mask (if specified)'
                                   'and then multiply by the value specified by'
                                   'this parameter. By using the default (1000)' \
                                   'output DVARS will be expressed in ' \
                                   'x10 % BOLD units compatible with Power et al.' \
                                   '2012. Set this to 0 to disable intensity' \
                                   'normalization altogether.')



class ComputeDVARSOutputSpec(TraitedSpec):
    out_std = File(exists=True, desc='output text file')
    out_nstd = File(exists=True, desc='output text file')
    out_vxstd = File(exists=True, desc='output text file')
    out_all = File(exists=True, desc='output text file')
    avg_std = traits.Float()
    avg_nstd = traits.Float()
    avg_vxstd = traits.Float()
    fig_std = File(exists=True, desc='output DVARS plot')
    fig_nstd = File(exists=True, desc='output DVARS plot')
    fig_vxstd = File(exists=True, desc='output DVARS plot')


class ComputeDVARS(BaseInterface):
    """
    Computes the DVARS.
    """
    input_spec = ComputeDVARSInputSpec
    output_spec = ComputeDVARSOutputSpec
    references_ = [{
        'entry': BibTeX("""\
@techreport{nichols_notes_2013,
    address = {Coventry, UK},
    title = {Notes on {Creating} a {Standardized} {Version} of {DVARS}},
    url = {http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-\
research/nichols/scripts/fsl/standardizeddvars.pdf},
    urldate = {2016-08-16},
    institution = {University of Warwick},
    author = {<NAME>},
    year = {2013}
}"""),
        'tags': ['method']
    }, {
        'entry': BibTeX("""\
@article{power_spurious_2012,
    title = {Spurious but systematic correlations in functional connectivity {MRI} networks \
arise from subject motion},
    volume = {59},
    doi = {10.1016/j.neuroimage.2011.10.018},
    number = {3},
    urldate = {2016-08-16},
    journal = {NeuroImage},
    author = {Power, <NAME>. and Barnes, <NAME>. and Snyder, <NAME>. and Schlaggar, \
<NAME>. and Petersen, <NAME>.},
    year = {2012},
    pages = {2142--2154},
}
"""),
        'tags': ['method']
    }]

    def __init__(self, **inputs):
        self._results = {}
        super(ComputeDVARS, self).__init__(**inputs)

    def _gen_fname(self, suffix, ext=None):
        fname, in_ext = op.splitext(op.basename(
            self.inputs.in_file))

        if in_ext == '.gz':
            fname, in_ext2 = op.splitext(fname)
            in_ext = in_ext2 + in_ext

        if ext is None:
            ext = in_ext

        if ext.startswith('.'):
            ext = ext[1:]

        return op.abspath('{}_{}.{}'.format(fname, suffix, ext))

    def _run_interface(self, runtime):
        dvars = compute_dvars(self.inputs.in_file, self.inputs.in_mask,
                              remove_zerovariance=self.inputs.remove_zerovariance,
                              intensity_normalization=self.inputs.intensity_normalization)

        (self._results['avg_std'],
         self._results['avg_nstd'],
         self._results['avg_vxstd']) = np.mean(dvars, axis=1).astype(float)

        tr = None
        if isdefined(self.inputs.series_tr):
            tr = self.inputs.series_tr

        if self.inputs.save_std:
            out_file = self._gen_fname('dvars_std', ext='tsv')
            np.savetxt(out_file, dvars[0], fmt=b'%0.6f')
            self._results['out_std'] = out_file

            if self.inputs.save_plot:
                self._results['fig_std'] = self._gen_fname(
                    'dvars_std', ext=self.inputs.figformat)
                fig = plot_confound(dvars[0], self.inputs.figsize, 'Standardized DVARS',
                                    series_tr=tr)
                fig.savefig(self._results['fig_std'], dpi=float(self.inputs.figdpi),
                        format=self.inputs.figformat,
                        bbox_inches='tight')
                fig.clf()

        if self.inputs.save_nstd:
            out_file = self._gen_fname('dvars_nstd', ext='tsv')
            np.savetxt(out_file, dvars[1], fmt=b'%0.6f')
            self._results['out_nstd'] = out_file

            if self.inputs.save_plot:
                self._results['fig_nstd'] = self._gen_fname(
                    'dvars_nstd', ext=self.inputs.figformat)
                fig = plot_confound(dvars[1], self.inputs.figsize, 'DVARS', series_tr=tr)
                fig.savefig(self._results['fig_nstd'], dpi=float(self.inputs.figdpi),
                            format=self.inputs.figformat,
                            bbox_inches='tight')
                fig.clf()

        if self.inputs.save_vxstd:
            out_file = self._gen_fname('dvars_vxstd', ext='tsv')
            np.savetxt(out_file, dvars[2], fmt=b'%0.6f')
            self._results['out_vxstd'] = out_file

            if self.inputs.save_plot:
                self._results['fig_vxstd'] = self._gen_fname(
                    'dvars_vxstd', ext=self.inputs.figformat)
                fig = plot_confound(dvars[2], self.inputs.figsize, 'Voxelwise std DVARS',
                                    series_tr=tr)
                fig.savefig(self._results['fig_vxstd'], dpi=float(self.inputs.figdpi),
                            format=self.inputs.figformat,
                            bbox_inches='tight')
                fig.clf()

        if self.inputs.save_all:
            out_file = self._gen_fname('dvars', ext='tsv')
            np.savetxt(out_file, np.vstack(dvars).T, fmt=b'%0.8f', delimiter=b'\t',
                       header='std DVARS\tnon-std DVARS\tvx-wise std DVARS', comments='')
            self._results['out_all'] = out_file

        return runtime

    def _list_outputs(self):
        return self._results


class FramewiseDisplacementInputSpec(BaseInterfaceInputSpec):
    in_file = File(exists=True, mandatory=True, desc='motion parameters')
    parameter_source = traits.Enum("FSL", "AFNI", "SPM", "FSFAST", "NIPY",
                                   desc="Source of movement parameters",
                                   mandatory=True)
    radius = traits.Float(50, usedefault=True,
                          desc='radius in mm to calculate angular FDs, 50mm is the '
                               'default since it is used in Power et al. 2012')
    out_file = File('fd_power_2012.txt', usedefault=True, desc='output file name')
    out_figure = File('fd_power_2012.pdf', usedefault=True, desc='output figure name')
    series_tr = traits.Float(desc='repetition time in sec.')
    save_plot = traits.Bool(False, usedefault=True, desc='write FD plot')
    normalize = traits.Bool(False, usedefault=True, desc='calculate FD in mm/s')
    figdpi = traits.Int(100, usedefault=True, desc='output dpi for the FD plot')
    figsize = traits.Tuple(traits.Float(11.7), traits.Float(2.3), usedefault=True,
                           desc='output figure size')


class FramewiseDisplacementOutputSpec(TraitedSpec):
    out_file = File(desc='calculated FD per timestep')
    out_figure = File(desc='output image file')
    fd_average = traits.Float(desc='average FD')


class FramewiseDisplacement(BaseInterface):
    """
    Calculate the :abbr:`FD (framewise displacement)` as in [Power2012]_.
    This implementation reproduces the calculation in fsl_motion_outliers

    .. [Power2012] Power et al., Spurious but systematic correlations in functional
         connectivity MRI networks arise from subject motion, NeuroImage 59(3),
         2012. doi:`10.1016/j.neuroimage.2011.10.018
         <http://dx.doi.org/10.1016/j.neuroimage.2011.10.018>`_.


    """

    input_spec = FramewiseDisplacementInputSpec
    output_spec = FramewiseDisplacementOutputSpec

    references_ = [{
        'entry': BibTeX("""\
@article{power_spurious_2012,
    title = {Spurious but systematic correlations in functional connectivity {MRI} networks \
arise from subject motion},
    volume = {59},
    doi = {10.1016/j.neuroimage.2011.10.018},
    number = {3},
    urldate = {2016-08-16},
    journal = {NeuroImage},
    author = {Power, <NAME>. and Barnes, <NAME>. and Snyder, <NAME>. and Schlaggar, \
<NAME>. and Petersen, <NAME>.},
    year = {2012},
    pages = {2142--2154},
}
"""),
        'tags': ['method']
    }]

    def _run_interface(self, runtime):
        mpars = np.loadtxt(self.inputs.in_file)  # mpars is N_t x 6
        mpars = np.apply_along_axis(func1d=normalize_mc_params,
                                    axis=1, arr=mpars,
                                    source=self.inputs.parameter_source)
        diff = mpars[:-1, :6] - mpars[1:, :6]
        diff[:, 3:6] *= self.inputs.radius
        fd_res = np.abs(diff).sum(axis=1)

        self._results = {
            'out_file': op.abspath(self.inputs.out_file),
            'fd_average': float(fd_res.mean())
        }
        np.savetxt(self.inputs.out_file, fd_res, header='FramewiseDisplacement', comments='')

        if self.inputs.save_plot:
            tr = None
            if isdefined(self.inputs.series_tr):
                tr = self.inputs.series_tr

            if self.inputs.normalize and tr is None:
                IFLOG.warn('FD plot cannot be normalized if TR is not set')

            self._results['out_figure'] = op.abspath(self.inputs.out_figure)
            fig = plot_confound(fd_res, self.inputs.figsize, 'FD', units='mm',
                                series_tr=tr, normalize=self.inputs.normalize)
            fig.savefig(self._results['out_figure'], dpi=float(self.inputs.figdpi),
                        format=self.inputs.out_figure[-3:],
                        bbox_inches='tight')
            fig.clf()

        return runtime

    def _list_outputs(self):
        return self._results


class CompCorInputSpec(BaseInterfaceInputSpec):
    realigned_file = File(exists=True, mandatory=True,
                          desc='already realigned brain image (4D)')
    mask_files = InputMultiPath(File(exists=True),
                                desc=('One or more mask files that determines '
                                      'ROI (3D). When more that one file is '
                                      'provided `merge_method` or '
                                      '`merge_index` must be provided'))
    merge_method = traits.Enum('union', 'intersect', 'none', xor=['mask_index'],
                               requires=['mask_files'],
                               desc=('Merge method if multiple masks are '
                                     'present - `union` uses voxels included in'
                                     ' at least one input mask, `intersect` '
                                     'uses only voxels present in all input '
                                     'masks, `none` performs CompCor on '
                                     'each mask individually'))
    mask_index = traits.Range(low=0, xor=['merge_method'],
                              requires=['mask_files'],
                              desc=('Position of mask in `mask_files` to use - '
                                    'first is the default.'))
    components_file = traits.Str('components_file.txt', usedefault=True,
                           desc='Filename to store physiological components')
    num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL
    use_regress_poly = traits.Bool(True, usedefault=True,
                                   desc=('use polynomial regression '
                                         'pre-component extraction'))
    regress_poly_degree = traits.Range(low=1, default=1, usedefault=True,
                                       desc='the degree polynomial to use')
    header_prefix = traits.Str(desc=('the desired header for the output tsv '
                                     'file (one column). If undefined, will '
                                     'default to "CompCor"'))


class CompCorOutputSpec(TraitedSpec):
    components_file = File(exists=True,
                           desc='text file containing the noise components')


class CompCor(BaseInterface):
    """
    Interface with core CompCor computation, used in aCompCor and tCompCor

    Example
    -------

    >>> ccinterface = CompCor()
    >>> ccinterface.inputs.realigned_file = 'functional.nii'
    >>> ccinterface.inputs.mask_files = 'mask.nii'
    >>> ccinterface.inputs.num_components = 1
    >>> ccinterface.inputs.use_regress_poly = True
    >>> ccinterface.inputs.regress_poly_degree = 2

    """
    input_spec = CompCorInputSpec
    output_spec = CompCorOutputSpec
    references_ = [{'entry': BibTeX("@article{compcor_2007,"
                                    "title = {A component based noise correction method (CompCor) for BOLD and perfusion based},"
                                    "volume = {37},"
                                    "number = {1},"
                                    "doi = {10.1016/j.neuroimage.2007.04.042},"
                                    "urldate = {2016-08-13},"
                                    "journal = {NeuroImage},"
                                    "author = {<NAME> and <NAME> <NAME>.},"
                                    "year = {2007},"
                                    "pages = {90-101},}"
                                   ),
                    'tags': ['method', 'implementation']
                   }]

    def __init__(self, *args, **kwargs):
        ''' exactly the same as compcor except the header '''
        super(CompCor, self).__init__(*args, **kwargs)
        self._header = 'CompCor'

    def _run_interface(self, runtime):
        mask_images = []
        if isdefined(self.inputs.mask_files):
            mask_images = combine_mask_files(self.inputs.mask_files,
                                             self.inputs.merge_method,
                                             self.inputs.mask_index)

        degree = (self.inputs.regress_poly_degree if
                  self.inputs.use_regress_poly else 0)

        imgseries = nb.load(self.inputs.realigned_file,
                            mmap=NUMPY_MMAP)

        if len(imgseries.shape) != 4:
            raise ValueError('tCompCor expected a 4-D nifti file. Input {} has '
                            '{} dimensions (shape {})'.format(
                            self.inputs.realigned_file, len(imgseries.shape),
                            imgseries.shape))

        if len(mask_images) == 0:
            img = nb.Nifti1Image(np.ones(imgseries.shape[:3], dtype=np.bool),
                                 affine=imgseries.affine,
                                 header=imgseries.header)
            mask_images = [img]

        mask_images = self._process_masks(mask_images, imgseries.get_data())

        components = compute_noise_components(imgseries.get_data(),
                                              mask_images, degree,
                                              self.inputs.num_components)

        components_file = os.path.join(os.getcwd(), self.inputs.components_file)
        np.savetxt(components_file, components, fmt=b"%.10f", delimiter='\t',
                   header=self._make_headers(components.shape[1]), comments='')
        return runtime

    def _process_masks(self, mask_images, timeseries=None):
        return mask_images

    def _list_outputs(self):
        outputs = self._outputs().get()
        outputs['components_file'] = os.path.abspath(self.inputs.components_file)
        return outputs

    def _make_headers(self, num_col):
        headers = []
        header = self.inputs.header_prefix if \
            isdefined(self.inputs.header_prefix) else self._header
        for i in range(num_col):
            headers.append(header + '{:02d}'.format(i))
        return '\t'.join(headers)


class ACompCor(CompCor):
    """
    Anatomical compcor: for inputs and outputs, see CompCor.
    When the mask provided is an anatomical mask, then CompCor
    is equivalent to ACompCor.
    """

    def __init__(self, *args, **kwargs):
        ''' exactly the same as compcor except the header '''
        super(ACompCor, self).__init__(*args, **kwargs)
        self._header = 'aCompCor'


class TCompCorInputSpec(CompCorInputSpec):
    # and all the fields in CompCorInputSpec
    percentile_threshold = traits.Range(low=0., high=1., value=.02,
                                        exclude_low=True, exclude_high=True,
                                        usedefault=True, desc='the percentile '
                                        'used to select highest-variance '
                                        'voxels, represented by a number '
                                        'between 0 and 1, exclusive. By '
                                        'default, this value is set to .02. '
                                        'That is, the 2% of voxels '
                                        'with the highest variance are used.')


class TCompCorOutputSpec(CompCorOutputSpec):
    # and all the fields in CompCorOutputSpec
    high_variance_masks = OutputMultiPath(File(exists=True),
                                          desc=(("voxels exceeding the variance"
                                                 " threshold")))


class TCompCor(CompCor):
    """
    Interface for tCompCor. Computes a ROI mask based on variance of voxels.

    Example
    -------

    >>> ccinterface = TCompCor()
    >>> ccinterface.inputs.realigned_file = 'functional.nii'
    >>> ccinterface.inputs.mask_files = 'mask.nii'
    >>> ccinterface.inputs.num_components = 1
    >>> ccinterface.inputs.use_regress_poly = True
    >>> ccinterface.inputs.regress_poly_degree = 2
    >>> ccinterface.inputs.percentile_threshold = .03

    """

    input_spec = TCompCorInputSpec
    output_spec = TCompCorOutputSpec

    def __init__(self, *args, **kwargs):
        ''' exactly the same as compcor except the header '''
        super(TCompCor, self).__init__(*args, **kwargs)
        self._header = 'tCompCor'
        self._mask_files = []

    def _process_masks(self, mask_images, timeseries=None):
        out_images = []
        self._mask_files = []
        for i, img in enumerate(mask_images):
            mask = img.get_data().astype(np.bool)
            imgseries = timeseries[mask, :]
            imgseries = regress_poly(2, imgseries)
            tSTD = _compute_tSTD(imgseries, 0, axis=-1)
            threshold_std = np.percentile(tSTD, np.round(100. *
                           (1. - self.inputs.percentile_threshold)).astype(int))
            mask_data = np.zeros_like(mask)
            mask_data[mask != 0] = tSTD >= threshold_std
            out_image = nb.Nifti1Image(mask_data, affine=img.affine,
                                       header=img.header)

            # save mask
            mask_file = os.path.abspath('mask_{:03d}.nii.gz'.format(i))
            out_image.to_filename(mask_file)
            IFLOG.debug('tCompcor computed and saved mask of shape {} to '
                        'mask_file {}'.format(mask.shape, mask_file))
            self._mask_files.append(mask_file)
            out_images.append(out_image)
        return out_images

    def _list_outputs(self):
        outputs = super(TCompCor, self)._list_outputs()
        outputs['high_variance_masks'] = self._mask_files
        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')
    tsnr_file = File('tsnr.nii.gz', usedefault=True, hash_files=False,
                     desc='output tSNR file')
    mean_file = File('mean.nii.gz', usedefault=True, hash_files=False,
                     desc='output mean file')
    stddev_file = File('stdev.nii.gz', usedefault=True, hash_files=False,
                       desc='output tSNR file')
    detrended_file = File('detrend.nii.gz', usedefault=True, hash_files=False,
                          desc='input file after detrending')


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 _run_interface(self, runtime):
        img = nb.load(self.inputs.in_file[0], mmap=NUMPY_MMAP)
        header = img.header.copy()
        vollist = [nb.load(filename, mmap=NUMPY_MMAP) for filename in self.inputs.in_file]
        data = np.concatenate([vol.get_data().reshape(
            vol.shape[:3] + (-1,)) for vol in vollist], axis=3)
        data = np.nan_to_num(data)

        if data.dtype.kind == 'i':
            header.set_data_dtype(np.float32)
            data = data.astype(np.float32)

        if isdefined(self.inputs.regress_poly):
            data = regress_poly(self.inputs.regress_poly, data, remove_mean=False)
            img = nb.Nifti1Image(data, img.affine, header)
            nb.save(img, op.abspath(self.inputs.detrended_file))

        meanimg = np.mean(data, axis=3)
        stddevimg = np.std(data, axis=3)
        tsnr = np.zeros_like(meanimg)
        tsnr[stddevimg > 1.e-3] = meanimg[stddevimg > 1.e-3] / stddevimg[stddevimg > 1.e-3]
        img = nb.Nifti1Image(tsnr, img.affine, header)
        nb.save(img, op.abspath(self.inputs.tsnr_file))
        img = nb.Nifti1Image(meanimg, img.affine, header)
        nb.save(img, op.abspath(self.inputs.mean_file))
        img = nb.Nifti1Image(stddevimg, img.affine, header)
        nb.save(img, op.abspath(self.inputs.stddev_file))
        return runtime

    def _list_outputs(self):
        outputs = self._outputs().get()
        for k in ['tsnr_file', 'mean_file', 'stddev_file']:
            outputs[k] = op.abspath(getattr(self.inputs, k))

        if isdefined(self.inputs.regress_poly):
            outputs['detrended_file'] = op.abspath(self.inputs.detrended_file)
        return outputs


class NonSteadyStateDetectorInputSpec(BaseInterfaceInputSpec):
    in_file = File(exists=True, mandatory=True, desc='4D NIFTI EPI file')


class NonSteadyStateDetectorOutputSpec(TraitedSpec):
    n_volumes_to_discard = traits.Int(desc='Number of non-steady state volumes'
                                           'detected in the beginning of the scan.')


class NonSteadyStateDetector(BaseInterface):
    """
    Returns the number of non-steady state volumes detected at the beginning
    of the scan.
    """

    input_spec = NonSteadyStateDetectorInputSpec
    output_spec = NonSteadyStateDetectorOutputSpec

    def _run_interface(self, runtime):
        in_nii = nb.load(self.inputs.in_file)
        global_signal = in_nii.get_data()[:,:,:,:50].mean(axis=0).mean(axis=0).mean(axis=0)

        self._results = {
            'n_volumes_to_discard': _is_outlier(global_signal)
        }

        return runtime

    def _list_outputs(self):
        return self._results


def compute_dvars(in_file, in_mask, remove_zerovariance=False,
                  intensity_normalization=1000):
    """
    Compute the :abbr:`DVARS (D referring to temporal
    derivative of timecourses, VARS referring to RMS variance over voxels)`
    [Power2012]_.

    Particularly, the *standardized* :abbr:`DVARS (D referring to temporal
    derivative of timecourses, VARS referring to RMS variance over voxels)`
    [Nichols2013]_ are computed.

    .. [Nichols2013] <NAME>, `Notes on creating a standardized version of
         DVARS <http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-\
research/nichols/scripts/fsl/standardizeddvars.pdf>`_, 2013.

    .. note:: Implementation details

      Uses the implementation of the `Yule-Walker equations
      from nitime
      <http://nipy.org/nitime/api/generated/nitime.algorithms.autoregressive.html\
#nitime.algorithms.autoregressive.AR_est_YW>`_
      for the :abbr:`AR (auto-regressive)` filtering of the fMRI signal.

    :param numpy.ndarray func: functional data, after head-motion-correction.
    :param numpy.ndarray mask: a 3D mask of the brain
    :param bool output_all: write out all dvars
    :param str out_file: a path to which the standardized dvars should be saved.
    :return: the standardized DVARS

    """
    import numpy as np
    import nibabel as nb
    from nitime.algorithms import AR_est_YW
    import warnings

    func = nb.load(in_file, mmap=NUMPY_MMAP).get_data().astype(np.float32)
    mask = nb.load(in_mask, mmap=NUMPY_MMAP).get_data().astype(np.uint8)

    if len(func.shape) != 4:
        raise RuntimeError(
            "Input fMRI dataset should be 4-dimensional")

    idx = np.where(mask > 0)
    mfunc = func[idx[0], idx[1], idx[2], :]

    if intensity_normalization != 0:
        mfunc = (mfunc / np.median(mfunc)) * intensity_normalization

    # Robust standard deviation (we are using "lower" interpolation
    # because this is what FSL is doing
    func_sd = (np.percentile(mfunc, 75, axis=1, interpolation="lower") -
               np.percentile(mfunc, 25, axis=1, interpolation="lower")) / 1.349

    if remove_zerovariance:
        mfunc = mfunc[func_sd != 0, :]
        func_sd = func_sd[func_sd != 0]

    # Compute (non-robust) estimate of lag-1 autocorrelation
    ar1 = np.apply_along_axis(AR_est_YW, 1,
                              regress_poly(0, mfunc, remove_mean=True).astype(
                                  np.float32), 1)[:, 0]

    # Compute (predicted) standard deviation of temporal difference time series
    diff_sdhat = np.squeeze(np.sqrt(((1 - ar1) * 2).tolist())) * func_sd
    diff_sd_mean = diff_sdhat.mean()

    # Compute temporal difference time series
    func_diff = np.diff(mfunc, axis=1)

    # DVARS (no standardization)
    dvars_nstd = np.sqrt(np.square(func_diff).mean(axis=0))

    # standardization
    dvars_stdz = dvars_nstd / diff_sd_mean

    with warnings.catch_warnings():  # catch, e.g., divide by zero errors
        warnings.filterwarnings('error')

        # voxelwise standardization
        diff_vx_stdz = np.square(
            func_diff / np.array([diff_sdhat] * func_diff.shape[-1]).T)
        dvars_vx_stdz = np.sqrt(diff_vx_stdz.mean(axis=0))

    return (dvars_stdz, dvars_nstd, dvars_vx_stdz)


def plot_confound(tseries, figsize, name, units=None,
                  series_tr=None, normalize=False):
    """
    A helper function to plot :abbr:`fMRI (functional MRI)` confounds.

    """
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.gridspec import GridSpec
    from matplotlib.backends.backend_pdf import FigureCanvasPdf as FigureCanvas
    import seaborn as sns

    fig = plt.Figure(figsize=figsize)
    FigureCanvas(fig)
    grid = GridSpec(1, 2, width_ratios=[3, 1], wspace=0.025)
    grid.update(hspace=1.0, right=0.95, left=0.1, bottom=0.2)

    ax = fig.add_subplot(grid[0, :-1])
    if normalize and series_tr is not None:
        tseries /= series_tr

    ax.plot(tseries)
    ax.set_xlim((0, len(tseries)))
    ylabel = name
    if units is not None:
        ylabel += (' speed [{}/s]' if normalize else ' [{}]').format(units)
    ax.set_ylabel(ylabel)

    xlabel = 'Frame #'
    if series_tr is not None:
        xlabel = 'Frame # ({} sec TR)'.format(series_tr)
    ax.set_xlabel(xlabel)
    ylim = ax.get_ylim()

    ax = fig.add_subplot(grid[0, -1])
    sns.distplot(tseries, vertical=True, ax=ax)
    ax.set_xlabel('Frames')
    ax.set_ylim(ylim)
    ax.set_yticklabels([])
    return fig


def is_outlier(points, thresh=3.5):
    """
    Returns a boolean array with True if points are outliers and False
    otherwise.

    :param nparray points: an numobservations by numdimensions numpy array of observations
    :param float thresh: the modified z-score to use as a threshold. Observations with
        a modified z-score (based on the median absolute deviation) greater
        than this value will be classified as outliers.

    :return: A bolean mask, of size numobservations-length array.

    .. note:: References

        <NAME> and <NAME> (1993), "Volume 16: How to Detect and
        Handle Outliers", The ASQC Basic References in Quality Control:
        Statistical Techniques, <NAME>, Ph.D., Editor.

    """
    if len(points.shape) == 1:
        points = points[:, None]
    median = np.median(points, axis=0)
    diff = np.sum((points - median) ** 2, axis=-1)
    diff = np.sqrt(diff)
    med_abs_deviation = np.median(diff)

    modified_z_score = 0.6745 * diff / med_abs_deviation

    timepoints_to_discard = 0
    for i in range(len(modified_z_score)):
        if modified_z_score[i] <= thresh:
            break
        else:
            timepoints_to_discard += 1

    return timepoints_to_discard


def regress_poly(degree, data, remove_mean=True, axis=-1):
    """
    Returns data with degree polynomial regressed out.

    :param bool remove_mean: whether or not demean data (i.e. degree 0),
    :param int axis: numpy array axes along which regression is performed

    """
    IFLOG.debug('Performing polynomial regression on data of shape ' + str(data.shape))

    datashape = data.shape
    timepoints = datashape[axis]

    # Rearrange all voxel-wise time-series in rows
    data = data.reshape((-1, timepoints))

    # Generate design matrix
    X = np.ones((timepoints, 1)) # quick way to calc degree 0
    for i in range(degree):
        polynomial_func = Legendre.basis(i + 1)
        value_array = np.linspace(-1, 1, timepoints)
        X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis]))

    # Calculate coefficients
    betas = np.linalg.pinv(X).dot(data.T)

    # Estimation
    if remove_mean:
        datahat = X.dot(betas).T
    else: # disregard the first layer of X, which is degree 0
        datahat = X[:, 1:].dot(betas[1:, ...]).T
    regressed_data = data - datahat

    # Back to original shape
    return regressed_data.reshape(datashape)


def combine_mask_files(mask_files, mask_method=None, mask_index=None):
    """Combines input mask files into a single nibabel image

    A helper function for CompCor

    mask_files: a list
        one or more binary mask files
    mask_method: enum ('union', 'intersect', 'none')
        determines how to combine masks
    mask_index: an integer
        determines which file to return (mutually exclusive with mask_method)

    returns: a list of nibabel images
    """

    if isdefined(mask_index) or not isdefined(mask_method):
        if not isdefined(mask_index):
            if len(mask_files) == 1:
                mask_index = 0
            else:
                raise ValueError(('When more than one mask file is provided, '
                                  'one of merge_method or mask_index must be '
                                  'set'))
        if mask_index < len(mask_files):
            mask = nb.load(mask_files[mask_index], mmap=NUMPY_MMAP)
            return [mask]
        raise ValueError(('mask_index {0} must be less than number of mask '
                          'files {1}').format(mask_index, len(mask_files)))
    masks = []
    if mask_method == 'none':
        for filename in mask_files:
            masks.append(nb.load(filename, mmap=NUMPY_MMAP))
        return masks

    if mask_method == 'union':
        mask = None
        for filename in mask_files:
            img = nb.load(filename, mmap=NUMPY_MMAP)
            if mask is None:
                mask = img.get_data() > 0
            np.logical_or(mask, img.get_data() > 0, mask)
        img = nb.Nifti1Image(mask, img.affine, header=img.header)
        return [img]

    if mask_method == 'intersect':
        mask = None
        for filename in mask_files:
            img = nb.load(filename, mmap=NUMPY_MMAP)
            if mask is None:
                mask = img.get_data() > 0
            np.logical_and(mask, img.get_data() > 0, mask)
        img = nb.Nifti1Image(mask, img.affine, header=img.header)
        return [img]


def compute_noise_components(imgseries, mask_images, degree, num_components):
    """Compute the noise components from the imgseries for each mask

    imgseries: a nibabel img
    mask_images: a list of nibabel images
    degree: order of polynomial used to remove trends from the timeseries
    num_components: number of noise components to return

    returns:

    components: a numpy array

    """
    components = None
    for img in mask_images:
        mask = img.get_data().astype(np.bool)
        if imgseries.shape[:3] != mask.shape:
            raise ValueError('Inputs for CompCor, timeseries and mask, '
                             'do not have matching spatial dimensions '
                             '({} and {}, respectively)'.format(
                imgseries.shape[:3], mask.shape))

        voxel_timecourses = imgseries[mask, :]

        # Zero-out any bad values
        voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0

        # from paper:
        # "The constant and linear trends of the columns in the matrix M were
        # removed [prior to ...]"
        voxel_timecourses = regress_poly(degree, voxel_timecourses)

        # "Voxel time series from the noise ROI (either anatomical or tSTD) were
        # placed in a matrix M of size Nxm, with time along the row dimension
        # and voxels along the column dimension."
        M = voxel_timecourses.T

        # "[... were removed] prior to column-wise variance normalization."
        M = M / _compute_tSTD(M, 1.)

        # "The covariance matrix C = MMT was constructed and decomposed into its
        # principal components using a singular value decomposition."
        u, _, _ = linalg.svd(M, full_matrices=False)
        if components is None:
            components = u[:, :num_components]
        else:
            components = np.hstack((components,
                                    u[:, :num_components]))
    if components is None and num_components > 0:
        raise ValueError('No components found')
    return components


def _compute_tSTD(M, x, axis=0):
    stdM = np.std(M, axis=axis)
    # set bad values to x
    stdM[stdM == 0] = x
    stdM[np.isnan(stdM)] = x
    return stdM
# -*- coding: utf-8 -*- # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: Algorithms to compute confounds in :abbr:`fMRI (functional MRI)`      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) Computes the DVARS. \ @techreport{nichols_notes_2013,     address = {Coventry, UK},     title = {Notes on {Creating} a {Standardized} {Version} of {DVARS}},     url = {http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-\ research/nichols/scripts/fsl/standardizeddvars.pdf},     urldate = {2016-08-16},     institution = {University of Warwick},     author = {<NAME>},     year = {2013} } \ @article{power_spurious_2012,     title = {Spurious but systematic correlations in functional connectivity {MRI} networks \ arise from subject motion},     volume = {59},     doi = {10.1016/j.neuroimage.2011.10.018},     number = {3},     urldate = {2016-08-16},     journal = {NeuroImage},     author = {Power, <NAME>. and Barnes, <NAME>. and Snyder, <NAME>. and Schlaggar, \ <NAME>. and Petersen, <NAME>.},     year = {2012},     pages = {2142--2154}, } Calculate the :abbr:`FD (framewise displacement)` as in [Power2012]_.     This implementation reproduces the calculation in fsl_motion_outliers      .. [Power2012] Power et al., Spurious but systematic correlations in functional          connectivity MRI networks arise from subject motion, NeuroImage 59(3),          2012. doi:`10.1016/j.neuroimage.2011.10.018          <http://dx.doi.org/10.1016/j.neuroimage.2011.10.018>`_. \ @article{power_spurious_2012,     title = {Spurious but systematic correlations in functional connectivity {MRI} networks \ arise from subject motion},     volume = {59},     doi = {10.1016/j.neuroimage.2011.10.018},     number = {3},     urldate = {2016-08-16},     journal = {NeuroImage},     author = {Power, <NAME>. and Barnes, <NAME>. and Snyder, <NAME>. and Schlaggar, \ <NAME>. and Petersen, <NAME>.},     year = {2012},     pages = {2142--2154}, } # mpars is N_t x 6 # 6 for BOLD, 4 for ASL Interface with core CompCor computation, used in aCompCor and tCompCor      Example     -------      >>> ccinterface = CompCor()     >>> ccinterface.inputs.realigned_file = 'functional.nii'     >>> ccinterface.inputs.mask_files = 'mask.nii'     >>> ccinterface.inputs.num_components = 1     >>> ccinterface.inputs.use_regress_poly = True     >>> ccinterface.inputs.regress_poly_degree = 2 exactly the same as compcor except the header Anatomical compcor: for inputs and outputs, see CompCor.     When the mask provided is an anatomical mask, then CompCor     is equivalent to ACompCor. exactly the same as compcor except the header # and all the fields in CompCorInputSpec # and all the fields in CompCorOutputSpec Interface for tCompCor. Computes a ROI mask based on variance of voxels.      Example     -------      >>> ccinterface = TCompCor()     >>> ccinterface.inputs.realigned_file = 'functional.nii'     >>> ccinterface.inputs.mask_files = 'mask.nii'     >>> ccinterface.inputs.num_components = 1     >>> ccinterface.inputs.use_regress_poly = True     >>> ccinterface.inputs.regress_poly_degree = 2     >>> ccinterface.inputs.percentile_threshold = .03 exactly the same as compcor except the header # save mask 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 Returns the number of non-steady state volumes detected at the beginning     of the scan. Compute the :abbr:`DVARS (D referring to temporal     derivative of timecourses, VARS referring to RMS variance over voxels)`     [Power2012]_.      Particularly, the *standardized* :abbr:`DVARS (D referring to temporal     derivative of timecourses, VARS referring to RMS variance over voxels)`     [Nichols2013]_ are computed.      .. [Nichols2013] <NAME>, `Notes on creating a standardized version of          DVARS <http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-\ research/nichols/scripts/fsl/standardizeddvars.pdf>`_, 2013.      .. note:: Implementation details        Uses the implementation of the `Yule-Walker equations       from nitime       <http://nipy.org/nitime/api/generated/nitime.algorithms.autoregressive.html\ #nitime.algorithms.autoregressive.AR_est_YW>`_       for the :abbr:`AR (auto-regressive)` filtering of the fMRI signal.      :param numpy.ndarray func: functional data, after head-motion-correction.     :param numpy.ndarray mask: a 3D mask of the brain     :param bool output_all: write out all dvars     :param str out_file: a path to which the standardized dvars should be saved.     :return: the standardized DVARS # Robust standard deviation (we are using "lower" interpolation # because this is what FSL is doing # Compute (non-robust) estimate of lag-1 autocorrelation # Compute (predicted) standard deviation of temporal difference time series # Compute temporal difference time series # DVARS (no standardization) # standardization # catch, e.g., divide by zero errors # voxelwise standardization A helper function to plot :abbr:`fMRI (functional MRI)` confounds. #' # ({} sec TR)'.format(series_tr) Returns a boolean array with True if points are outliers and False     otherwise.      :param nparray points: an numobservations by numdimensions numpy array of observations     :param float thresh: the modified z-score to use as a threshold. Observations with         a modified z-score (based on the median absolute deviation) greater         than this value will be classified as outliers.      :return: A bolean mask, of size numobservations-length array.      .. note:: References          <NAME> and <NAME> (1993), "Volume 16: How to Detect and         Handle Outliers", The ASQC Basic References in Quality Control:         Statistical Techniques, <NAME>, Ph.D., Editor. Returns data with degree polynomial regressed out.      :param bool remove_mean: whether or not demean data (i.e. degree 0),     :param int axis: numpy array axes along which regression is performed # Rearrange all voxel-wise time-series in rows # Generate design matrix # quick way to calc degree 0 # Calculate coefficients # Estimation # disregard the first layer of X, which is degree 0 # Back to original shape Combines input mask files into a single nibabel image      A helper function for CompCor      mask_files: a list         one or more binary mask files     mask_method: enum ('union', 'intersect', 'none')         determines how to combine masks     mask_index: an integer         determines which file to return (mutually exclusive with mask_method)      returns: a list of nibabel images Compute the noise components from the imgseries for each mask      imgseries: a nibabel img     mask_images: a list of nibabel images     degree: order of polynomial used to remove trends from the timeseries     num_components: number of noise components to return      returns:      components: a numpy array # Zero-out any bad values # from paper: # "The constant and linear trends of the columns in the matrix M were # removed [prior to ...]" # "Voxel time series from the noise ROI (either anatomical or tSTD) were # placed in a matrix M of size Nxm, with time along the row dimension # and voxels along the column dimension." # "[... were removed] prior to column-wise variance normalization." # "The covariance matrix C = MMT was constructed and decomposed into its # principal components using a singular value decomposition." # set bad values to x