EventArray

class BaseEventArray(ABC):
    """Base class for EventArray and BinnedEventArray.

    NOTE: This class can't really be instantiated, almost like a pseudo
    abstract class. In particular, during initialization it might fail
    because it checks the n_series of its derived classes to validate
    input to series_ids and series_labels. If NoneTypes are used, then you
    may actually succeed in creating an instance of this class, but it
    will be pretty useless.

    This docstring only applies to this base class, as subclasses may have
    different behavior. Therefore read the particular subclass' docstring
    for the most accurate information.

    Parameters
    ----------
    fs: float, optional
        Sampling rate / frequency (Hz).
    series_ids : list of int, optional
        Unit IDs
    series_labels : list of str, optional
        Labels corresponding to series. Default casts series_ids to str.
    series_tags : optional
        Tags correponding to series.
        NOTE: Currently we do not do any input validation so these can
        be any type. We also don't use these for anything yet.
    label : str, optional
        Information pertaining to the source of the event series.
        Default is None.
    empty : bool, optional
        Whether to create an empty class instance (no data).
        Default is False
    abscissa : optional
        Object for the abscissa (x-axis) coordinate
    ordinate : optional
        Object for the ordinate (y-axis) coordinate
    kwargs : optional
        Other keyword arguments. NOTE: Currently we do not do anything
        with these.

    Attributes
    ----------
    n_series : int
        Number of series in event series.
    issempty : bool
        Whether the class instance is empty (no data)
    series_ids : list of int
        Unit IDs
    series_labels : list of str
        Labels corresponding to series. Default casts series_ids to str.
    series_tags :
        Tags corresponding to series.
        NOTE: Currently we do not do any input validation so these can
        be any type. We also don't use these for anything yet.
    n_intervals : int
        The number of underlying intervals.
    support : nelpy.IntervalArray
        The support of the EventArray.
    fs: float
        Sampling frequency (Hz).
    label : str or None
        Information pertaining to the source of the event series.
    """

    __aliases__ = {}

    __attributes__ = ["_fs", "_series_ids", "_series_labels", "_series_tags", "_label"]

    def __init__(
        self,
        *,
        fs=None,
        series_ids=None,
        series_labels=None,
        series_tags=None,
        label=None,
        empty=False,
        abscissa=None,
        ordinate=None,
        **kwargs,
    ):
        self.__version__ = version.__version__
        self.type_name = self.__class__.__name__
        if abscissa is None:
            abscissa = core.Abscissa()  # TODO: integrate into constructor?
        if ordinate is None:
            ordinate = core.Ordinate()  # TODO: integrate into constructor?
        self._abscissa = abscissa
        self._ordinate = ordinate

        series_label = kwargs.pop("series_label", None)
        if series_label is None:
            series_label = "series"
        self._series_label = series_label

        # if an empty object is requested, return it:
        if empty:
            for attr in self.__attributes__:
                exec("self." + attr + " = None")
            self._abscissa.support = type(self._abscissa.support)(empty=True)
            self.loc = _accessors.ItemGetterLoc(self)
            self.iloc = _accessors.ItemGetterIloc(self)
            return

        # set initial fs to None
        self._fs = None
        # then attempt to update the fs; this does input validation:
        self.fs = fs

        # WARNING! we need to ensure that self.n_series can work BEFORE
        # we can set self.series_ids or self.series_labels, since those
        # setters check that the lengths of the inputs are consistent
        # with self.n_series.

        # inherit series IDs if available, otherwise initialize to default
        if series_ids is None:
            series_ids = list(range(1, self.n_series + 1))

        series_ids = np.array(series_ids, ndmin=1)  # standardize series_ids

        # if series_labels is empty, default to series_ids
        if series_labels is None:
            series_labels = series_ids

        series_labels = np.array(series_labels, ndmin=1)  # standardize

        self.series_ids = series_ids
        self.series_labels = series_labels
        self._series_tags = series_tags  # no input validation yet
        self.label = label

        self.loc = _accessors.ItemGetterLoc(self)
        self.iloc = _accessors.ItemGetterIloc(self)

    def __renew__(self):
        """Re-attach slicers and indexers."""
        self.loc = _accessors.ItemGetterLoc(self)
        self.iloc = _accessors.ItemGetterIloc(self)

    def __repr__(self):
        address_str = " at " + str(hex(id(self)))
        return "<BaseEventArray" + address_str + ">"

    @abstractmethod
    @keyword_equivalence(this_or_that={"n_intervals": "n_epochs"})
    def partition(self, ds=None, n_intervals=None):
        """
        Returns a BaseEventArray whose support has been partitioned.

        Parameters
        ----------
        ds : float, optional
            Maximum duration (in seconds), for each interval.
        n_intervals : int, optional
            Number of intervals. If ds is None and n_intervals is None, then
            default is to use n_intervals = 100

        Returns
        -------
        out : BaseEventArray
            BaseEventArray that has been partitioned.

        Notes
        -----
        Irrespective of whether 'ds' or 'n_intervals' are used, the exact
        underlying support is propagated, and the first and last points
        of the supports are always included, even if this would cause
        n_points or ds to be violated.
        """
        return

    @abstractmethod
    def isempty(self):
        """(bool) Empty BaseEventArray."""
        return

    @abstractmethod
    def n_series(self):
        """(int) The number of series."""
        return

    @property
    def n_intervals(self):
        if self.isempty:
            return 0
        """(int) The number of underlying intervals."""
        return self._abscissa.support.n_intervals

    @property
    def series_ids(self):
        """Unit IDs contained in the BaseEventArray."""
        return self._series_ids

    @series_ids.setter
    def series_ids(self, val):
        if len(val) != self.n_series:
            raise TypeError("series_ids must be of length n_series")
        elif len(set(val)) < len(val):
            raise TypeError("duplicate series_ids are not allowed")
        else:
            try:
                # cast to int:
                series_ids = [int(id) for id in val]
            except TypeError:
                raise TypeError("series_ids must be int-like")
        self._series_ids = series_ids

    @property
    def series_labels(self):
        """Labels corresponding to series contained in the BaseEventArray."""
        if self._series_labels is None:
            logging.warning("series labels have not yet been specified")
            return self.series_ids
        return self._series_labels

    @series_labels.setter
    def series_labels(self, val):
        if len(val) != self.n_series:
            raise TypeError("labels must be of length n_series")
        else:
            try:
                # cast to str:
                labels = [str(label) for label in val]
            except TypeError:
                raise TypeError("labels must be string-like")
        self._series_labels = labels

    @property
    def series_tags(self):
        """Tags corresponding to series contained in the BaseEventArray"""
        if self._series_tags is None:
            logging.warning("series tags have not yet been specified")
        return self._series_tags

    @property
    def support(self):
        """(nelpy.IntervalArray) The support of the EventArray."""
        return self._abscissa.support

    @support.setter
    def support(self, val):
        """(nelpy.IntervalArray) The support of the EventArray."""
        # modify support
        if isinstance(val, type(self._abscissa.support)):
            self._abscissa.support = val
        elif isinstance(val, (tuple, list)):
            prev_domain = self._abscissa.domain
            self._abscissa.support = type(self._abscissa.support)([val[0], val[1]])
            self._abscissa.domain = prev_domain
        else:
            raise TypeError(
                "support must be of type {}".format(str(type(self._abscissa.support)))
            )
        # restrict data to new support
        self._restrict_to_interval(self._abscissa.support)

    @property
    def domain(self):
        """(nelpy.IntervalArray) The domain of the EventArray."""
        return self._abscissa.domain

    @domain.setter
    def domain(self, val):
        """(nelpy.IntervalArray) The domain of the EventArray."""
        # modify domain
        if isinstance(val, type(self._abscissa.support)):
            self._abscissa.domain = val
        elif isinstance(val, (tuple, list)):
            self._abscissa.domain = type(self._abscissa.support)([val[0], val[1]])
        else:
            raise TypeError(
                "support must be of type {}".format(str(type(self._abscissa.support)))
            )
        # restrict data to new support
        self._restrict_to_interval(self._abscissa.support)

    @property
    def fs(self):
        """(float) Sampling rate."""
        return self._fs

    @fs.setter
    def fs(self, val):
        """(float) Sampling rate."""
        if self._fs == val:
            return
        try:
            if val <= 0:
                raise ValueError("sampling rate must be positive")
        except TypeError:
            raise TypeError("sampling rate must be a scalar")
        self._fs = val

    @property
    def label(self):
        """Label pertaining to the source of the event series."""
        if self._label is None:
            logging.warning("label has not yet been specified")
        return self._label

    @label.setter
    def label(self, val):
        if val is not None:
            try:  # cast to str:
                label = str(val)
            except TypeError:
                raise TypeError("cannot convert label to string")
        else:
            label = val
        self._label = label

    def _series_subset(self, series_list):
        """Return a BaseEventArray restricted to a subset of series.

        Parameters
        ----------
        series_list : array-like
            Array or list of series_ids.
        """
        return self.loc[:, series_list]

    def __setattr__(self, name, value):
        # https://stackoverflow.com/questions/4017572/how-can-i-make-an-alias-to-a-non-function-member-attribute-in-a-python-class
        name = self.__aliases__.get(name, name)
        object.__setattr__(self, name, value)

    def __getattr__(self, name):
        # https://stackoverflow.com/questions/4017572/how-can-i-make-an-alias-to-a-non-function-member-attribute-in-a-python-class
        if name == "aliases":
            raise AttributeError  # http://nedbatchelder.com/blog/201010/surprising_getattr_recursion.html
        name = self.__aliases__.get(name, name)
        # return getattr(self, name) #Causes infinite recursion on non-existent attribute
        return object.__getattribute__(self, name)

class BinnedEventArray(BaseEventArray):
    """BinnedEventArray.

    Parameters
    ----------
    eventarray : nelpy.EventArray or nelpy.RegularlySampledAnalogSignalArray
        Input data.
    ds : float
        The bin width, in seconds.
        Default is 0.0625 (62.5 ms)
    empty : bool, optional
        Whether an empty BinnedEventArray should be constructed (no data).
    fs : float, optional
        Sampling rate in Hz. If fs is passed as a parameter, then data
        is assumed to be in sample numbers instead of actual data.
    kwargs : optional
        Additional keyword arguments to forward along to the BaseEventArray
        constructor.

    Attributes
    ----------
    Note : Read the docstring for the BaseEventArray superclass for additional
    attributes that are defined there.
    isempty : bool
        Whether the BinnedEventArray is empty (no data).
    n_series : int
        The number of series.
    bin_centers : np.ndarray
        The bin centers, in seconds.
    event_centers : np.ndarray
        The centers of each event, in seconds.
    data : np.array, with shape (n_series, n_bins)
        Event counts in all bins.
    bins : np.ndarray
        The bin edges, in seconds.
    binned_support : np.ndarray, with shape (n_intervals, 2)
        The binned support of the BinnedEventArray (in
        bin IDs).
    lengths : np.ndarray
        Lengths of contiguous segments, in number of bins.
    eventarray : nelpy.EventArray
        The original eventarray associated with the binned data.
    n_bins : int
        The number of bins.
    ds : float
        Bin width, in seconds.
    n_active : int
        The number of active series. A series is considered active if
        it fired at least one event.
    n_active_per_bin : np.ndarray, with shape (n_bins, )
        Number of active series per data bin.
    n_events : np.ndarray
        The number of events in each series.
    support : nelpy.IntervalArray
        The support of the BinnedEventArray.
    """

    __attributes__ = [
        "_ds",
        "_bins",
        "_data",
        "_bin_centers",
        "_binned_support",
        "_eventarray",
    ]
    __attributes__.extend(BaseEventArray.__attributes__)

    def __init__(self, eventarray=None, *, ds=None, empty=False, **kwargs):
        super().__init__(empty=True)

        # if an empty object is requested, return it:
        if empty:
            # super().__init__(empty=True)
            for attr in self.__attributes__:
                exec("self." + attr + " = None")
            self._abscissa.support = type(self._abscissa.support)(empty=True)
            self._event_centers = None
            return

        # handle casting other nelpy objects to BinnedEventArray:
        if isinstance(eventarray, core.RegularlySampledAnalogSignalArray):
            if eventarray.isempty:
                for attr in self.__attributes__:
                    exec("self." + attr + " = None")
                self._abscissa.support = type(eventarray._abscissa.support)(empty=True)
                self._event_centers = None
                return
            eventarray = eventarray.copy()  # Note: this is a deep copy
            n_empty_epochs = np.sum(eventarray.support.lengths == 0)
            if n_empty_epochs > 0:
                logging.warning(
                    "Detected {} empty epochs. Removing these in the cast object".format(
                        n_empty_epochs
                    )
                )
                eventarray.support = eventarray.support._drop_empty_intervals()
            if not eventarray.support.ismerged:
                logging.warning(
                    "Detected overlapping epochs. Merging these in the cast object"
                )
                eventarray.support = eventarray.support.merge()

            self._eventarray = None
            self._ds = 1 / eventarray.fs
            self._series_labels = eventarray._series_labels
            self._bin_centers = eventarray.abscissa_vals
            tmp = np.insert(np.cumsum(eventarray.lengths), 0, 0)
            self._binned_support = np.array((tmp[:-1], tmp[1:] - 1)).T
            self._abscissa.support = eventarray.support
            try:
                self._series_ids = (
                    np.array(eventarray.series_labels).astype(int)
                ).tolist()
            except (ValueError, TypeError):
                self._series_ids = (np.arange(eventarray.n_signals) + 1).tolist()
            self._data = eventarray._ydata_rowsig

            bins = []
            for starti, stopi in self._binned_support:
                bins_edges_in_interval = (
                    self._bin_centers[starti : stopi + 1] - self._ds / 2
                ).tolist()
                bins_edges_in_interval.append(self._bin_centers[stopi] + self._ds / 2)
                bins.extend(bins_edges_in_interval)
            self._bins = np.array(bins)
            return

        if type(eventarray).__name__ == "BinnedSpikeTrainArray":
            # old-style nelpy BinnedSpikeTrainArray object?
            try:
                self._eventarray = eventarray._spiketrainarray
                self._ds = eventarray.ds
                self._series_labels = eventarray.unit_labels
                self._bin_centers = eventarray.bin_centers
                self._binned_support = eventarray.binned_support
                try:
                    self._abscissa.support = core.EpochArray(eventarray.support.data)
                except AttributeError:
                    self._abscissa.support = core.EpochArray(eventarray.support.time)
                self._series_ids = eventarray.unit_ids
                self._data = eventarray.data
                return
            except Exception:
                pass

        if not isinstance(eventarray, EventArray):
            raise TypeError("eventarray must be a nelpy.EventArray object.")

        self._ds = None
        self._bin_centers = np.array([])
        self._event_centers = None

        logging.disable(logging.CRITICAL)
        kwargs = {
            "fs": eventarray.fs,
            "series_ids": eventarray.series_ids,
            "series_labels": eventarray.series_labels,
            "series_tags": eventarray.series_tags,
            "label": eventarray.label,
        }
        logging.disable(0)

        # initialize super so that self.fs is set:
        self._data = np.zeros((eventarray.n_series, 0))
        # the above is necessary so that super() can determine
        # self.n_series when initializing. self.data will
        # be updated later in __init__ to reflect subsequent changes
        super().__init__(**kwargs)

        if ds is None:
            logging.warning("no bin size was given, assuming 62.5 ms")
            ds = 0.0625

        self._eventarray = eventarray  # TODO: remove this if we don't need it, or decide that it's too wasteful
        self._abscissa = copy.deepcopy(eventarray._abscissa)
        self.ds = ds

        self._bin_events(eventarray=eventarray, intervalArray=eventarray.support, ds=ds)

    def __mul__(self, other):
        """Overloaded * operator"""

        if isinstance(other, numbers.Number):
            neweva = self.copy()
            neweva._data = self.data * other
            return neweva
        elif isinstance(other, np.ndarray):
            neweva = self.copy()
            neweva._data = (self.data.T * other).T
            return neweva
        else:
            raise TypeError(
                "unsupported operand type(s) for *: '{}' and '{}'".format(
                    str(type(self)), str(type(other))
                )
            )

    def __rmul__(self, other):
        """Overloaded * operator"""
        return self.__mul__(other)

    def __sub__(self, other):
        """Overloaded - operator"""
        if isinstance(other, numbers.Number):
            neweva = self.copy()
            neweva._data = self.data - other
            return neweva
        elif isinstance(other, np.ndarray):
            neweva = self.copy()
            neweva._data = (self.data.T - other).T
            return neweva
        else:
            raise TypeError(
                "unsupported operand type(s) for -: '{}' and '{}'".format(
                    str(type(self)), str(type(other))
                )
            )

    def __add__(self, other):
        """Overloaded + operator"""

        if isinstance(other, numbers.Number):
            neweva = self.copy()
            neweva._data = self.data + other
            return neweva
        elif isinstance(other, np.ndarray):
            neweva = self.copy()
            neweva._data = (self.data.T + other).T
            return neweva
        elif isinstance(other, type(self)):
            # TODO: additional checks need to be done, e.g., same series ids...
            assert self.n_series == other.n_series
            support = self._abscissa.support + other.support

            newdata = []
            for series in range(self.n_series):
                newdata.append(np.append(self.data[series], other.data[series]))

            fs = self.fs
            if self.fs != other.fs:
                fs = None
            return type(self)(newdata, support=support, fs=fs)
        else:
            raise TypeError(
                "unsupported operand type(s) for +: '{}' and '{}'".format(
                    str(type(self)), str(type(other))
                )
            )

    def __truediv__(self, other):
        """Overloaded / operator"""

        if isinstance(other, numbers.Number):
            neweva = self.copy()
            neweva._data = self.data / other
            return neweva
        elif isinstance(other, np.ndarray):
            neweva = self.copy()
            neweva._data = (self.data.T / other).T
            return neweva
        else:
            raise TypeError(
                "unsupported operand type(s) for /: '{}' and '{}'".format(
                    str(type(self)), str(type(other))
                )
            )

    def median(self, *, axis=1):
        """Returns the median of each series in BinnedEventArray."""
        try:
            medians = np.nanmedian(self.data, axis=axis).squeeze()
            if medians.size == 1:
                return medians.item()
            return medians
        except IndexError:
            raise IndexError("Empty BinnedEventArray; cannot calculate median.")

    def mean(self, *, axis=1):
        """Returns the mean of each series in BinnedEventArray."""
        try:
            means = np.nanmean(self.data, axis=axis).squeeze()
            if means.size == 1:
                return means.item()
            return means
        except IndexError:
            raise IndexError("Empty BinnedEventArray; cannot calculate mean.")

    def std(self, *, axis=1):
        """Returns the standard deviation of each series in BinnedEventArray."""
        try:
            stds = np.nanstd(self.data, axis=axis).squeeze()
            if stds.size == 1:
                return stds.item()
            return stds
        except IndexError:
            raise IndexError(
                "Empty BinnedEventArray; cannot calculate standard deviation"
            )

    def center(self, inplace=False):
        """Center data (zero mean)."""
        if inplace:
            out = self
        else:
            out = self.copy()
        out._data = (out._data.T - out.mean()).T
        return out

    def normalize(self, inplace=False):
        """Normalize data (unit standard deviation)."""
        if inplace:
            out = self
        else:
            out = self.copy()
        std = out.std()
        std[std == 0] = 1
        out._data = (out._data.T / std).T
        return out

    def standardize(self, inplace=False):
        """Standardize data (zero mean and unit std deviation)."""
        if inplace:
            out = self
        else:
            out = self.copy()
        out._data = (out._data.T - out.mean()).T
        std = out.std()
        std[std == 0] = 1
        out._data = (out._data.T / std).T

        return out

    @keyword_equivalence(this_or_that={"n_intervals": "n_epochs"})
    def partition(self, ds=None, n_intervals=None):
        """Returns a BaseEventArray whose support has been partitioned.

        Parameters
        ----------
        ds : float, optional
            Maximum duration (in seconds), for each interval.
        n_points : int, optional
            Number of intervals. If ds is None and n_intervals is None, then
            default is to use n_intervals = 100

        Returns
        -------
        out : BaseEventArray
            BaseEventArray that has been partitioned.

        Notes
        -----
        Irrespective of whether 'ds' or 'n_intervals' are used, the exact
        of the supports are always included, even if this would cause
        n_points or ds to be violated.
        """

        partitioned = type(self)(
            core.RegularlySampledAnalogSignalArray(self).partition(
                ds=ds, n_intervals=n_intervals
            )
        )
        # partitioned.loc = ItemGetter_loc(partitioned)
        # partitioned.iloc = ItemGetter_iloc(partitioned)
        return partitioned

        # raise NotImplementedError('workaround: cast to AnalogSignalArray, partition, and cast back to BinnedEventArray')

    def _copy_without_data(self):
        """Returns a copy of the BinnedEventArray, without data.
        Note: the support is left unchanged, but the binned_support is removed.
        """
        out = copy.copy(self)  # shallow copy
        out._bin_centers = None
        out._binned_support = None
        out._bins = None
        out._data = np.zeros((self.n_series, 0))
        out._eventarray = out._eventarray._copy_without_data()
        out = copy.deepcopy(
            out
        )  # just to be on the safe side, but at least now we are not copying the data!
        out.__renew__()
        return out

    def copy(self):
        """Returns a copy of the BinnedEventArray."""
        newcopy = copy.deepcopy(self)
        newcopy.__renew__()
        return newcopy

    def __repr__(self):
        address_str = " at " + str(hex(id(self)))
        if self.isempty:
            return "<empty " + self.type_name + address_str + ">"
        # Summarize labels if too many
        max_display = 5
        labels = self._series_labels
        n = len(labels)
        if n <= max_display:
            label_str = str(labels)
        else:
            label_str = f"[{', '.join(map(str, labels[:3]))}, ..., {', '.join(map(str, labels[-2:]))}]"
        ustr = f" {self.n_series} {label_str}"
        if self._abscissa.support.n_intervals > 1:
            epstr = f" ({self._abscissa.support.n_intervals} segments) in"
        else:
            epstr = " in"
        if self.n_bins == 1:
            bstr = f" {self.n_bins} bin of width {utils.PrettyDuration(self.ds)}"
            dstr = ""
        else:
            bstr = f" {self.n_bins} bins of width {utils.PrettyDuration(self.ds)}"
            dstr = f" for a total of {utils.PrettyDuration(self.n_bins * self.ds)}"
        return f"<{self.type_name}{address_str}:{ustr}{epstr}{bstr}>{dstr}"

    def __iter__(self):
        """BinnedEventArray iterator initialization."""
        # initialize the internal index to zero when used as iterator
        self._index = 0
        return self

    def __next__(self):
        """BinnedEventArray iterator advancer."""
        index = self._index

        if index > self._abscissa.support.n_intervals - 1:
            raise StopIteration

        # TODO: return self.loc[index], and make sure that __getitem__ is updated
        logging.disable(logging.CRITICAL)
        support = self._abscissa.support[index]
        bsupport = self.binned_support[[index], :]

        binnedeventarray = type(self)(empty=True)
        exclude = ["_bins", "_data", "_support", "_bin_centers", "_binned_support"]
        attrs = (x for x in self.__attributes__ if x not in exclude)
        for attr in attrs:
            exec("binnedeventarray." + attr + " = self." + attr)
        binindices = np.insert(0, 1, np.cumsum(self.lengths + 1))  # indices of bins
        binstart = binindices[index]
        binstop = binindices[index + 1]
        binnedeventarray._bins = self._bins[binstart:binstop]
        binnedeventarray._data = self._data[:, bsupport[0][0] : bsupport[0][1] + 1]
        binnedeventarray._abscissa.support = support
        binnedeventarray._bin_centers = self._bin_centers[
            bsupport[0][0] : bsupport[0][1] + 1
        ]
        binnedeventarray._binned_support = bsupport - bsupport[0, 0]
        logging.disable(0)
        self._index += 1
        binnedeventarray.__renew__()
        return binnedeventarray

    def empty(self, *, inplace=False):
        """Remove data (but not metadata) from BinnedEventArray.

        Attributes 'data', and 'support' 'binned_support' are all emptied.

        Note: n_series, series_ids, etc. are all preserved.
        """
        if not inplace:
            out = self._copy_without_data()
            out._abscissa.support = type(self._abscissa.support)(empty=True)
            return out
        out = self
        out._data = np.zeros((self.n_series, 0))
        out._abscissa.support = type(self._abscissa.support)(empty=True)
        out._binned_support = None
        out._bin_centers = None
        out._bins = None
        out._eventarray.empty(inplace=True)
        out.__renew__()
        return out

    def __getitem__(self, idx):
        """BinnedEventArray index access.

        By default, this method is bound to .loc
        """
        return self.loc[idx]

    def _restrict(self, intervalslice, seriesslice):
        # This function should be called only by an itemgetter
        # because it mutates data.
        # The itemgetter is responsible for creating copies
        # of objects

        self._restrict_to_series_subset(seriesslice)
        self._eventarray._restrict_to_series_subset(seriesslice)

        self._restrict_to_interval(intervalslice)
        self._eventarray._restrict_to_interval(intervalslice)
        return self

    def _restrict_to_series_subset(self, idx):
        # Warning: This function can mutate data

        if isinstance(idx, core.IntervalArray):
            raise IndexError(
                "Slicing is [intervals, signal]; perhaps you have the order reversed?"
            )

        # TODO: update tags
        try:
            self._data = np.atleast_2d(self.data[idx, :])
            self._series_ids = list(np.atleast_1d(np.atleast_1d(self._series_ids)[idx]))
            self._series_labels = list(
                np.atleast_1d(np.atleast_1d(self._series_labels)[idx])
            )
        except IndexError:
            raise IndexError(
                "One of more indices were out of bounds for n_series with size {}".format(
                    self.n_series
                )
            )
        except Exception:
            raise TypeError("Unsupported indexing type {}".format(type(idx)))

    def _restrict_to_interval(self, intervalslice):
        # Warning: This function can mutate data. It should only be called from
        # _restrict

        if isinstance(intervalslice, slice):
            if (
                intervalslice.start is None
                and intervalslice.stop is None
                and intervalslice.step is None
            ):
                # no restriction on interval
                return self

        newintervals = self._abscissa.support[intervalslice].merge()
        if newintervals.isempty:
            logging.warning("Index resulted in empty interval array")
            return self.empty(inplace=True)

        bcenter_inds = []
        bin_inds = []
        start = 0
        bsupport = np.zeros((newintervals.n_intervals, 2), dtype=int)
        support_intervals = np.zeros((newintervals.n_intervals, 2))

        if not self.isempty:
            for ii, interval in enumerate(newintervals.data):
                a_start = interval[0]
                a_stop = interval[1]
                frm, to = np.searchsorted(self._bins, (a_start, a_stop))
                # If bin edges equal a_stop, they should still be included
                if self._bins[to] <= a_stop:
                    bin_inds.extend(np.arange(frm, to + 1, step=1))
                else:
                    bin_inds.extend(np.arange(frm, to, step=1))
                    to -= 1
                support_intervals[ii] = [self._bins[frm], self._bins[to]]

                lind, rind = np.searchsorted(
                    self._bin_centers, (self._bins[frm], self._bins[to])
                )
                # We don't have to worry about an if-else block here unlike
                # for the bin_inds because the bin_centers can NEVER equal
                # the bins. Therefore we know every interval looks like
                # the following:
                #  first desired bin         last desired bin
                # |------------------|......|-------------------|
                #          ^                                         ^
                #          |                                         |
                #        lind                                      rind
                # Since arange is half-open, the indices we actually take
                # will be such that all bin centers fall within the desired
                # bin edges.
                bcenter_inds.extend(np.arange(lind, rind, step=1))

                bsupport[ii] = [start, start + (to - frm - 1)]
                start += to - frm

            self._bins = self._bins[bin_inds]
            self._bin_centers = self._bin_centers[bcenter_inds]
            self._data = np.atleast_2d(self._data[:, bcenter_inds])
            self._binned_support = bsupport

        self._abscissa.support = type(self._abscissa.support)(support_intervals)

    @property
    def isempty(self):
        """(bool) Empty BinnedEventArray."""
        try:
            return len(self.bin_centers) == 0
        except TypeError:
            return True  # this happens when self.bin_centers is None

    @property
    def n_series(self):
        """(int) The number of series."""
        try:
            return utils.PrettyInt(self.data.shape[0])
        except AttributeError:
            return 0

    @property
    def centers(self):
        """(np.array) The bin centers (in seconds)."""
        logging.warning("centers is deprecated. Use bin_centers instead.")
        return self.bin_centers

    @property
    def _abscissa_vals(self):
        """(np.array) The bin centers (in seconds)."""
        return self._bin_centers

    @property
    def bin_centers(self):
        """(np.array) The bin centers (in seconds)."""
        return self._bin_centers

    @property
    def event_centers(self):
        """(np.array) The centers (in seconds) of each event."""
        if self._event_centers is None:
            raise NotImplementedError("event_centers not yet implemented")
            # self._event_centers = midpoints
        return self._event_centers

    @property
    def _midpoints(self):
        """(np.array) The centers (in index space) of all events.

        Examples
        -------
        ax, img = npl.imagesc(bst.data) # data is of shape (n_series, n_bins)
        # then _midpoints correspond to the xvals at the center of
        # each event.
        ax.plot(bst.event_centers, np.repeat(1, self.n_intervals), marker='o', color='w')

        """
        if self._event_centers is None:
            midpoints = np.zeros(len(self.lengths))
            for idx, length in enumerate(self.lengths):
                midpoints[idx] = np.sum(self.lengths[:idx]) + length / 2
            self._event_centers = midpoints
        return self._event_centers

    @property
    def data(self):
        """(np.array) Event counts in all bins, with shape (n_series, n_bins)."""
        return self._data

    @property
    def bins(self):
        """(np.array) The bin edges (in seconds)."""
        return self._bins

    @property
    def binnedSupport(self):
        """(np.array) The binned support of the BinnedEventArray (in
        bin IDs) of shape (n_intervals, 2).
        """
        logging.warning("binnedSupport is deprecated. Use bined_support instead.")
        return self._binned_support

    @property
    def binned_support(self):
        """(np.array) The binned support of the BinnedEventArray (in
        bin IDs) of shape (n_intervals, 2).
        """
        return self._binned_support

    @property
    def lengths(self):
        """Lengths of contiguous segments, in number of bins."""
        if self.isempty:
            return 0
        return np.atleast_1d(
            (self.binned_support[:, 1] - self.binned_support[:, 0] + 1).squeeze()
        )

    @property
    def eventarray(self):
        """(nelpy.EventArray) The original EventArray associated with
        the binned data.
        """
        return self._eventarray

    @property
    def n_bins(self):
        """(int) The number of bins."""
        if self.isempty:
            return 0
        return utils.PrettyInt(len(self.bin_centers))

    @property
    def ds(self):
        """(float) Bin width in seconds."""
        return self._ds

    @ds.setter
    def ds(self, val):
        if self._ds is not None:
            raise AttributeError("can't set attribute")
        else:
            try:
                if val <= 0:
                    pass
            except ValueError:
                raise TypeError("bin width must be a scalar")
            if val <= 0:
                raise ValueError("bin width must be positive")
            self._ds = val

    @staticmethod
    def _get_bins_inside_interval(interval, ds):
        """(np.array) Return bin edges entirely contained inside an interval.

        Bin edges always start at interval.start, and continue for as many
        bins as would fit entirely inside the interval.

        NOTE 1: there are (n+1) bin edges associated with n bins.

        WARNING: if an interval is smaller than ds, then no bin will be
                associated with the particular interval.

        NOTE 2: nelpy uses half-open intervals [a,b), but if the bin
                width divides b-a, then the bins will cover the entire
                range. For example, if interval = [0,2) and ds = 1, then
                bins = [0,1,2], even though [0,2] is not contained in
                [0,2).

        Parameters
        ----------
        interval : IntervalArray
            IntervalArray containing a single interval with a start, and stop
        ds : float
            Time bin width, in seconds.

        Returns
        -------
        bins : array
            Bin edges in an array of shape (n+1,) where n is the number
            of bins
        centers : array
            Bin centers in an array of shape (n,) where n is the number
            of bins
        """

        if interval.length < ds:
            logging.warning("interval duration is less than bin size: ignoring...")
            return None, None

        n = int(np.floor(interval.length / ds))  # number of bins

        # linspace is better than arange for non-integral steps
        bins = np.linspace(interval.start, interval.start + n * ds, n + 1)
        centers = bins[:-1] + (ds / 2)
        return bins, centers

    def _bin_events(self, eventarray, intervalArray, ds):
        """
        Docstring goes here. TBD. For use with bins that are contained
        wholly inside the intervals.

        """
        b = []  # bin list
        c = []  # centers list
        s = []  # data list
        for nn in range(eventarray.n_series):
            s.append([])
        left_edges = []
        right_edges = []
        counter = 0
        for interval in intervalArray:
            bins, centers = self._get_bins_inside_interval(interval, ds)
            if bins is not None:
                for uu, eventarraydatas in enumerate(eventarray.data):
                    event_counts, _ = np.histogram(
                        eventarraydatas,
                        bins=bins,
                        density=False,
                        range=(interval.start, interval.stop),
                    )  # TODO: is it faster to limit range, or to cut out events?
                    s[uu].extend(event_counts.tolist())
                left_edges.append(counter)
                counter += len(centers) - 1
                right_edges.append(counter)
                counter += 1
                b.extend(bins.tolist())
                c.extend(centers.tolist())
        self._bins = np.array(b)
        self._bin_centers = np.array(c)
        self._data = np.array(s)
        le = np.array(left_edges)
        le = le[:, np.newaxis]
        re = np.array(right_edges)
        re = re[:, np.newaxis]
        self._binned_support = np.hstack((le, re))
        support_starts = self.bins[np.insert(np.cumsum(self.lengths + 1), 0, 0)[:-1]]
        support_stops = self.bins[np.insert(np.cumsum(self.lengths + 1) - 1, 0, 0)[1:]]
        supportdata = np.vstack([support_starts, support_stops]).T
        self._abscissa.support = type(self._abscissa.support)(
            supportdata
        )  # set support to TRUE bin support

    @keyword_deprecation(replace_x_with_y={"bw": "truncate"})
    def smooth(
        self, *, sigma=None, inplace=False, truncate=None, within_intervals=False
    ):
        """Smooth BinnedEventArray by convolving with a Gaussian kernel.

        Smoothing is applied in data, and the same smoothing is applied
        to each series in a BinnedEventArray.

        Smoothing is applied within each interval.

        Parameters
        ----------
        sigma : float, optional
            Standard deviation of Gaussian kernel, in seconds. Default is 0.01 (10 ms)
        truncate : float, optional
            Bandwidth outside of which the filter value will be zero. Default is 4.0
        inplace : bool
            If True the data will be replaced with the smoothed data.
            Default is False.

        Returns
        -------
        out : BinnedEventArray
            New BinnedEventArray with smoothed data.
        """

        if truncate is None:
            truncate = 4
        if sigma is None:
            sigma = 0.01  # 10 ms default

        fs = 1 / self.ds

        return utils.gaussian_filter(
            self,
            fs=fs,
            sigma=sigma,
            truncate=truncate,
            inplace=inplace,
            within_intervals=within_intervals,
        )

    @staticmethod
    def _smooth_array(arr, w=None):
        """Smooth an array by convolving a boxcar, row-wise.

        Parameters
        ----------
        w : int, optional
            Number of bins to include in boxcar window. Default is 10.

        Returns
        -------
        smoothed: array
            Smoothed array with same shape as arr.
        """

        if w is None:
            w = 10

        if w == 1:  # perform no smoothing
            return arr

        w = np.min((w, arr.shape[1]))

        smoothed = arr.astype(float)  # copy array and cast to float
        window = np.ones((w,)) / w

        # smooth per row
        for rowi, row in enumerate(smoothed):
            smoothed[rowi, :] = np.convolve(row, window, mode="same")

        if arr.shape[1] != smoothed.shape[1]:
            raise TypeError("Incompatible shape returned!")

        return smoothed

    @staticmethod
    def _rebin_array(arr, w):
        """Rebin an array of shape (n_signals, n_bins) into a
        coarser bin size.

        Parameters
        ----------
        arr : array
            Array with shape (n_signals, n_bins) to re-bin. A copy
            is returned.
        w : int
            Number of original bins to combine into each new bin.

        Returns
        -------
        out : array
            Bnned array with shape (n_signals, n_new_bins)
        bin_idx : array
            Array of shape (n_new_bins,) with the indices of the new
            binned array, relative to the original array.
        """
        cs = np.cumsum(arr, axis=1)
        binidx = np.arange(start=w, stop=cs.shape[1] + 1, step=w) - 1

        rebinned = np.hstack(
            (np.array(cs[:, w - 1], ndmin=2).T, cs[:, binidx[1:]] - cs[:, binidx[:-1]])
        )
        # bins = bins[np.insert(binidx+1, 0, 0)]
        return rebinned, binidx

    def rebin(self, w=None):
        """Rebin the BinnedEventArray into a coarser bin size.

        Parameters
        ----------
        w : int, optional
            number of bins of width bst.ds to bin into new bin of
            width bst.ds*w. Default is w=1 (no re-binning).

        Returns
        -------
        out : BinnedEventArray
            New BinnedEventArray with coarser resolution.
        """

        if w is None:
            w = 1

        if not float(w).is_integer:
            raise ValueError("w has to be an integer!")

        w = int(w)

        bst = self
        return self._rebin_binnedeventarray(bst, w=w)

    @staticmethod
    def _rebin_binnedeventarray(bst, w=None):
        """Rebin a BinnedEventArray into a coarser bin size.

        Parameters
        ----------
        bst : BinnedEventArray
            BinnedEventArray to re-bin into a coarser resolution.
        w : int, optional
            number of bins of width bst.ds to bin into new bin of
            width bst.ds*w. Default is w=1 (no re-binning).

        Returns
        -------
        out : BinnedEventArray
            New BinnedEventArray with coarser resolution.

        # FFB! TODO: if w is longer than some event size,
        # an exception will occur. Handle it! Although I may already
        # implicitly do that.
        """

        if w is None:
            w = 1

        if w == 1:
            return bst

        edges = np.insert(np.cumsum(bst.lengths), 0, 0)
        newlengths = [0]
        binedges = np.insert(np.cumsum(bst.lengths + 1), 0, 0)
        n_events = bst.support.n_intervals
        newdata = None

        for ii in range(n_events):
            data = bst.data[:, edges[ii] : edges[ii + 1]]
            bins = bst.bins[binedges[ii] : binedges[ii + 1]]

            datalen = data.shape[1]
            if w <= datalen:
                rebinned, binidx = bst._rebin_array(data, w=w)
                bins = bins[np.insert(binidx + 1, 0, 0)]

                newlengths.append(rebinned.shape[1])

                if newdata is None:
                    newdata = rebinned
                    newbins = bins
                    newcenters = bins[:-1] + np.diff(bins) / 2
                    newsupport = np.array([bins[0], bins[-1]])
                else:
                    newdata = np.hstack((newdata, rebinned))
                    newbins = np.hstack((newbins, bins))
                    newcenters = np.hstack((newcenters, bins[:-1] + np.diff(bins) / 2))
                    newsupport = np.vstack((newsupport, np.array([bins[0], bins[-1]])))
            else:
                pass

        # assemble new binned event series array:
        newedges = np.cumsum(newlengths)
        newbst = bst._copy_without_data()
        abscissa = copy.copy(bst._abscissa)
        if newdata is not None:
            newbst._data = newdata
            newbst._abscissa = abscissa
            newbst._abscissa.support = type(bst.support)(newsupport)
            newbst._bins = newbins
            newbst._bin_centers = newcenters
            newbst._ds = bst.ds * w
            newbst._binned_support = np.array((newedges[:-1], newedges[1:] - 1)).T
        else:
            logging.warning(
                "No events are long enough to contain any bins of width {}".format(
                    utils.PrettyDuration(bst.ds)
                )
            )
            newbst._data = None
            newbst._abscissa = abscissa
            newbst._abscissa.support = None
            newbst._binned_support = None
            newbst._bin_centers = None
            newbst._bins = None

        newbst.__renew__()

        return newbst

    def bst_from_indices(self, idx):
        """
        Return a BinnedEventArray from a list of indices.

        bst : BinnedEventArray
        idx : list of sample (bin) numbers with shape (n_intervals, 2) INCLUSIVE

        Examples
        --------
        idx = [[10, 20]
            [25, 50]]
        bst_from_indices(bst, idx=idx)
        """

        idx = np.atleast_2d(idx)

        newbst = self._copy_without_data()
        ds = self.ds
        bin_centers_ = []
        bins_ = []
        binned_support_ = []
        support_ = []
        all_abscissa_vals = []

        n_preceding_bins = 0

        for frm, to in idx:
            idx_array = np.arange(frm, to + 1).astype(int)
            all_abscissa_vals.append(idx_array)
            bin_centers = self.bin_centers[idx_array]
            bins = np.append(bin_centers - ds / 2, bin_centers[-1] + ds / 2)

            binned_support = [n_preceding_bins, n_preceding_bins + len(bins) - 2]
            n_preceding_bins += len(bins) - 1
            support = type(self._abscissa.support)((bins[0], bins[-1]))

            bin_centers_.append(bin_centers)
            bins_.append(bins)
            binned_support_.append(binned_support)
            support_.append(support)

        bin_centers = np.concatenate(bin_centers_)
        bins = np.concatenate(bins_)
        binned_support = np.array(binned_support_)
        support = np.sum(support_)
        all_abscissa_vals = np.concatenate(all_abscissa_vals)

        newbst._bin_centers = bin_centers
        newbst._bins = bins
        newbst._binned_support = binned_support
        newbst._abscissa.support = support
        newbst._data = newbst.data[:, all_abscissa_vals]

        newbst.__renew__()

        return newbst

    @property
    def n_active(self):
        """Number of active series.

        An active series is any series that fired at least one event.
        """
        if self.isempty:
            return 0
        return utils.PrettyInt(np.count_nonzero(self.n_events))

    @property
    def n_active_per_bin(self):
        """Number of active series per data bin with shape (n_bins,)."""
        if self.isempty:
            return 0
        # TODO: profile several alternatves. Could use data > 0, or
        # other numpy methods to get a more efficient implementation:
        return self.data.clip(max=1).sum(axis=0)

    @property
    def n_events(self):
        """(np.array) The number of events in each series."""
        if self.isempty:
            return 0
        return self.data.sum(axis=1)

    def flatten(self, *, series_id=None, series_label=None):
        """Collapse events across series.

        WARNING! series_tags are thrown away when flattening.

        Parameters
        ----------
        series_id: (int)
            (series) ID to assign to flattened event series, default is 0.
        series_label (str)
            (series) Label for event series, default is 'flattened'.
        """
        if self.n_series < 2:  # already flattened
            return self

        # default args:
        if series_id is None:
            series_id = 0
        if series_label is None:
            series_label = "flattened"

        binnedeventarray = self._copy_without_data()

        binnedeventarray._data = np.array(self.data.sum(axis=0), ndmin=2)

        binnedeventarray._bins = self.bins
        binnedeventarray._abscissa.support = self.support
        binnedeventarray._bin_centers = self.bin_centers
        binnedeventarray._binned_support = self.binned_support

        binnedeventarray._series_ids = [series_id]
        binnedeventarray._series_labels = [series_label]
        binnedeventarray._series_tags = None
        binnedeventarray.__renew__()

        return binnedeventarray

    @property
    def support(self):
        """(nelpy.IntervalArray) The support of the underlying BinnedEventArray."""
        return self._abscissa.support

    @support.setter
    def support(self, val):
        """(nelpy.IntervalArray) The support of the underlying BinnedEventArray."""
        # modify support
        if isinstance(val, type(self._abscissa.support)):
            self._abscissa.support = val
        elif isinstance(val, (tuple, list)):
            prev_domain = self._abscissa.domain
            self._abscissa.support = type(self._abscissa.support)([val[0], val[1]])
            self._abscissa.domain = prev_domain
        else:
            raise TypeError(
                "support must be of type {}".format(str(type(self._abscissa.support)))
            )
        # restrict data to new support
        self._restrict_to_interval(self._abscissa.support)