<gh_stars>0
import os

import numpy as np


def _deposit_or_erode(layers, n_layers, dz):
    """Update the array that contains layers with deposition or erosion.

    This function operates on the entire array that contain the layers (active,
    and allocated but not yet active). The number of active layers includes the
    layer that is currently being added. Thus the row with the index
    ``n_layers - 1`` is the layer that is currently being added as an active
    layer.

    Note that in EventLayers, layers represent an event, and not necessarily
    material.

    This means that if only erosion occurs, the array elements in the row with
    the index ``n_layers - 1`` will be set to zero and thickness will be
    removed from lower layers. Note that lower layers have smaller row indicies
    as the bottom of the layer stack has row index zero.

    If deposition occurs, the array elements in the row with index
    ``n_layers - 1`` will be set to the value of dz.

    Parameters
    ----------
    layers : ndarray of shape `(n_layers, n_nodes)`
        Array of layer thicknesses. This array is the datastructure that
        contains all allocated layers, active or inactive.
    n_layers : int
        Number of active layers.
    dz : ndarray of shape `(n_nodes, )`
        Thickness of the new layer. Negative thicknesses mean
        erode the top-most layers.

    Examples
    --------
    >>> import numpy as np
    >>> from landlab.layers.eventlayers import _deposit_or_erode

    First we create a numpy array allocated to contain layers. We fill it with
    -1. These -1.'s do not represent negative layer thickness. In EventLayers
    this array is created with np.empty, but that creates different numbers
    every time and doesn't work for testing.

    >>> allocated_layers_array = np.full((4, 3), 0.)

    Next we add a layer with spatially variable thickness. We specify that the
    number of active layers (including the one being added) is 1.

    >>> dz = np.array([1., 2., 3.])
    >>> _deposit_or_erode(allocated_layers_array, 1, dz)
    >>> allocated_layers_array
    array([[ 1.,  2.,  3.],
           [ 0.,  0.,  0.],
           [ 0.,  0.,  0.],
           [ 0.,  0.,  0.]])

    As you can see, this changes the value of the first row in the array. The
    remainder of the array represents space in the datatastructure that has
    been allocated to contain layers, but does not yet contain active layers.

    Next we add a layer of thickness 1. To do this, we now need to specify that
    the number of active layers is 2.

    >>> dz = np.array([1., 1., 1.])
    >>> _deposit_or_erode(allocated_layers_array, 2, dz)
    >>> allocated_layers_array
    array([[ 1.,  2.,  3.],
           [ 1.,  1.,  1.],
           [ 0.,  0.,  0.],
           [ 0.,  0.,  0.]])

    Finally, we do some erosion. We specify that the number of active layers is
    3 and give a spatially variable field of erosion and deposition.

    >>> _deposit_or_erode(allocated_layers_array, 3, [1., -1., -2.])
    >>> allocated_layers_array
    array([[ 1.,  2.,  2.],
           [ 1.,  0.,  0.],
           [ 1.,  0.,  0.],
           [ 0.,  0.,  0.]])

    >>> _deposit_or_erode(allocated_layers_array, 3, [1., -1., -2.])
    >>> allocated_layers_array
    array([[ 1.,  1.,  0.],
           [ 1.,  0.,  0.],
           [ 2.,  0.,  0.],
           [ 0.,  0.,  0.]])
    """
    from .ext.eventlayers import deposit_or_erode

    layers = layers.reshape((layers.shape[0], -1))
    try:
        dz = dz.reshape((layers.shape[1],))
    except (AttributeError, ValueError):
        dz = np.broadcast_to(dz, (layers.shape[1],))
    finally:
        dz = np.asfarray(dz)

    deposit_or_erode(layers, n_layers, dz)


def _get_surface_index(layers, n_layers, surface_index):
    """Get index within each stack of the layer at the topographic surface.

    Parameters
    ----------
    layers : ndarray of shape `(n_layers, n_nodes)`
        Array of layer thicknesses. This array is the datastructure that
        contains all allocated layers, active or inactive.
    n_layers : int
        Number of active layers.
    surface_index : ndarray of shape `(n_nodes, )`

    Examples
    --------
    >>> import numpy as np
    >>> from landlab.layers.eventlayers import (_deposit_or_erode,
    ...                                         _get_surface_index)

    >>> layers = np.full((5, 3), 1.)
    >>> dz = np.array([-1., -2., -3.])

    Note here, if you are very confused by the use of ``_deposit_or_erode``
    we recommend you read the docstring associated with that function.

    >>> _deposit_or_erode(layers, 5, dz)
    >>> layers
    array([[ 1.,  1.,  1.],
           [ 1.,  1.,  1.],
           [ 1.,  1.,  0.],
           [ 1.,  0.,  0.],
           [ 0.,  0.,  0.]])

    >>> surface_index = np.empty(3, dtype=int)
    >>> _get_surface_index(layers, 5, surface_index)
    >>> surface_index
    array([3, 2, 1])
    """
    from .ext.eventlayers import get_surface_index

    layers = layers.reshape((layers.shape[0], -1))

    get_surface_index(layers, n_layers, surface_index)


def _reduce_matrix(array, step, reducer):
    """Combine rows of a 2D matrix.

    Parameters
    ----------
    array : ndarray, shape (m, n)
        Matrix to reduce.
    step : int
        Number of rows in each block to reduce.
    reducer : ufunc
        Function to use for combining rows.

    Returns
    -------
    ndarray
        Matrix with rows combined.

    Examples
    --------
    >>> import numpy as np
    >>> from landlab.layers.eventlayers import _reduce_matrix
    >>> array = np.arange(12).reshape((4, 3))
    >>> array
    array([[ 0,  1,  2],
           [ 3,  4,  5],
           [ 6,  7,  8],
           [ 9, 10, 11]])
    >>> _reduce_matrix(array, 4, np.sum)
    array([[18, 22, 26]])
    >>> _reduce_matrix(array, 2, np.sum)
    array([[ 3,  5,  7],
           [15, 17, 19]])
    """
    return reducer(array.reshape((-1, step) + array.shape[1:]), axis=1).reshape(
        (-1,) + array.shape[1:]
    )


class _BlockSlice:
    """Slices that divide a matrix into equally sized blocks."""

    def __init__(self, *args):
        """_BlockSlice([start], stop, [step])"""
        if len(args) > 3:
            raise TypeError(
                "_BlockSlice expected at most 3 arguments, got {0}".format(len(args))
            )

        self._args = tuple(args)

        self._start = 0
        self._stop = None
        self._step = None

        if len(args) == 1:
            self._stop = args[0]
        elif len(args) == 2:
            self._start, self._stop = args
        elif len(args) == 3:
            self._start, self._stop, self._step = args
            self._step = min(self._stop - self._start, self._step)

        if self._stop is not None and self._stop < self._start:
            raise ValueError(
                "stop ({0}) must be greater than start ({1})".format(
                    self._stop, self._start
                )
            )

    def __repr__(self):
        return "_BlockSlice({0})".format(", ".join([repr(arg) for arg in self._args]))

    @property
    def start(self):
        return self._start

    @property
    def stop(self):
        return self._stop

    @property
    def step(self):
        return self._step

    def indices(self, n_rows):
        """Row indices to blocks within a matrix.

        Parameters
        ----------
        n_rows : int
            The number of rows in the matrix.

        Returns
        -------
        (start, stop, step)
            Tuple of (int* that gives the row of the first block, row of the
            last block, and the number of rows in each block.

        Examples
        --------
        >>> from landlab.layers.eventlayers import _BlockSlice

        The default is one single block that encomapses all the rows.

        >>> _BlockSlice().indices(4)
        (0, 4, 4)

        >>> _BlockSlice(3).indices(4)
        (0, 3, 3)

        >>> _BlockSlice(1, 3).indices(4)
        (1, 3, 2)

        >>> _BlockSlice(1, 7, 2).indices(8)
        (1, 7, 2)
        """
        start, stop, step = self.start, self.stop, self.step
        if stop is None:
            stop = n_rows

        start, stop, _ = slice(start, stop).indices(n_rows)

        if step is None:
            step = stop - start

        if step != 0 and (stop - start) % step != 0:
            stop = (stop - start) // step * step + start

        return start, stop, step


def _valid_keywords_or_raise(kwds, required=(), optional=()):
    """Check for valid keyword arguments.

    Parameters
    ----------
    kwds : iterable of str
        Keywords to check for validity.
    required : iterable of str, optional
        Keywords that are required.
    optional : iterable of str, optional
        Keywords that are optional.

    Examples
    --------
    >>> from landlab.layers.eventlayers import _valid_keywords_or_raise
    >>> _valid_keywords_or_raise(["foo"], optional=["foo", "bar"])
    >>> _valid_keywords_or_raise(["foo"], required=["foo", "bar"])
    Traceback (most recent call last):
        ...
    TypeError: missing keyword arguments ('bar')
    >>> _valid_keywords_or_raise(["baz"], optional=["foo", "bar"])
    Traceback (most recent call last):
        ...
    TypeError: invalid keyword arguments ('baz' not in {'bar', 'foo'})
    """
    keys = set(kwds)
    required = set(required)
    optional = required | set(optional)

    unknown = keys - optional
    if unknown:
        raise TypeError(
            "invalid keyword arguments ({0} not in {{{1}}})".format(
                ", ".join(sorted([repr(name) for name in unknown])),
                ", ".join(sorted([repr(name) for name in optional])),
            )
        )

    missing = required - keys
    if missing:
        raise TypeError(
            "missing keyword arguments ({0})".format(
                ", ".join(sorted([repr(name) for name in missing]))
            )
        )


def resize_array(array, newsize, exact=False):
    """Increase the size of an array, leaving room to grow.

    Parameters
    ----------
    array : ndarray
        The array to resize.
    newsize : int
        Size of the zero-th dimension of the resized array.
    exact : bool, optional
        Should the new array have the exact size provided or
        at least that size.

    Returns
    -------
    ndarray
        Copy of the input array resized.

    Examples
    --------
    >>> import numpy as np
    >>> from landlab.layers.eventlayers import resize_array

    >>> x = np.arange(6)
    >>> bigger_x = resize_array(x, 10)
    >>> bigger_x.size
    17
    >>> np.all(x[:6] == bigger_x[:6])
    True
    >>> x is bigger_x
    False

    >>> x = np.arange(6).reshape((2, 3))
    >>> bigger_x = resize_array(x, 4)
    >>> bigger_x.shape == (10, 3)
    True

    >>> bigger_x = resize_array(x, 4, exact=True)
    >>> bigger_x.shape == (4, 3)
    True

    :meta private:
    """
    newsize = int(newsize)
    allocated = array.shape[0]

    if newsize <= allocated:
        return array

    if exact:
        new_allocated = newsize
    else:
        new_allocated = (newsize >> 3) + 6 + newsize

    larger_array = np.empty((new_allocated,) + array.shape[1:], dtype=array.dtype)
    larger_array[:allocated] = array

    return larger_array


def _allocate_layers_for(array, number_of_layers, number_of_stacks):
    """Allocate a layer matrix.

    Parameters
    ----------
    array : number or ndarray
        Array of layer properties to track.
    number_of_layers : int
        Number of layers to allocate.
    number_of_stacks : int
        Number of stacks to allocate.

    Returns
    -------
    ndarray of size `(number_of_layers, number_of_stacks, values_per_stack)`
        Newly allocated matrix for storing layer properties.

    Examples
    --------
    >>> import numpy as np
    >>> from landlab.layers.eventlayers import _allocate_layers_for

    >>> layers = _allocate_layers_for(3, 2, 4)
    >>> layers.shape == (2, 4)
    True
    >>> layers.dtype.kind == 'i'
    True

    >>> layers = _allocate_layers_for(np.zeros(4), 2, 4)
    >>> layers.shape == (2, 4)
    True
    >>> layers.dtype.kind == 'f'
    True

    >>> layers = _allocate_layers_for(np.zeros(2), 2, 4)
    >>> layers.shape == (2, 4, 2)
    True
    >>> layers.dtype.kind == 'f'
    True
    """
    array = np.asarray(array)

    if array.ndim > 0 and len(array) != number_of_stacks:
        values_per_stack = array.shape
    else:
        values_per_stack = array.shape[1:]

    return np.empty(
        (number_of_layers, number_of_stacks) + values_per_stack, dtype=array.dtype
    )


class EventLayersMixIn:

    """MixIn that adds a EventLayers attribute to a ModelGrid."""

    @property
    def event_layers(self):
        """EventLayers for each cell."""
        try:
            self._event_layers
        except AttributeError:
            self._event_layers = EventLayers(self.number_of_cells)
        finally:
            return self._event_layers

    @property
    def at_layer(self):
        """EventLayers for each cell."""
        return self.event_layers


class EventLayers:

    """Track EventLayers where each event is its own layer.

    EventLayers are meant to represent a layered object in which each layer
    represents a event. Thus they are likely the most appropriate tool to use
    if the user is interested in chronostratigraphy. If erosion occurs, a new
    layer with zero thickness is created. Thus, EventLayers may not be the most
    memory efficent layers datastructure.

    EventLayers exists in contrast to the MaterialLayers object which does not
    make a new layer if only erosion occurs and if the attributes of the new
    layer are equivalent to the attributes of the material at the surface of the
    layer stack.

    Parameters
    ----------
    number_of_stacks : int
        Number of layer stacks to track.

    Examples
    --------
    >>> from landlab.layers.eventlayers import EventLayers

    Create an empty layer stack with 5 stacks.

    >>> layers = EventLayers(5)
    >>> layers.number_of_stacks
    5
    >>> layers.number_of_layers
    0

    Add a layer with a uniform thickness.

    >>> layers.add(1.5)
    >>> layers.number_of_layers
    1
    >>> layers.dz
    array([[ 1.5,  1.5,  1.5,  1.5,  1.5]])

    Add a second layer with uneven thickness.

    >>> layers.add([1., 2., .5, 5., 0.])
    >>> layers.dz
    array([[ 1.5,  1.5,  1.5,  1.5,  1.5],
           [ 1. ,  2. ,  0.5,  5. ,  0. ]])

    Adding a layer with negative thickness will remove
    existing layers for the top of the stack. Note that
    this will create a new layer with thickness zero
    that represents this 'event'. If instead your
    application would prefer that no new row is added to
    the layers datastructure, you may want to consider
    the MaterialLayers object.

    >>> layers.add(-1)
    >>> layers.dz
    array([[ 1.5,  1.5,  1. ,  1.5,  0.5],
           [ 0. ,  1. ,  0. ,  4. ,  0. ],
           [ 0. ,  0. ,  0. ,  0. ,  0. ]])

    Get the index value of the layer within each stack
    at the topographic surface.

    >>> layers.surface_index
    array([0, 1, 0, 1, 0])
    """

    def __init__(self, number_of_stacks, allocated=0):
        self._number_of_layers = 0
        self._number_of_stacks = number_of_stacks
        self._surface_index = np.zeros(number_of_stacks, dtype=int)
        self._attrs = dict()

        dims = (self.number_of_layers, self.number_of_stacks)
        self._attrs["_dz"] = np.empty(dims, dtype=float)
        self._resize(allocated, exact=True)

    def __getitem__(self, name):
        return self._attrs[name][: self.number_of_layers]

    def __setitem__(self, name, values):
        dims = (self.allocated, self.number_of_stacks)
        values = np.asarray(values)
        if values.ndim == 1:
            values = np.expand_dims(values, 1)
        values = np.broadcast_to(values, (self.number_of_layers, self.number_of_stacks))
        self._attrs[name] = _allocate_layers_for(values.flatten()[0], *dims)
        self._attrs[name][: self.number_of_layers] = values

    def __iter__(self):
        return (name for name in self._attrs if not name.startswith("_"))

    def __str__(self):
        lines = [
            "number_of_layers: {number_of_layers}",
            "number_of_stacks: {number_of_stacks}",
            "tracking: {attrs}",
        ]
        return os.linesep.join(lines).format(
            number_of_layers=self.number_of_layers,
            number_of_stacks=self.number_of_stacks,
            attrs=", ".join(self.tracking) or "null",
        )

    def __repr__(self):
        return self.__class__.__name__ + "({number_of_stacks})".format(
            number_of_stacks=self.number_of_stacks
        )

    @property
    def tracking(self):
        """Layer properties being tracked.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers
        >>> layers = EventLayers(3)
        >>> layers.tracking
        []
        >>> layers.add(1., age=1.)
        >>> layers.tracking
        ['age']
        """
        return [name for name in self._attrs if not name.startswith("_")]

    def _setup_layers(self, **kwds):
        dims = (self.allocated, self.number_of_stacks)
        for name, array in kwds.items():
            self._attrs[name] = _allocate_layers_for(array, *dims)

    @property
    def number_of_stacks(self):
        """Number of stacks."""
        return self._number_of_stacks

    @property
    def thickness(self):
        """Total thickness of the columns.

        The sum of all layer thicknesses for each stack as an array
        of shape `(number_of_stacks, )`.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        Initially there are no layers so the total thickness is 0.

        >>> layers = EventLayers(3)
        >>> layers.thickness
        array([ 0.,  0.,  0.])

        After adding some layers, the stacks have varying thicknesses.

        >>> layers.add(15.)
        >>> layers.add([1., -1., 2.])
        >>> layers.thickness
        array([ 16.,  14.,  17.])
        """
        return np.sum(self.dz, axis=0)

    @property
    def z(self):
        """Thickness to top of each layer.

        Thickness from the bottom of each stack to the top of each layer
        as an array of shape `(number_of_layers, number_of_stacks)`.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        Initially there are no layers so the elevation to the top
        is 0.

        >>> layers = EventLayers(3)
        >>> layers.z.shape == (0, 3)
        True

        After adding some layers, elevations are to the top of each layer.

        >>> layers.add(15.)
        >>> layers.add([1., -1., 2.])
        >>> layers.dz
        array([[ 15.,  14.,  15.],
               [  1.,   0.,   2.]])
        >>> layers.z
        array([[ 15.,  14.,  15.],
               [ 16.,  14.,  17.]])
        """
        return np.cumsum(self.dz, axis=0)

    @property
    def dz(self):
        """Thickness of each layer.

        The thickness of each layer at each stack as an array of shape
        `(number_of_layers, number_of_stacks)`.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        Initially there are no layers so there are not thicknesses.

        >>> layers = EventLayers(3)
        >>> layers.dz.shape == (0, 3)
        True

        Now add two layers, the first of uniform thickness and the
        second non-uniform and with some erosion.

        >>> layers.add(15.)
        >>> layers.add([1., -1., 2.])
        >>> layers.dz
        array([[ 15.,  14.,  15.],
               [  1.,   0.,   2.]])
        """
        return self._attrs["_dz"][: self.number_of_layers]

    @property
    def number_of_layers(self):
        """Total number of layers.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        >>> layers = EventLayers(3)
        >>> layers.number_of_layers
        0

        >>> layers.add(15.)
        >>> layers.add([1., -1., 2.])
        >>> layers.number_of_layers
        2
        """
        return self._number_of_layers

    @property
    def allocated(self):
        """Total number of allocated layers.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        >>> layers = EventLayers(3)
        >>> layers.number_of_layers
        0
        >>> layers.allocated == 0
        True

        >>> layers.add(15.)
        >>> layers.number_of_layers
        1
        >>> layers.allocated == 7
        True
        >>> for _ in range(layers.allocated): layers.add(0.)
        >>> layers.number_of_layers
        8
        >>> layers.allocated == 15
        True

        If you know how many layers you will ultimately have, you
        can allocated enough memory for them when you create your
        layer stacks.

        >>> layers = EventLayers(3, allocated=15)
        >>> layers.number_of_layers
        0
        >>> layers.allocated == 15
        True

        >>> layers.add(15.)
        >>> layers.number_of_layers
        1
        >>> layers.allocated == 15
        True
        >>> for _ in range(layers.allocated): layers.add(0.)
        >>> layers.number_of_layers
        16
        >>> layers.allocated == 24
        True
        """
        return self._attrs["_dz"].shape[0]

    def add(self, dz, **kwds):
        """Add a layer to the stacks.

        Parameters
        ----------
        dz : float or array_like
            Thickness to add to each stack.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        Create an empty layer stack with 3 stacks.

        >>> layers = EventLayers(3)
        >>> layers.number_of_layers
        0

        To add a layer of uniform thickness to every stack.

        >>> layers.add(1.5)
        >>> layers.number_of_layers
        1
        >>> layers.dz
        array([[ 1.5,  1.5,  1.5]])

        Add a second layer with uneven thickness.

        >>> layers.add([1., 2., .5])
        >>> layers.dz
        array([[ 1.5,  1.5,  1.5],
               [ 1. ,  2. ,  0.5]])

        Adding a layer with negative thickness will remove
        existing layers for the top of the stack.

        >>> layers.add(-1)
        >>> layers.dz
        array([[ 1.5,  1.5,  1. ],
               [ 0. ,  1. ,  0. ],
               [ 0. ,  0. ,  0. ]])

        Use keywords to track properties of each layer. For instance,
        here we create a new stack and add a layer with a particular
        *age*. You can access the layer properties as if the object
        were a dictionary.

        >>> layers = EventLayers(3)
        >>> layers.add(1., age=3.)
        >>> layers.dz
        array([[ 1.,  1.,  1.]])
        >>> layers['age']
        array([[ 3.,  3.,  3.]])
        >>> layers.add(2., age=6.)
        >>> layers['age']
        array([[ 3.,  3.,  3.],
               [ 6.,  6.,  6.]])

        Attributes for each layer will exist even if the the layer is
        associated with erosion.

        >>> layers.add([-2, -1, 1], age=8.)
        >>> layers.dz
        array([[ 1.,  1.,  1.],
               [ 0.,  1.,  2.],
               [ 0.,  0.,  1.]])
        >>> layers['age']
        array([[ 3.,  3.,  3.],
               [ 6.,  6.,  6.],
               [ 8.,  8.,  8.]])

        To get the values at the surface of the layer stack:

        >>> layers.get_surface_values('age')
        array([ 3.,  6.,  8.])
        """
        if self.number_of_layers == 0:
            self._setup_layers(**kwds)

        self._add_empty_layer()

        _deposit_or_erode(self._attrs["_dz"], self.number_of_layers, dz)
        _get_surface_index(
            self._attrs["_dz"], self.number_of_layers, self._surface_index
        )

        for name in kwds:
            try:
                self[name][-1] = kwds[name]
            except KeyError:
                raise ValueError(
                    "EventLayers: {0} is not being tracked. Error in adding.".format(
                        name
                    )
                )

    def reduce(self, *args, **kwds):
        """reduce([start], stop, [step])
        Combine layers.

        Reduce adjacent layers into a single layer.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        Create an empty layer stack with 3 stacks.

        >>> layers = EventLayers(3)
        >>> layers.number_of_layers
        0

        To add a layer of uniform thickness to every stack.

        >>> layers.add(1.5)
        >>> layers.number_of_layers
        1
        >>> layers.dz
        array([[ 1.5,  1.5,  1.5]])

        Add a second layer with uneven thickness.

        >>> layers.add([1., 2., .5])
        >>> layers.dz
        array([[ 1.5,  1.5,  1.5],
               [ 1. ,  2. ,  0.5]])

        Combine all of the layers into a single layer.

        >>> layers.reduce()
        >>> layers.dz
        array([[ 2.5,  3.5,  2. ]])

        Add two additional layers to the top. The bottom-most layer is row
        0, and the two new layers are rows 1 and 2.

        >>> layers.add([1., 2., .5])
        >>> layers.add([1., 2., .5])
        >>> layers.dz
        array([[ 2.5,  3.5,  2. ],
               [ 1. ,  2. ,  0.5],
               [ 1. ,  2. ,  0.5]])

        Combine the two new layers (layers 1 and 2) into a single layer.

        >>> layers.reduce(1, 3)
        >>> layers.dz
        array([[ 2.5,  3.5,  2. ],
               [ 2. ,  4. ,  1. ]])

        >>> layers.add([1., 2., .5])
        >>> layers.add([1., 2., .5])
        >>> layers.dz
        array([[ 2.5,  3.5,  2. ],
               [ 2. ,  4. ,  1. ],
               [ 1. ,  2. ,  0.5],
               [ 1. ,  2. ,  0.5]])

        Combine the middle two layers.

        >>> layers.reduce(1, 3)
        >>> layers.dz
        array([[ 2.5,  3.5,  2. ],
               [ 3. ,  6. ,  1.5],
               [ 1. ,  2. ,  0.5]])
        >>> layers.add([1., 1., 1.])
        >>> layers.dz
        array([[ 2.5,  3.5,  2. ],
               [ 3. ,  6. ,  1.5],
               [ 1. ,  2. ,  0.5],
               [ 1. ,  1. ,  1. ]])

        Combine every two layers (layers 0 and 1 and combined, and layers
        1 and 2 are combined).

        >>> layers.reduce(0, 4, 2)
        >>> layers.dz
        array([[ 5.5,  9.5,  3.5],
               [ 2. ,  3. ,  1.5]])

        When layers are combined, thicknesses are summed but layer attributes
        can be combined in other ways (e.g. max, or mean)

        >>> layers = EventLayers(3)
        >>> layers.add([1, 1, 1], age=0.)
        >>> layers.add([1, 2, 5], age=1.)
        >>> layers.add([2, 2, 2], age=2.)
        >>> layers.reduce(age=np.max)
        >>> layers["age"]
        array([[ 2.,  2.,  2.]])

        >>> layers.add([2, 2, 2], age=3.)
        >>> layers.add([2, 2, 2], age=4.)
        >>> layers.reduce(1, 3, age=np.mean)
        >>> layers["age"]
        array([[ 2. ,  2. ,  2. ],
               [ 3.5,  3.5,  3.5]])
        """
        _valid_keywords_or_raise(kwds, required=self.tracking, optional=self._attrs)

        start, stop, step = _BlockSlice(*args).indices(self._number_of_layers)

        if step <= 1:
            return

        n_blocks = (stop - start) // step
        n_removed = n_blocks * (step - 1)
        for name, array in self._attrs.items():
            middle = _reduce_matrix(array[start:stop, :], step, kwds.get(name, np.sum))
            top = array[stop : self._number_of_layers, :]

            array[start : start + n_blocks, :] = middle
            array[start + n_blocks : start + n_blocks + len(top)] = top

        self._number_of_layers -= n_removed
        self._surface_index[:] -= n_removed

    @property
    def surface_index(self):
        """Index to the top non-empty layer.

        Examples
        --------
        >>> from landlab.layers.eventlayers import EventLayers

        Create an empty layer stack with 5 stacks.

        >>> layers = EventLayers(3)
        >>> layers.surface_index
        array([0, 0, 0])

        Add a layer with a uniform thickness.

        >>> for _ in range(5): layers.add(1.0)
        >>> layers.surface_index
        array([4, 4, 4])

        Add a layer with varying thickness. Negative thickness
        removes thickness from underlying layers, zero thickness adds a
        layer but doesn't change the surface index.

        >>> layers.add([-1.0, 0.0, 1.0])
        >>> layers.surface_index
        array([3, 4, 5])
        """
        return self._surface_index

    def get_surface_values(self, name):
        """Values of a field on the surface layer."""
        return self._attrs[name][self.surface_index, np.arange(self._number_of_stacks)]

    def _add_empty_layer(self):
        """Add a new empty layer to the stacks."""
        if self.number_of_layers >= self.allocated:
            self._resize(self.allocated + 1)

        self._number_of_layers += 1
        self._attrs["_dz"][self.number_of_layers - 1, :] = 0.0
        for name in self._attrs:
            self._attrs[name][self.number_of_layers - 1] = 0.0

    def _resize(self, newsize, exact=False):
        """Allocate more memory for the layers."""
        for name in self._attrs:
            self._attrs[name] = resize_array(self._attrs[name], newsize, exact=exact)
