<commit_before># Authors: <NAME> <<EMAIL>>
#          <NAME> <<EMAIL>>
#          <NAME> <<EMAIL>>
#
# License: BSD (3-clause)

import os
import copy
from math import ceil
import numpy as np
from scipy import linalg, sparse
from scipy.sparse import csr_matrix, coo_matrix
import warnings

from .filter import resample
from .parallel import parallel_func
from .surface import (read_surface, _get_ico_surface, read_morph_map,
                      _compute_nearest)
from .utils import (get_subjects_dir, _check_subject,
                    _check_pandas_index_arguments, _check_pandas_installed,
                    logger, verbose)
from .viz import plot_source_estimates
from . fixes import in1d


def _read_stc(filename):
    """ Aux Function
    """
    fid = open(filename, 'rb')

    stc = dict()

    fid.seek(0, 2)  # go to end of file
    file_length = fid.tell()
    fid.seek(0, 0)  # go to beginning of file

    # read tmin in ms
    stc['tmin'] = float(np.fromfile(fid, dtype=">f4", count=1))
    stc['tmin'] /= 1000.0

    # read sampling rate in ms
    stc['tstep'] = float(np.fromfile(fid, dtype=">f4", count=1))
    stc['tstep'] /= 1000.0

    # read number of vertices/sources
    vertices_n = int(np.fromfile(fid, dtype=">u4", count=1))

    # read the source vector
    stc['vertices'] = np.fromfile(fid, dtype=">u4", count=vertices_n)

    # read the number of timepts
    data_n = int(np.fromfile(fid, dtype=">u4", count=1))

    if (vertices_n and  # vertices_n can be 0 (empty stc)
            ((file_length / 4 - 4 - vertices_n) % (data_n * vertices_n)) != 0):
        raise ValueError('incorrect stc file size')

    # read the data matrix
    stc['data'] = np.fromfile(fid, dtype=">f4", count=vertices_n * data_n)
    stc['data'] = stc['data'].reshape([data_n, vertices_n]).T

    # close the file
    fid.close()
    return stc


def _write_stc(filename, tmin, tstep, vertices, data):
    """Write an STC file

    Parameters
    ----------
    filename : string
        The name of the STC file.
    tmin : float
        The first time point of the data in seconds.
    tstep : float
        Time between frames in seconds.
    vertices : array of integers
        Vertex indices (0 based).
    data : 2D array
        The data matrix (nvert * ntime).
    """
    fid = open(filename, 'wb')

    # write start time in ms
    fid.write(np.array(1000 * tmin, dtype='>f4').tostring())
    # write sampling rate in ms
    fid.write(np.array(1000 * tstep, dtype='>f4').tostring())
    # write number of vertices
    fid.write(np.array(vertices.shape[0], dtype='>u4').tostring())
    # write the vertex indices
    fid.write(np.array(vertices, dtype='>u4').tostring())

    # write the number of timepts
    fid.write(np.array(data.shape[1], dtype='>u4').tostring())
    #
    # write the data
    #
    fid.write(np.array(data.T, dtype='>f4').tostring())

    # close the file
    fid.close()


def _read_3(fid):
    """ Read 3 byte integer from file
    """
    data = np.fromfile(fid, dtype=np.uint8, count=3).astype(np.int32)

    out = np.left_shift(data[0], 16) + np.left_shift(data[1], 8) + data[2]

    return out


def _read_w(filename):
    """Read a w file and return as dict

    w files contain activations or source reconstructions for a single time
    point.

    Parameters
    ----------
    filename : string
        The name of the w file.

    Returns
    -------
    data: dict
        The w structure. It has the following keys:
           vertices       vertex indices (0 based)
           data           The data matrix (nvert long)
    """

    fid = open(filename, 'rb')

    # skip first 2 bytes
    fid.read(2)

    # read number of vertices/sources (3 byte integer)
    vertices_n = int(_read_3(fid))

    vertices = np.zeros((vertices_n), dtype=np.int32)
    data = np.zeros((vertices_n), dtype=np.float32)

    # read the vertices and data
    for i in range(vertices_n):
        vertices[i] = _read_3(fid)
        data[i] = np.fromfile(fid, dtype='>f4', count=1)

    w = dict()
    w['vertices'] = vertices
    w['data'] = data

    # close the file
    fid.close()
    return w


def _write_3(fid, val):
    """ Write 3 byte integer to file
    """

    f_bytes = np.zeros((3), dtype=np.uint8)

    f_bytes[0] = (val >> 16) & 255
    f_bytes[1] = (val >> 8) & 255
    f_bytes[2] = val & 255

    fid.write(f_bytes.tostring())


def _write_w(filename, vertices, data):
    """Read a w file

    w files contain activations or source reconstructions for a single time
    point.

    Parameters
    ----------
    filename: string
        The name of the w file.
    vertices: array of int
        Vertex indices (0 based).
    data: 1D array
        The data array (nvert).
    """

    assert(len(vertices) == len(data))

    fid = open(filename, 'wb')

    # write 2 zero bytes
    fid.write(np.zeros((2), dtype=np.uint8).tostring())

    # write number of vertices/sources (3 byte integer)
    vertices_n = len(vertices)
    _write_3(fid, vertices_n)

    # write the vertices and data
    for i in range(vertices_n):
        _write_3(fid, vertices[i])
        #XXX: without float() endianness is wrong, not sure why
        fid.write(np.array(float(data[i]), dtype='>f4').tostring())

    # close the file
    fid.close()


def read_source_estimate(fname, subject=None):
    """Read a soure estimate object

    Parameters
    ----------
    fname : str
        Path to (a) source-estimate file(s).
    subject : str | None
        Name of the subject the source estimate(s) is (are) from.
        It is good practice to set this attribute to avoid combining
        incompatible labels and SourceEstimates (e.g., ones from other
        subjects). Note that due to file specification limitations, the
        subject name isn't saved to or loaded from files written to disk.

    Returns
    -------
    stc : SourceEstimate | VolSourceEstimate
        The soure estimate object loaded from file.

    Notes
    -----
     - for volume source estimates, ``fname`` should provide the path to a
       single file named '*-vl.stc` or '*-vol.stc'
     - for surface source estimates, ``fname`` should either provide the
       path to the file corresponding to a single hemisphere ('*-lh.stc',
       '*-rh.stc') or only specify the asterisk part in these patterns. In any
       case, the function expects files for both hemisphere with names
       following this pattern.
     - for single time point .w files, ``fname`` should follow the same
       pattern as for surface estimates, except that files are named
       '*-lh.w' and '*-rh.w'.
    """
    fname_arg = fname

    # make sure corresponding file(s) can be found
    ftype = None
    if os.path.exists(fname):
        if fname.endswith('-vl.stc') or fname.endswith('-vol.stc') or \
                fname.endswith('-vl.w') or fname.endswith('-vol.w'):
            ftype = 'volume'
        elif fname.endswith('.stc'):
            ftype = 'surface'
            if fname.endswith(('-lh.stc', '-rh.stc')):
                fname = fname[:-7]
            else:
                err = ("Invalid .stc filename: %r; needs to end with "
                       "hemisphere tag ('...-lh.stc' or '...-rh.stc')"
                       % fname)
                raise IOError(err)
        elif fname.endswith('.w'):
            ftype = 'w'
            if fname.endswith(('-lh.w', '-rh.w')):
                fname = fname[:-5]
            else:
                err = ("Invalid .w filename: %r; needs to end with "
                       "hemisphere tag ('...-lh.w' or '...-rh.w')"
                       % fname)
                raise IOError(err)

    if ftype is not 'volume':
        stc_exist = map(os.path.exists, (fname + '-rh.stc', fname + '-lh.stc'))
        w_exist = map(os.path.exists, (fname + '-rh.w', fname + '-lh.w'))
        if all(stc_exist) and (ftype is not 'w'):
            ftype = 'surface'
        elif all(w_exist):
            ftype = 'w'
        elif any(stc_exist) or any(w_exist):
            raise IOError("Hemisphere missing for %r" % fname_arg)
        else:
            raise IOError("SourceEstimate File(s) not found for: %r"
                          % fname_arg)

    # read the files
    if ftype == 'volume':  # volume source space
        if fname.endswith('.stc'):
            kwargs = _read_stc(fname)
        elif fname.endswith('.w'):
            kwargs = _read_w(fname)
            kwargs['data'] = kwargs['data'][:, np.newaxis]
            kwargs['tmin'] = 0.0
            kwargs['tstep'] = 0.0
        else:
            raise IOError('Volume source estimate must end with .stc or .w')
    elif ftype == 'surface':  # stc file with surface source spaces
        lh = _read_stc(fname + '-lh.stc')
        rh = _read_stc(fname + '-rh.stc')
        assert lh['tmin'] == rh['tmin']
        assert lh['tstep'] == rh['tstep']
        kwargs = lh.copy()
        kwargs['data'] = np.r_[lh['data'], rh['data']]
        kwargs['vertices'] = [lh['vertices'], rh['vertices']]
    elif ftype == 'w':  # w file with surface source spaces
        lh = _read_w(fname + '-lh.w')
        rh = _read_w(fname + '-rh.w')
        kwargs = lh.copy()
        kwargs['data'] = np.atleast_2d(np.r_[lh['data'], rh['data']]).T
        kwargs['vertices'] = [lh['vertices'], rh['vertices']]
        # w files only have a single time point
        kwargs['tmin'] = 0.0
        kwargs['tstep'] = 1.0

    if ftype != 'volume':
        # Make sure the vertices are ordered
        vertices = kwargs['vertices']
        if any([np.any(np.diff(v.astype(int)) <= 0) for v in vertices]):
            sidx = [np.argsort(verts) for verts in vertices]
            vertices = [verts[idx] for verts, idx in zip(vertices, sidx)]
            data = kwargs['data'][np.r_[sidx[0], len(sidx[0]) + sidx[1]]]
            kwargs['vertices'] = vertices
            kwargs['data'] = data

    kwargs['subject'] = subject

    if ftype == 'volume':
        stc = VolSourceEstimate(**kwargs)
    else:
        stc = SourceEstimate(**kwargs)

    return stc


def _make_stc(data, vertices, tmin=None, tstep=None, subject=None):
    """Helper function to generate either a surface or volume source estimate
    """

    if isinstance(vertices, list) and len(vertices) == 2:
        # make a surface source estimate
        stc = SourceEstimate(data, vertices=vertices, tmin=tmin, tstep=tstep,
                             subject=subject)
    elif isinstance(vertices, np.ndarray) or isinstance(vertices, list)\
            and len(vertices) == 1:
        stc = VolSourceEstimate(data, vertices=vertices, tmin=tmin,
                                tstep=tstep, subject=subject)
    else:
        raise ValueError('vertices has to be either a list with one or two '
                         'arrays or an array')
    return stc


class _NotifyArray(np.ndarray):
    """Array class that executes a callback when it is modified
    """
    def __new__(cls, input_array, modify_callback=None):
        obj = np.asarray(input_array).view(cls)
        obj.modify_callback = modify_callback
        return obj

    def __array_finalize__(self, obj):
        if obj is None:
            # an empty constructor was used
            return

        # try to copy the callback
        self.modify_callback = getattr(obj, 'modify_callback', None)

    def _modified_(self):
        """Execute the callback if it is set"""
        if self.modify_callback is not None:
            self.modify_callback()

    def __getattribute__(self, name):
        # catch ndarray methods that modify the array inplace
        if name in ['fill', 'itemset', 'resize', 'sort']:
            self._modified_()

        return object.__getattribute__(self, name)

    def __setitem__(self, item, value):
        self._modified_()
        np.ndarray.__setitem__(self, item, value)

    def __array_wrap__(self, out_arr, context=None):
        # this method is called whenever a numpy ufunc (+, +=..) is called
        # the last entry in context is the array that receives the result
        if (context is not None and len(context[1]) == 3
                and context[1][2] is self):
            self._modified_()

        return np.ndarray.__array_wrap__(self, out_arr, context)


def _verify_source_estimate_compat(a, b):
    """Make sure two SourceEstimates are compatible for arith. operations"""
    compat = False
    if len(a.vertno) == len(b.vertno):
        if all([np.array_equal(av, vv) for av, vv in zip(a.vertno, b.vertno)]):
            compat = True
    if not compat:
        raise ValueError('Cannot combine SourceEstimates that do not have the '
                         'same vertices. Consider using stc.expand().')
    if a.subject != b.subject:
        raise ValueError('source estimates do not have the same subject '
                         'names, "%s" and "%s"' % (a.name, b.name))


class _BaseSourceEstimate(object):
    """Abstract base class for source estimates

    Parameters
    ----------
    data : array of shape (n_dipoles, n_times) | 2-tuple (kernel, sens_data)
        The data in source space. The data can either be a single array or
        a tuple with two arrays: "kernel" shape (n_vertices, n_sensors) and
        "sens_data" shape (n_sensors, n_times). In this case, the source
        space data corresponds to "numpy.dot(kernel, sens_data)".
    vertices : array | list of two arrays
        Vertex numbers corresponding to the data.
    tmin : scalar
        Time point of the first sample in data.
    tstep : scalar
        Time step between successive samples in data.
    subject : str | None
        The subject name. While not necessary, it is safer to set the
        subject parameter to avoid analysis errors.
    verbose : bool, str, int, or None
        If not None, override default verbose level (see mne.verbose).

    Attributes
    ----------
    subject : str | None
        The subject name.
    times : array of shape (n_times,)
        The time vector.
    vertices : array or list of arrays of shape (n_dipoles,)
        The indices of the dipoles in the different source spaces. Can
        be an array if there is only one source space (e.g., for volumes).
    data : array of shape (n_dipoles, n_times)
        The data in source space.
    shape : tuple
        The shape of the data. A tuple of int (n_dipoles, n_times).
    """
    @verbose
    def __init__(self, data, vertices=None, tmin=None, tstep=None,
                 subject=None, verbose=None):
        kernel, sens_data = None, None
        if isinstance(data, tuple):
            if len(data) != 2:
                raise ValueError('If data is a tuple it has to be length 2')
            kernel, sens_data = data
            data = None
            if kernel.shape[1] != sens_data.shape[0]:
                raise ValueError('kernel and sens_data have invalid '
                                 'dimensions')

        if isinstance(vertices, list):
            if not (len(vertices) == 2 or len(vertices) == 1) or \
                    not all([isinstance(v, np.ndarray) for v in vertices]):
                raise ValueError('Vertices, if a list, must contain one or '
                                 'two numpy arrays')

            if any([np.any(np.diff(v.astype(int)) <= 0) for v in vertices]):
                raise ValueError('Vertices must be ordered in increasing '
                                 'order.')

            n_src = sum([len(v) for v in vertices])

            if len(vertices) == 1:
                vertices = vertices[0]
        elif isinstance(vertices, np.ndarray):
            n_src = len(vertices)
        else:
            raise ValueError('Vertices must be a list or numpy array')

        # safeguard the user against doing something silly
        if data is not None and data.shape[0] != n_src:
            raise ValueError('Number of vertices (%i) and stc.shape[0] (%i) '
                             'must match' % (n_src, data.shape[0]))

        self._data = data
        self.tmin = tmin
        self.tstep = tstep
        self.vertno = vertices
        self.verbose = verbose
        self._kernel = kernel
        self._sens_data = sens_data
        self.times = None
        self._update_times()
        self.subject = _check_subject(None, subject, False)

    def _remove_kernel_sens_data_(self):
        """Remove kernel and sensor space data

        Note: self._data is also computed if it is None
        """
        if self._kernel is not None or self._sens_data is not None:
            # we can no longer use the kernel and sens_data
            logger.info('STC data modified: removing kernel and sensor data')
            if self._data is None:
                self._data = np.dot(self._kernel, self._sens_data)
            self._kernel = None
            self._sens_data = None

    def crop(self, tmin=None, tmax=None):
        """Restrict SourceEstimate to a time interval

        Parameters
        ----------
        tmin : float or None
            The first time point in seconds. If None the first present is used.
        tmax : float or None
            The last time point in seconds. If None the last present is used.
        """
        mask = np.ones(len(self.times), dtype=np.bool)
        if tmax is not None:
            mask = mask & (self.times <= tmax)
        if tmin is not None:
            mask = mask & (self.times >= tmin)
            self.tmin = tmin

        if self._kernel is not None and self._sens_data is not None:
            self._sens_data = self._sens_data[:, mask]
            self._data = None  # will be recomputed when data is accessed
        else:
            self._data = self._data[:, mask]

        self._update_times()
        return self  # return self for chaining methods

    @verbose
    def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1,
                 verbose=None):
        """Resample data

        Parameters
        ----------
        sfreq : float
            New sample rate to use.
        npad : int
            Amount to pad the start and end of the data.
        window : string or tuple
            Window to use in resampling. See scipy.signal.resample.
        n_jobs : int
            Number of jobs to run in parallel.
        verbose : bool, str, int, or None
            If not None, override default verbose level (see mne.verbose).
            Defaults to self.verbose.

        Notes
        -----
        For some data, it may be more accurate to use npad=0 to reduce
        artifacts. This is dataset dependent -- check your data!

        Note that the sample rate of the original data is inferred from tstep.
        """
        # resampling in sensor instead of source space gives a somewhat
        # different result, so we don't allow it
        self._remove_kernel_sens_data_()

        o_sfreq = 1.0 / self.tstep
        self._data = resample(self._data, sfreq, o_sfreq, npad, n_jobs=n_jobs)

        # adjust indirectly affected variables
        self.tstep = 1.0 / sfreq
        self._update_times()

    @property
    def data(self):
        if self._data is None:
            # compute the solution the first time the data is accessed
            # return a "notify array", so we can later remove the kernel
            # and sensor data if the user modifies self._data
            self._data = _NotifyArray(np.dot(self._kernel, self._sens_data),
                modify_callback=self._remove_kernel_sens_data_)
        return self._data

    @property
    def shape(self):
        if self._data is not None:
            return self._data.shape
        return (self._kernel.shape[0], self._sens_data.shape[1])

    def _update_times(self):
        """Update the times attribute after changing tmin, tmax, or tstep"""
        self.times = self.tmin + (self.tstep * np.arange(self.shape[1]))

    def __add__(self, a):
        stc = copy.deepcopy(self)
        stc += a
        return stc

    def __iadd__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data += a.data
        else:
            self._data += a
        return self

    def mean(self):
        """Make a summary stc file with mean power between tmin and tmax.

        Returns
        -------
        stc : instance of SourceEstimate
            The modified stc (note: method operates inplace).
        """
        data = self.data
        tmax = self.tmin + self.tstep * data.shape[1]
        tmin = (self.tmin + tmax) / 2.
        tstep = tmax - self.tmin
        mean_stc = SourceEstimate(self.data.mean(axis=1)[:, np.newaxis],
                                  vertices=self.vertno, tmin=tmin,
                                  tstep=tstep, subject=self.subject)
        return mean_stc

    def __sub__(self, a):
        stc = copy.deepcopy(self)
        stc -= a
        return stc

    def __isub__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data -= a.data
        else:
            self._data -= a
        return self

    def __div__(self, a):
        stc = copy.deepcopy(self)
        stc /= a
        return stc

    def __idiv__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data /= a.data
        else:
            self._data /= a
        return self

    def __mul__(self, a):
        stc = copy.deepcopy(self)
        stc *= a
        return stc

    def __imul__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data *= a.data
        else:
            self._data *= a
        return self

    def __pow__(self, a):
        stc = copy.deepcopy(self)
        stc **= a
        return stc

    def __ipow__(self, a):
        self._remove_kernel_sens_data_()
        self._data **= a
        return self

    def __radd__(self, a):
        return self + a

    def __rsub__(self, a):
        return self - a

    def __rmul__(self, a):
        return self * a

    def __rdiv__(self, a):
        return self / a

    def __neg__(self):
        stc = copy.deepcopy(self)
        stc._remove_kernel_sens_data_()
        stc._data *= -1
        return stc

    def __pos__(self):
        return self

    def sqrt(self):
        """Return copy of SourceEstimate with sqrt(data)."""
        return self ** (0.5)

    def copy(self):
        """Return copy of SourceEstimate instance"""
        return copy.deepcopy(self)

    def bin(self, width, tstart=None, tstop=None, func=np.mean):
        """Returns a SourceEstimate object with data summarized over time bins

        Time bins of ``width`` seconds. This method is intended for
        visualization only. No filter is applied to the data before binning,
        making the method inappropriate as a tool for downsampling data.

        Parameters
        ----------
        width : scalar
            Width of the individual bins in seconds.
        func : callable
            Function that is applied to summarize the data. Needs to accept a
            numpy.array as first input and an ``axis`` keyword argument.
        tstart : scalar | None
            Time point where the first bin starts. The default is the first
            time point of the stc.
        tstop : scalar | None
            Last possible time point contained in a bin (if the last bin would
            be shorter than width it is dropped). The default is the last time
            point of the stc.

        Returns
        -------
        stc : instance of SourceEstimate
            The binned SourceEstimate.
        """
        if tstart is None:
            tstart = self.tmin
        if tstop is None:
            tstop = self.times[-1]

        times = np.arange(tstart, tstop + self.tstep, width)
        nv, _ = self.shape
        nt = len(times) - 1
        data = np.empty((nv, nt), dtype=self.data.dtype)
        for i in xrange(nt):
            idx = (self.times >= times[i]) & (self.times < times[i + 1])
            data[:, i] = func(self.data[:, idx], axis=1)

        tmin = times[0] + width / 2.
        stc = _make_stc(data, vertices=self.vertno,
                        tmin=tmin, tstep=width, subject=self.subject)
        return stc

    def transform_data(self, transform_fun, fun_args=None,
                       idx=None, tmin_idx=None, tmax_idx=None, **kwargs):
        """Get data after a linear (time) transform has been applied

        The transorm is applied to each source time course independently.


        Parameters
        ----------
        transform_fun : callable
            The transform to be applied. The first parameter of the function
            is the input data. The first return value is the transformed
            data, remaining outputs are ignored. The first dimension of the
            transformed data has to be the same as the first dimension of the
            input data.
        fun_args : tuple | None
            Additional parameters to be passed to transform_fun.
        idx : array | None
            Indicices of source time courses for which to compute transform.
            If None, all time courses are used.
        tmin_idx : int | None
            Index of first time point to include. If None, the index of the
            first time point is used.
        tmax_idx : int | None
            Index of the first time point not to include. If None, time points
            up to (and including) the last time point are included.
        **kwargs : dict
            Keyword arguments to be passed to transform_fun.

        Returns
        -------
        data_t : ndarray
            The transformed data.

        .. note::
            Applying transforms can be significantly faster if the
            SourceEstimate object was created using "(kernel, sens_data)", for
            the "data" parameter as the transform is applied in sensor space.
            Inverse methods, e.g., "apply_inverse_epochs", or "lcmv_epochs" do
            this automatically (if possible).
        """

        if idx is None:
            # use all time courses by default
            idx = slice(None, None)

        if fun_args is None:
            fun_args = tuple()

        if self._kernel is None and self._sens_data is None:
            # transform source space data directly
            data_t = transform_fun(self.data[idx, tmin_idx:tmax_idx],
                                   *fun_args, **kwargs)

            if isinstance(data_t, tuple):
                # use only first return value
                data_t = data_t[0]
        else:
            # apply transform in sensor space
            sens_data_t = transform_fun(self._sens_data[:, tmin_idx:tmax_idx],
                                        *fun_args, **kwargs)

            if isinstance(sens_data_t, tuple):
                # use only first return value
                sens_data_t = sens_data_t[0]

            # apply inverse
            data_shape = sens_data_t.shape
            if len(data_shape) > 2:
                # flatten the last dimensions
                sens_data_t = sens_data_t.reshape(data_shape[0],
                                                  np.prod(data_shape[1:]))

            data_t = np.dot(self._kernel[idx, :], sens_data_t)

            # restore original shape if necessary
            if len(data_shape) > 2:
                data_t = data_t.reshape(data_t.shape[0], *data_shape[1:])

        return data_t

    def transform(self, transform_fun, fun_args=None,
                  idx=None, tmin=None, tmax=None, copy=False, **kwargs):
        """Apply linear transform

        The transorm is applied to each source time course independently.


        Parameters
        ----------
        transform_fun : callable
            The transform to be applied. The first parameter of the function
            is the input data. The first return value is the transformed
            data, remaining outputs are ignored. The first dimension of the
            transformed data has to be the same as the first dimension of
            the input data.
        fun_args : tuple | None
            Additional parameters to be passed to transform_fun.
        idx : array | None
            Indicices of source time courses for which to compute transform.
            If None, all time courses are used.
        tmin : int | None
            First time point to include (ms). If None, self.tmin is used.
        tmax : int | None
            Last time point to include (ms). If None, self.tmax is used.
        copy: bool
            If True, returns copy instead of modifying inplace the input stc.
        **kwargs : dict
            Keyword arguments to be passed to transform_fun.

        Returns
        -------
        stcs : SourceEstimate | list
            The transformed stc or, in the case of transforms which yield
            N-dimensional output (where N > 2), a list of stcs.  For a list,
            copy must be True.


        .. note::
            Applying transforms can be significantly faster if the
            SourceEstimate object was created using "(kernel, sens_data)", for
            the "data" parameter as the transform is applied in sensor space.
            Inverse methods, e.g., "apply_inverse_epochs", or "lcmv_epochs" do
            this automatically (if possible).
        """

        if fun_args is None:
            fun_args = tuple()

        if idx is None:
            # use all time courses by default
            idx = slice(None, None)

        if tmin is None:
            tmin_idx = None
        else:
            tmin_idx = np.where(self.times >= tmin/1000.)[0][0]

        if tmax is None:
            tmax_idx = None
        else:
            tmax_idx = np.where(self.times <= tmax/1000.)[0][-1]

        if self._kernel is None and self._sens_data is None:
            # transform source space data directly
            data_t = transform_fun(self.data[idx, tmin_idx:tmax_idx],
                                   *fun_args, **kwargs)

            if isinstance(data_t, tuple):
                # use only first return value
                data_t = data_t[0]
        else:
            # apply transform in sensor space
            sens_data_t = transform_fun(self._sens_data[:, tmin_idx:tmax_idx],
                                        *fun_args, **kwargs)

            if isinstance(sens_data_t, tuple):
                # use only first return value
                sens_data_t = sens_data_t[0]

            # apply inverse
            data_shape = sens_data_t.shape
            if len(data_shape) > 2:
                # flatten the last dimensions
                sens_data_t = sens_data_t.reshape(data_shape[0],
                                                  np.prod(data_shape[1:]))

            data_t = np.dot(self._kernel[idx, :], sens_data_t)

            # restore original shape if necessary
            if len(data_shape) > 2:
                data_t = data_t.reshape(data_t.shape[0], *data_shape[1:])

        tmin = tmin if tmin is not None else self.tmin
        if len(data_t.shape) > 2:
            if copy:
                stcs = [SourceEstimate(data_t[:, :, a], self.vertno, tmin,
                                       self.tstep, self.subject)
                                       for a in range(data_t.shape[-1])]
                return stcs
            else:
                raise ValueError('copy must be True if transformed data has '
                                 'more than 2 dimensions')
        else:
            if copy:
                stc = SourceEstimate(data_t, self.vertno, tmin, self.tstep,
                                     self.subject)
                return stc
            else:
                self._data = data_t
                return self

    def as_data_frame(self, index=None, scale_time=1e3, copy=True):
        """Represent source estimates as Pandas DataFrame

        Export source estimates in tabular structure with vertices as columns
        and two additional info columns 'subject' and 'time'.
        This function is useful to visualize and analyse source time courses
        with external statistical software such as statsmodels or R.

        Parameters
        ----------
        index : tuple of str | None
            Column to be used as index for the data. Valid string options
            are 'subject' and 'time'. If None, both info
            columns will be included in the table as categorial data.
            If stc.subject is None, only time will be included.
        scale_time : float
            Scaling to be applied to time units.
        copy : bool
            If true, data will be copied. Else data may be modified in place.

        Returns
        -------
        df : instance of DataFrame
            Source estimates exported into tabular data structure.
        """
        pd = _check_pandas_installed()

        default_index = ['subject', 'time']
        if index is not None:<commit_msg>WIP: add transform method to SourceEstimate class<commit_after># Authors: <NAME> <<EMAIL>>
#          <NAME> <<EMAIL>>
#          <NAME> <<EMAIL>>
#
# License: BSD (3-clause)

import os
import copy
from math import ceil
import numpy as np
from scipy import linalg, sparse
from scipy.sparse import csr_matrix, coo_matrix
import warnings

from .filter import resample
from .parallel import parallel_func
from .surface import (read_surface, _get_ico_surface, read_morph_map,
                      _compute_nearest)
from .utils import (get_subjects_dir, _check_subject,
                    _check_pandas_index_arguments, _check_pandas_installed,
                    logger, verbose)
from .viz import plot_source_estimates
from . fixes import in1d


def _read_stc(filename):
    """ Aux Function
    """
    fid = open(filename, 'rb')

    stc = dict()

    fid.seek(0, 2)  # go to end of file
    file_length = fid.tell()
    fid.seek(0, 0)  # go to beginning of file

    # read tmin in ms
    stc['tmin'] = float(np.fromfile(fid, dtype=">f4", count=1))
    stc['tmin'] /= 1000.0

    # read sampling rate in ms
    stc['tstep'] = float(np.fromfile(fid, dtype=">f4", count=1))
    stc['tstep'] /= 1000.0

    # read number of vertices/sources
    vertices_n = int(np.fromfile(fid, dtype=">u4", count=1))

    # read the source vector
    stc['vertices'] = np.fromfile(fid, dtype=">u4", count=vertices_n)

    # read the number of timepts
    data_n = int(np.fromfile(fid, dtype=">u4", count=1))

    if (vertices_n and  # vertices_n can be 0 (empty stc)
            ((file_length / 4 - 4 - vertices_n) % (data_n * vertices_n)) != 0):
        raise ValueError('incorrect stc file size')

    # read the data matrix
    stc['data'] = np.fromfile(fid, dtype=">f4", count=vertices_n * data_n)
    stc['data'] = stc['data'].reshape([data_n, vertices_n]).T

    # close the file
    fid.close()
    return stc


def _write_stc(filename, tmin, tstep, vertices, data):
    """Write an STC file

    Parameters
    ----------
    filename : string
        The name of the STC file.
    tmin : float
        The first time point of the data in seconds.
    tstep : float
        Time between frames in seconds.
    vertices : array of integers
        Vertex indices (0 based).
    data : 2D array
        The data matrix (nvert * ntime).
    """
    fid = open(filename, 'wb')

    # write start time in ms
    fid.write(np.array(1000 * tmin, dtype='>f4').tostring())
    # write sampling rate in ms
    fid.write(np.array(1000 * tstep, dtype='>f4').tostring())
    # write number of vertices
    fid.write(np.array(vertices.shape[0], dtype='>u4').tostring())
    # write the vertex indices
    fid.write(np.array(vertices, dtype='>u4').tostring())

    # write the number of timepts
    fid.write(np.array(data.shape[1], dtype='>u4').tostring())
    #
    # write the data
    #
    fid.write(np.array(data.T, dtype='>f4').tostring())

    # close the file
    fid.close()


def _read_3(fid):
    """ Read 3 byte integer from file
    """
    data = np.fromfile(fid, dtype=np.uint8, count=3).astype(np.int32)

    out = np.left_shift(data[0], 16) + np.left_shift(data[1], 8) + data[2]

    return out


def _read_w(filename):
    """Read a w file and return as dict

    w files contain activations or source reconstructions for a single time
    point.

    Parameters
    ----------
    filename : string
        The name of the w file.

    Returns
    -------
    data: dict
        The w structure. It has the following keys:
           vertices       vertex indices (0 based)
           data           The data matrix (nvert long)
    """

    fid = open(filename, 'rb')

    # skip first 2 bytes
    fid.read(2)

    # read number of vertices/sources (3 byte integer)
    vertices_n = int(_read_3(fid))

    vertices = np.zeros((vertices_n), dtype=np.int32)
    data = np.zeros((vertices_n), dtype=np.float32)

    # read the vertices and data
    for i in range(vertices_n):
        vertices[i] = _read_3(fid)
        data[i] = np.fromfile(fid, dtype='>f4', count=1)

    w = dict()
    w['vertices'] = vertices
    w['data'] = data

    # close the file
    fid.close()
    return w


def _write_3(fid, val):
    """ Write 3 byte integer to file
    """

    f_bytes = np.zeros((3), dtype=np.uint8)

    f_bytes[0] = (val >> 16) & 255
    f_bytes[1] = (val >> 8) & 255
    f_bytes[2] = val & 255

    fid.write(f_bytes.tostring())


def _write_w(filename, vertices, data):
    """Read a w file

    w files contain activations or source reconstructions for a single time
    point.

    Parameters
    ----------
    filename: string
        The name of the w file.
    vertices: array of int
        Vertex indices (0 based).
    data: 1D array
        The data array (nvert).
    """

    assert(len(vertices) == len(data))

    fid = open(filename, 'wb')

    # write 2 zero bytes
    fid.write(np.zeros((2), dtype=np.uint8).tostring())

    # write number of vertices/sources (3 byte integer)
    vertices_n = len(vertices)
    _write_3(fid, vertices_n)

    # write the vertices and data
    for i in range(vertices_n):
        _write_3(fid, vertices[i])
        #XXX: without float() endianness is wrong, not sure why
        fid.write(np.array(float(data[i]), dtype='>f4').tostring())

    # close the file
    fid.close()


def read_source_estimate(fname, subject=None):
    """Read a soure estimate object

    Parameters
    ----------
    fname : str
        Path to (a) source-estimate file(s).
    subject : str | None
        Name of the subject the source estimate(s) is (are) from.
        It is good practice to set this attribute to avoid combining
        incompatible labels and SourceEstimates (e.g., ones from other
        subjects). Note that due to file specification limitations, the
        subject name isn't saved to or loaded from files written to disk.

    Returns
    -------
    stc : SourceEstimate | VolSourceEstimate
        The soure estimate object loaded from file.

    Notes
    -----
     - for volume source estimates, ``fname`` should provide the path to a
       single file named '*-vl.stc` or '*-vol.stc'
     - for surface source estimates, ``fname`` should either provide the
       path to the file corresponding to a single hemisphere ('*-lh.stc',
       '*-rh.stc') or only specify the asterisk part in these patterns. In any
       case, the function expects files for both hemisphere with names
       following this pattern.
     - for single time point .w files, ``fname`` should follow the same
       pattern as for surface estimates, except that files are named
       '*-lh.w' and '*-rh.w'.
    """
    fname_arg = fname

    # make sure corresponding file(s) can be found
    ftype = None
    if os.path.exists(fname):
        if fname.endswith('-vl.stc') or fname.endswith('-vol.stc') or \
                fname.endswith('-vl.w') or fname.endswith('-vol.w'):
            ftype = 'volume'
        elif fname.endswith('.stc'):
            ftype = 'surface'
            if fname.endswith(('-lh.stc', '-rh.stc')):
                fname = fname[:-7]
            else:
                err = ("Invalid .stc filename: %r; needs to end with "
                       "hemisphere tag ('...-lh.stc' or '...-rh.stc')"
                       % fname)
                raise IOError(err)
        elif fname.endswith('.w'):
            ftype = 'w'
            if fname.endswith(('-lh.w', '-rh.w')):
                fname = fname[:-5]
            else:
                err = ("Invalid .w filename: %r; needs to end with "
                       "hemisphere tag ('...-lh.w' or '...-rh.w')"
                       % fname)
                raise IOError(err)

    if ftype is not 'volume':
        stc_exist = map(os.path.exists, (fname + '-rh.stc', fname + '-lh.stc'))
        w_exist = map(os.path.exists, (fname + '-rh.w', fname + '-lh.w'))
        if all(stc_exist) and (ftype is not 'w'):
            ftype = 'surface'
        elif all(w_exist):
            ftype = 'w'
        elif any(stc_exist) or any(w_exist):
            raise IOError("Hemisphere missing for %r" % fname_arg)
        else:
            raise IOError("SourceEstimate File(s) not found for: %r"
                          % fname_arg)

    # read the files
    if ftype == 'volume':  # volume source space
        if fname.endswith('.stc'):
            kwargs = _read_stc(fname)
        elif fname.endswith('.w'):
            kwargs = _read_w(fname)
            kwargs['data'] = kwargs['data'][:, np.newaxis]
            kwargs['tmin'] = 0.0
            kwargs['tstep'] = 0.0
        else:
            raise IOError('Volume source estimate must end with .stc or .w')
    elif ftype == 'surface':  # stc file with surface source spaces
        lh = _read_stc(fname + '-lh.stc')
        rh = _read_stc(fname + '-rh.stc')
        assert lh['tmin'] == rh['tmin']
        assert lh['tstep'] == rh['tstep']
        kwargs = lh.copy()
        kwargs['data'] = np.r_[lh['data'], rh['data']]
        kwargs['vertices'] = [lh['vertices'], rh['vertices']]
    elif ftype == 'w':  # w file with surface source spaces
        lh = _read_w(fname + '-lh.w')
        rh = _read_w(fname + '-rh.w')
        kwargs = lh.copy()
        kwargs['data'] = np.atleast_2d(np.r_[lh['data'], rh['data']]).T
        kwargs['vertices'] = [lh['vertices'], rh['vertices']]
        # w files only have a single time point
        kwargs['tmin'] = 0.0
        kwargs['tstep'] = 1.0

    if ftype != 'volume':
        # Make sure the vertices are ordered
        vertices = kwargs['vertices']
        if any([np.any(np.diff(v.astype(int)) <= 0) for v in vertices]):
            sidx = [np.argsort(verts) for verts in vertices]
            vertices = [verts[idx] for verts, idx in zip(vertices, sidx)]
            data = kwargs['data'][np.r_[sidx[0], len(sidx[0]) + sidx[1]]]
            kwargs['vertices'] = vertices
            kwargs['data'] = data

    kwargs['subject'] = subject

    if ftype == 'volume':
        stc = VolSourceEstimate(**kwargs)
    else:
        stc = SourceEstimate(**kwargs)

    return stc


def _make_stc(data, vertices, tmin=None, tstep=None, subject=None):
    """Helper function to generate either a surface or volume source estimate
    """

    if isinstance(vertices, list) and len(vertices) == 2:
        # make a surface source estimate
        stc = SourceEstimate(data, vertices=vertices, tmin=tmin, tstep=tstep,
                             subject=subject)
    elif isinstance(vertices, np.ndarray) or isinstance(vertices, list)\
            and len(vertices) == 1:
        stc = VolSourceEstimate(data, vertices=vertices, tmin=tmin,
                                tstep=tstep, subject=subject)
    else:
        raise ValueError('vertices has to be either a list with one or two '
                         'arrays or an array')
    return stc


class _NotifyArray(np.ndarray):
    """Array class that executes a callback when it is modified
    """
    def __new__(cls, input_array, modify_callback=None):
        obj = np.asarray(input_array).view(cls)
        obj.modify_callback = modify_callback
        return obj

    def __array_finalize__(self, obj):
        if obj is None:
            # an empty constructor was used
            return

        # try to copy the callback
        self.modify_callback = getattr(obj, 'modify_callback', None)

    def _modified_(self):
        """Execute the callback if it is set"""
        if self.modify_callback is not None:
            self.modify_callback()

    def __getattribute__(self, name):
        # catch ndarray methods that modify the array inplace
        if name in ['fill', 'itemset', 'resize', 'sort']:
            self._modified_()

        return object.__getattribute__(self, name)

    def __setitem__(self, item, value):
        self._modified_()
        np.ndarray.__setitem__(self, item, value)

    def __array_wrap__(self, out_arr, context=None):
        # this method is called whenever a numpy ufunc (+, +=..) is called
        # the last entry in context is the array that receives the result
        if (context is not None and len(context[1]) == 3
                and context[1][2] is self):
            self._modified_()

        return np.ndarray.__array_wrap__(self, out_arr, context)


def _verify_source_estimate_compat(a, b):
    """Make sure two SourceEstimates are compatible for arith. operations"""
    compat = False
    if len(a.vertno) == len(b.vertno):
        if all([np.array_equal(av, vv) for av, vv in zip(a.vertno, b.vertno)]):
            compat = True
    if not compat:
        raise ValueError('Cannot combine SourceEstimates that do not have the '
                         'same vertices. Consider using stc.expand().')
    if a.subject != b.subject:
        raise ValueError('source estimates do not have the same subject '
                         'names, "%s" and "%s"' % (a.name, b.name))


class _BaseSourceEstimate(object):
    """Abstract base class for source estimates

    Parameters
    ----------
    data : array of shape (n_dipoles, n_times) | 2-tuple (kernel, sens_data)
        The data in source space. The data can either be a single array or
        a tuple with two arrays: "kernel" shape (n_vertices, n_sensors) and
        "sens_data" shape (n_sensors, n_times). In this case, the source
        space data corresponds to "numpy.dot(kernel, sens_data)".
    vertices : array | list of two arrays
        Vertex numbers corresponding to the data.
    tmin : scalar
        Time point of the first sample in data.
    tstep : scalar
        Time step between successive samples in data.
    subject : str | None
        The subject name. While not necessary, it is safer to set the
        subject parameter to avoid analysis errors.
    verbose : bool, str, int, or None
        If not None, override default verbose level (see mne.verbose).

    Attributes
    ----------
    subject : str | None
        The subject name.
    times : array of shape (n_times,)
        The time vector.
    vertices : array or list of arrays of shape (n_dipoles,)
        The indices of the dipoles in the different source spaces. Can
        be an array if there is only one source space (e.g., for volumes).
    data : array of shape (n_dipoles, n_times)
        The data in source space.
    shape : tuple
        The shape of the data. A tuple of int (n_dipoles, n_times).
    """
    @verbose
    def __init__(self, data, vertices=None, tmin=None, tstep=None,
                 subject=None, verbose=None):
        kernel, sens_data = None, None
        if isinstance(data, tuple):
            if len(data) != 2:
                raise ValueError('If data is a tuple it has to be length 2')
            kernel, sens_data = data
            data = None
            if kernel.shape[1] != sens_data.shape[0]:
                raise ValueError('kernel and sens_data have invalid '
                                 'dimensions')

        if isinstance(vertices, list):
            if not (len(vertices) == 2 or len(vertices) == 1) or \
                    not all([isinstance(v, np.ndarray) for v in vertices]):
                raise ValueError('Vertices, if a list, must contain one or '
                                 'two numpy arrays')

            if any([np.any(np.diff(v.astype(int)) <= 0) for v in vertices]):
                raise ValueError('Vertices must be ordered in increasing '
                                 'order.')

            n_src = sum([len(v) for v in vertices])

            if len(vertices) == 1:
                vertices = vertices[0]
        elif isinstance(vertices, np.ndarray):
            n_src = len(vertices)
        else:
            raise ValueError('Vertices must be a list or numpy array')

        # safeguard the user against doing something silly
        if data is not None and data.shape[0] != n_src:
            raise ValueError('Number of vertices (%i) and stc.shape[0] (%i) '
                             'must match' % (n_src, data.shape[0]))

        self._data = data
        self.tmin = tmin
        self.tstep = tstep
        self.vertno = vertices
        self.verbose = verbose
        self._kernel = kernel
        self._sens_data = sens_data
        self.times = None
        self._update_times()
        self.subject = _check_subject(None, subject, False)

    def _remove_kernel_sens_data_(self):
        """Remove kernel and sensor space data

        Note: self._data is also computed if it is None
        """
        if self._kernel is not None or self._sens_data is not None:
            # we can no longer use the kernel and sens_data
            logger.info('STC data modified: removing kernel and sensor data')
            if self._data is None:
                self._data = np.dot(self._kernel, self._sens_data)
            self._kernel = None
            self._sens_data = None

    def crop(self, tmin=None, tmax=None):
        """Restrict SourceEstimate to a time interval

        Parameters
        ----------
        tmin : float or None
            The first time point in seconds. If None the first present is used.
        tmax : float or None
            The last time point in seconds. If None the last present is used.
        """
        mask = np.ones(len(self.times), dtype=np.bool)
        if tmax is not None:
            mask = mask & (self.times <= tmax)
        if tmin is not None:
            mask = mask & (self.times >= tmin)
            self.tmin = tmin

        if self._kernel is not None and self._sens_data is not None:
            self._sens_data = self._sens_data[:, mask]
            self._data = None  # will be recomputed when data is accessed
        else:
            self._data = self._data[:, mask]

        self._update_times()
        return self  # return self for chaining methods

    @verbose
    def resample(self, sfreq, npad=100, window='boxcar', n_jobs=1,
                 verbose=None):
        """Resample data

        Parameters
        ----------
        sfreq : float
            New sample rate to use.
        npad : int
            Amount to pad the start and end of the data.
        window : string or tuple
            Window to use in resampling. See scipy.signal.resample.
        n_jobs : int
            Number of jobs to run in parallel.
        verbose : bool, str, int, or None
            If not None, override default verbose level (see mne.verbose).
            Defaults to self.verbose.

        Notes
        -----
        For some data, it may be more accurate to use npad=0 to reduce
        artifacts. This is dataset dependent -- check your data!

        Note that the sample rate of the original data is inferred from tstep.
        """
        # resampling in sensor instead of source space gives a somewhat
        # different result, so we don't allow it
        self._remove_kernel_sens_data_()

        o_sfreq = 1.0 / self.tstep
        self._data = resample(self._data, sfreq, o_sfreq, npad, n_jobs=n_jobs)

        # adjust indirectly affected variables
        self.tstep = 1.0 / sfreq
        self._update_times()

    @property
    def data(self):
        if self._data is None:
            # compute the solution the first time the data is accessed
            # return a "notify array", so we can later remove the kernel
            # and sensor data if the user modifies self._data
            self._data = _NotifyArray(np.dot(self._kernel, self._sens_data),
                modify_callback=self._remove_kernel_sens_data_)
        return self._data

    @property
    def shape(self):
        if self._data is not None:
            return self._data.shape
        return (self._kernel.shape[0], self._sens_data.shape[1])

    def _update_times(self):
        """Update the times attribute after changing tmin, tmax, or tstep"""
        self.times = self.tmin + (self.tstep * np.arange(self.shape[1]))

    def __add__(self, a):
        stc = copy.deepcopy(self)
        stc += a
        return stc

    def __iadd__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data += a.data
        else:
            self._data += a
        return self

    def mean(self):
        """Make a summary stc file with mean power between tmin and tmax.

        Returns
        -------
        stc : instance of SourceEstimate
            The modified stc (note: method operates inplace).
        """
        data = self.data
        tmax = self.tmin + self.tstep * data.shape[1]
        tmin = (self.tmin + tmax) / 2.
        tstep = tmax - self.tmin
        mean_stc = SourceEstimate(self.data.mean(axis=1)[:, np.newaxis],
                                  vertices=self.vertno, tmin=tmin,
                                  tstep=tstep, subject=self.subject)
        return mean_stc

    def __sub__(self, a):
        stc = copy.deepcopy(self)
        stc -= a
        return stc

    def __isub__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data -= a.data
        else:
            self._data -= a
        return self

    def __div__(self, a):
        stc = copy.deepcopy(self)
        stc /= a
        return stc

    def __idiv__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data /= a.data
        else:
            self._data /= a
        return self

    def __mul__(self, a):
        stc = copy.deepcopy(self)
        stc *= a
        return stc

    def __imul__(self, a):
        self._remove_kernel_sens_data_()
        if isinstance(a, _BaseSourceEstimate):
            _verify_source_estimate_compat(self, a)
            self._data *= a.data
        else:
            self._data *= a
        return self

    def __pow__(self, a):
        stc = copy.deepcopy(self)
        stc **= a
        return stc

    def __ipow__(self, a):
        self._remove_kernel_sens_data_()
        self._data **= a
        return self

    def __radd__(self, a):
        return self + a

    def __rsub__(self, a):
        return self - a

    def __rmul__(self, a):
        return self * a

    def __rdiv__(self, a):
        return self / a

    def __neg__(self):
        stc = copy.deepcopy(self)
        stc._remove_kernel_sens_data_()
        stc._data *= -1
        return stc

    def __pos__(self):
        return self

    def sqrt(self):
        """Return copy of SourceEstimate with sqrt(data)."""
        return self ** (0.5)

    def copy(self):
        """Return copy of SourceEstimate instance"""
        return copy.deepcopy(self)

    def bin(self, width, tstart=None, tstop=None, func=np.mean):
        """Returns a SourceEstimate object with data summarized over time bins

        Time bins of ``width`` seconds. This method is intended for
        visualization only. No filter is applied to the data before binning,
        making the method inappropriate as a tool for downsampling data.

        Parameters
        ----------
        width : scalar
            Width of the individual bins in seconds.
        func : callable
            Function that is applied to summarize the data. Needs to accept a
            numpy.array as first input and an ``axis`` keyword argument.
        tstart : scalar | None
            Time point where the first bin starts. The default is the first
            time point of the stc.
        tstop : scalar | None
            Last possible time point contained in a bin (if the last bin would
            be shorter than width it is dropped). The default is the last time
            point of the stc.

        Returns
        -------
        stc : instance of SourceEstimate
            The binned SourceEstimate.
        """
        if tstart is None:
            tstart = self.tmin
        if tstop is None:
            tstop = self.times[-1]

        times = np.arange(tstart, tstop + self.tstep, width)
        nv, _ = self.shape
        nt = len(times) - 1
        data = np.empty((nv, nt), dtype=self.data.dtype)
        for i in xrange(nt):
            idx = (self.times >= times[i]) & (self.times < times[i + 1])
            data[:, i] = func(self.data[:, idx], axis=1)

        tmin = times[0] + width / 2.
        stc = _make_stc(data, vertices=self.vertno,
                        tmin=tmin, tstep=width, subject=self.subject)
        return stc

    def transform_data(self, transform_fun, fun_args=None,
                       idx=None, tmin_idx=None, tmax_idx=None, **kwargs):
        """Get data after a linear (time) transform has been applied

        The transorm is applied to each source time course independently.


        Parameters
        ----------
        transform_fun : callable
            The transform to be applied. The first parameter of the function
            is the input data. The first return value is the transformed
            data, remaining outputs are ignored. The first dimension of the
            transformed data has to be the same as the first dimension of the
            input data.
        fun_args : tuple | None
            Additional parameters to be passed to transform_fun.
        idx : array | None
            Indicices of source time courses for which to compute transform.
            If None, all time courses are used.
        tmin_idx : int | None
            Index of first time point to include. If None, the index of the
            first time point is used.
        tmax_idx : int | None
            Index of the first time point not to include. If None, time points
            up to (and including) the last time point are included.
        **kwargs : dict
            Keyword arguments to be passed to transform_fun.

        Returns
        -------
        data_t : ndarray
            The transformed data.

        .. note::
            Applying transforms can be significantly faster if the
            SourceEstimate object was created using "(kernel, sens_data)", for
            the "data" parameter as the transform is applied in sensor space.
            Inverse methods, e.g., "apply_inverse_epochs", or "lcmv_epochs" do
            this automatically (if possible).
        """

        if idx is None:
            # use all time courses by default
            idx = slice(None, None)

        if fun_args is None:
            fun_args = tuple()

        if self._kernel is None and self._sens_data is None:
            # transform source space data directly
            data_t = transform_fun(self.data[idx, tmin_idx:tmax_idx],
                                   *fun_args, **kwargs)

            if isinstance(data_t, tuple):
                # use only first return value
                data_t = data_t[0]
        else:
            # apply transform in sensor space
            sens_data_t = transform_fun(self._sens_data[:, tmin_idx:tmax_idx],
                                        *fun_args, **kwargs)

            if isinstance(sens_data_t, tuple):
                # use only first return value
                sens_data_t = sens_data_t[0]

            # apply inverse
            data_shape = sens_data_t.shape
            if len(data_shape) > 2:
                # flatten the last dimensions
                sens_data_t = sens_data_t.reshape(data_shape[0],
                                                  np.prod(data_shape[1:]))

            data_t = np.dot(self._kernel[idx, :], sens_data_t)

            # restore original shape if necessary
            if len(data_shape) > 2:
                data_t = data_t.reshape(data_t.shape[0], *data_shape[1:])

        return data_t

    def transform(self, func, func_args=None,
                  idx=None, tmin=None, tmax=None, copy=False, **kwargs):
        """Apply linear transform

        The transorm is applied to each source time course independently.

        Parameters
        ----------
        transform_fun : callable
            The transform to be applied. The first parameter of the function
            is the input data. The first return value is the transformed
            data, remaining outputs are ignored. The first dimension of the
            transformed data has to be the same as the first dimension of
            the input data.
        fun_args : tuple | None
            Additional parameters to be passed to transform_fun.
        idx : array | None
            Indicices of source time courses for which to compute transform.
            If None, all time courses are used.
        tmin : int | None
            First time point to include (ms). If None, self.tmin is used.
        tmax : int | None
            Last time point to include (ms). If None, self.tmax is used.
        copy: bool
            If True, returns copy instead of modifying inplace the input stc.
        **kwargs : dict
            Keyword arguments to be passed to transform_fun.

        Returns
        -------
        stcs : instance of SourceEstimate | list
            The transformed stc or, in the case of transforms which yield
            N-dimensional output (where N > 2), a list of stcs.  For a list,
            copy must be True.

        Notes
        -----
            Applying transforms can be significantly faster if the
            SourceEstimate object was created using "(kernel, sens_data)", for
            the "data" parameter as the transform is applied in sensor space.
            Inverse methods, e.g., "apply_inverse_epochs", or "lcmv_epochs" do
            this automatically (if possible).
        """

        if tmin is None:
            tmin_idx = None
        else:
            tmin_idx = np.where(self.times >= tmin / 1000.)[0][0]

        if tmax is None:
            tmax_idx = None
        else:
            tmax_idx = np.where(self.times <= tmax / 1000.)[0][-1]

        data_t = self.transform_data(func, fun_args=func_args, idx=idx,
                                     tmin_idx=tmin_idx, tmax_idx=tmax_idx,
                                     **kwargs)

        if idx is not None:
            idx_lh = idx[idx < len(self.lh_vertno)]
            idx_rh = idx[idx >= len(self.lh_vertno)] - len(self.lh_vertno)
            verts_lh = self.lh_vertno[idx_lh]
            verts_rh = self.rh_vertno[idx_rh]
        else:
            verts_lh = self.lh_vertno
            verts_rh = self.rh_vertno
        verts = [verts_lh, verts_rh]

        tmin = tmin if tmin is not None else self.tmin

        if len(data_t.shape) > 2:
            if copy:
                stcs = [SourceEstimate(data_t[:, :, a], verts, tmin,
                                       self.tstep, self.subject)
                                       for a in range(data_t.shape[-1])]
            else:
                raise ValueError('copy must be True if transformed data has '
                                 'more than 2 dimensions')
        else:
            if copy:
                stcs = SourceEstimate(data_t, verts, tmin, self.tstep,
                                      self.subject)
            else:
                self._data, self.vertno = data_t, verts
                stcs = self

        return stcs

    def as_data_frame(self, index=None, scale_time=1e3, copy=True):
        """Represent source estimates as Pandas DataFrame

        Export source estimates in tabular structure with vertices as columns
        and two additional info columns 'subject' and 'time'.
        This function is useful to visualize and analyse source time courses
        with external statistical software such as statsmodels or R.

        Parameters
        ----------
        index : tuple of str | None
            Column to be used as index for the data. Valid string options
            are 'subject' and 'time'. If None, both info
            columns will be included in the table as categorial data.
            If stc.subject is None, only time will be included.
        scale_time : float
            Scaling to be applied to time units.
        copy : bool
            If true, data will be copied. Else data may be modified in place.

        Returns
        -------
        df : instance of DataFrame
            Source estimates exported into tabular data structure.
        """
        pd = _check_pandas_installed()

        default_index = ['subject', 'time']
        if index is not None:
{'fraction_non_alphanumeric': 0.09734311278836204, 'fraction_numerical': 0.005721589871568568, 'mean_word_length': 1.8681098066599746, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 74, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7622310', 'n_tokens_mistral': 19310, 'n_tokens_neox': 18626, 'n_words': 6358}