ENH: improved, faster algorithm for array padding

New padding method which scales much better with dimensionality.
This new implementation is fully vectorized, builds each abstracted
n-dimensional padding block in a single step, and takes advantage
of separability. The API is completely preserved, and the old
algorithm is used if a vector function is input for `mode`.

The new algorithm is faster for all tested combinations of inputs,
and scales much better with dimensionality. Execution time reductions
from ~25% for small rank 1 arrays to >99% for rank 4+ arrays observed.

1 files changed, 1087 insertions, 414 deletions

diff --git a/numpy/lib/arraypad.py b/numpy/lib/arraypad.py
index 09ca99332..c80f6f54d 100644
--- a/numpy/lib/arraypad.py
+++ b/numpy/lib/arraypad.py
@@ -15,508 +15,1075 @@ __all__ = ['pad']
 # Private utility functions.
 
 
-def _create_vector(vector, pad_tuple, before_val, after_val):
+def _arange_ndarray(arr, shape, axis, reverse=False):
     """
-    Private function which creates the padded vector.
+    Create an ndarray of `shape` with increments along specified `axis`
 
     Parameters
     ----------
-    vector : ndarray of rank 1, length N + pad_tuple[0] + pad_tuple[1]
-        Input vector including blank padded values. `N` is the lenth of the
-        original vector.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding along
-        this particular iaxis.
-    before_val : scalar or ndarray of rank 1, length pad_tuple[0]
-        This is the value(s) that will pad the beginning of `vector`.
-    after_val : scalar or ndarray of rank 1, length pad_tuple[1]
-        This is the value(s) that will pad the end of the `vector`.
+    arr : ndarray
+        Input array of arbitrary shape.
+    shape : tuple of ints
+        Shape of desired array. Should be equivalent to `arr.shape` except
+        `shape[axis]` which may have any positive value.
+    axis : int
+        Axis to increment along.
+    reverse : bool
+        If False, increment in a positive fashion from 1 to `shape[axis]`,
+        inclusive. If True, the bounds are the same but the order reversed.
 
     Returns
     -------
-    _create_vector : ndarray
-        Vector with before_val and after_val replacing the blank pad values.
+    padarr : ndarray
+        Output array sized to pad `arr` along `axis`, with linear range from
+        1 to `shape[axis]` along specified `axis`.
+
+    Notes
+    -----
+    The range is deliberately 1-indexed for this specific use case. Think of
+    this algorithm as broadcasting `np.arange` to a single `axis` of an
+    arbitrarily shaped ndarray.
 
     """
-    vector[:pad_tuple[0]] = before_val
-    if pad_tuple[1] > 0:
-        vector[-pad_tuple[1]:] = after_val
-    return vector
+    initshape = tuple(1 if i != axis else shape[axis]
+                       for (i, x) in enumerate(arr.shape))
+    if not reverse:
+        padarr = np.arange(1, shape[axis] + 1)
+    else:
+        padarr = np.arange(shape[axis], 0, -1)
+    padarr = padarr.reshape(initshape)
+    for i, dim in enumerate(shape):
+        if padarr.shape[i] != dim:
+            padarr = padarr.repeat(dim, axis=i)
+    return padarr
 
 
-def _normalize_shape(narray, shape):
+def _round_ifneeded(arr, dtype):
     """
-    Private function which does some checks and normalizes the possibly
-    much simpler representations of 'pad_width', 'stat_length',
-    'constant_values', 'end_values'.
+    Rounds arr inplace if destination dtype is integer.
 
     Parameters
     ----------
-    narray : ndarray
-        Input ndarray
-    shape : {sequence, int}, optional
-        The width of padding (pad_width) or the number of elements on the
-        edge of the narray used for statistics (stat_length).
-        ((before_1, after_1), ... (before_N, after_N)) unique number of
-        elements for each axis where `N` is rank of `narray`.
-        ((before, after),) yields same before and after constants for each
-        axis.
-        (constant,) or int is a shortcut for before = after = constant for
-        all axes.
+    arr : ndarray
+        Input array.
+    dtype : dtype
+        The dtype of the destination array.
+
+    """
+    if np.issubdtype(dtype, np.integer):
+        arr.round(out=arr)
+
+
+def _prepend_const(arr, pad_amt, val, axis=-1):
+    """
+    Prepend constant `val` along `axis` of `arr`.
+
+    Parameters
+    ----------
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to prepend.
+    val : scalar
+        Constant value to use. For best results should be of type `arr.dtype`;
+        if not `arr.dtype` will be cast to `arr.dtype`.
+    axis : int
+        Axis along which to pad `arr`.
 
     Returns
     -------
-    _normalize_shape : tuple of tuples
-        int                               => ((int, int), (int, int), ...)
-        [[int1, int2], [int3, int4], ...] => ((int1, int2), (int3, int4), ...)
-        ((int1, int2), (int3, int4), ...) => no change
-        [[int1, int2], ]                  => ((int1, int2), (int1, int2), ...)
-        ((int1, int2), )                  => ((int1, int2), (int1, int2), ...)
-        [[int ,     ], ]                  => ((int, int), (int, int), ...)
-        ((int ,     ), )                  => ((int, int), (int, int), ...)
+    padarr : ndarray
+        Output array, with `pad_amt` constant `val` prepended along `axis`.
 
     """
-    normshp = None
-    shapelen = len(np.shape(narray))
-    if (isinstance(shape, int)):
-        normshp = ((shape, shape), ) * shapelen
-    elif (isinstance(shape, (tuple, list))
-            and isinstance(shape[0], (tuple, list))
-            and len(shape) == shapelen):
-        normshp = shape
-        for i in normshp:
-            if len(i) != 2:
-                fmt = "Unable to create correctly shaped tuple from %s"
-                raise ValueError(fmt % (normshp,))
-    elif (isinstance(shape, (tuple, list))
-            and isinstance(shape[0], (int, float, long))
-            and len(shape) == 1):
-        normshp = ((shape[0], shape[0]), ) * shapelen
-    elif (isinstance(shape, (tuple, list))
-            and isinstance(shape[0], (int, float, long))
-            and len(shape) == 2):
-        normshp = (shape, ) * shapelen
-    if normshp is None:
-        fmt = "Unable to create correctly shaped tuple from %s"
-        raise ValueError(fmt % (shape,))
-    return normshp
+    if pad_amt == 0:
+        return arr
+    padshape = tuple(x if i != axis else pad_amt
+                     for (i, x) in enumerate(arr.shape))
+    if val == 0:
+        return np.concatenate((np.zeros(padshape, dtype=arr.dtype), arr),
+                              axis=axis)
+    else:
+        return np.concatenate(((np.zeros(padshape) + val).astype(arr.dtype),
+                               arr), axis=axis)
 
 
-def _validate_lengths(narray, number_elements):
+def _append_const(arr, pad_amt, val, axis=-1):
     """
-    Private function which does some checks and reformats pad_width and
-    stat_length using _normalize_shape.
+    Append constant `val` along `axis` of `arr`.
 
     Parameters
     ----------
-    narray : ndarray
-        Input ndarray
-    number_elements : {sequence, int}, optional
-        The width of padding (pad_width) or the number of elements on the edge
-        of the narray used for statistics (stat_length).
-        ((before_1, after_1), ... (before_N, after_N)) unique number of
-        elements for each axis.
-        ((before, after),) yields same before and after constants for each
-        axis.
-        (constant,) or int is a shortcut for before = after = constant for all
-        axes.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to append.
+    val : scalar
+        Constant value to use. For best results should be of type `arr.dtype`;
+        if not `arr.dtype` will be cast to `arr.dtype`.
+    axis : int
+        Axis along which to pad `arr`.
 
     Returns
     -------
-    _validate_lengths : tuple of tuples
-        int                               => ((int, int), (int, int), ...)
-        [[int1, int2], [int3, int4], ...] => ((int1, int2), (int3, int4), ...)
-        ((int1, int2), (int3, int4), ...) => no change
-        [[int1, int2], ]                  => ((int1, int2), (int1, int2), ...)
-        ((int1, int2), )                  => ((int1, int2), (int1, int2), ...)
-        [[int ,     ], ]                  => ((int, int), (int, int), ...)
-        ((int ,     ), )                  => ((int, int), (int, int), ...)
+    padarr : ndarray
+        Output array, with `pad_amt` constant `val` appended along `axis`.
 
     """
-    normshp = _normalize_shape(narray, number_elements)
-    for i in normshp:
-        if i[0] < 0 or i[1] < 0:
-            fmt = "%s cannot contain negative values."
-            raise ValueError(fmt % (number_elements,))
-    return normshp
+    if pad_amt == 0:
+        return arr
+    padshape = tuple(x if i != axis else pad_amt
+                     for (i, x) in enumerate(arr.shape))
+    if val == 0:
+        return np.concatenate((arr, np.zeros(padshape, dtype=arr.dtype)),
+                              axis=axis)
+    else:
+        return np.concatenate(
+            (arr, (np.zeros(padshape) + val).astype(arr.dtype)), axis=axis)
 
 
-def _create_stat_vectors(vector, pad_tuple, iaxis, kwargs):
+def _prepend_edge(arr, pad_amt, axis=-1):
     """
-    Returns the portion of the vector required for any statistic.
+    Prepend `pad_amt` to `arr` along `axis` by extending edge values.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.
-    kwargs : keyword arguments
-        Keyword arguments. Only 'stat_length' is used. 'stat_length'
-        defaults to the entire vector if not supplied.
-
-    Return
-    ------
-    _create_stat_vectors : ndarray
-        The values from the original vector that will be used to calculate
-        the statistic.
-
-    """
-
-    # Can't have 0 represent the end if a slice... a[1:0] doesnt' work
-    pt1 = -pad_tuple[1]
-    if pt1 == 0:
-        pt1 = None
-
-    # Default is the entire vector from the original array.
-    sbvec = vector[pad_tuple[0]:pt1]
-    savec = vector[pad_tuple[0]:pt1]
-
-    if kwargs['stat_length']:
-        stat_length = kwargs['stat_length'][iaxis]
-        sl0 = min(stat_length[0], len(sbvec))
-        sl1 = min(stat_length[1], len(savec))
-        sbvec = np.arange(0)
-        savec = np.arange(0)
-        if pad_tuple[0] > 0:
-            sbvec = vector[pad_tuple[0]:pad_tuple[0] + sl0]
-        if pad_tuple[1] > 0:
-            savec = vector[-pad_tuple[1] - sl1:pt1]
-    return (sbvec, savec)
-
-
-def _maximum(vector, pad_tuple, iaxis, kwargs):
-    """
-    Private function to calculate the before/after vectors for pad_maximum.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to prepend.
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, extended by `pad_amt` edge values appended along `axis`.
+
+    """
+    if pad_amt == 0:
+        return arr
+
+    edge_slice = tuple(slice(None) if i != axis else 0
+                       for (i, x) in enumerate(arr.shape))
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+    edge_arr = arr[edge_slice].reshape(pad_singleton)
+    return np.concatenate((edge_arr.repeat(pad_amt, axis=axis), arr),
+                          axis=axis)
+
+
+def _append_edge(arr, pad_amt, axis=-1):
+    """
+    Append `pad_amt` to `arr` along `axis` by extending edge values.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.
-    kwargs : keyword arguments
-        Keyword arguments. Only 'stat_length' is used. 'stat_length'
-        defaults to the entire vector if not supplied.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to append.
+    axis : int
+        Axis along which to pad `arr`.
 
-    Return
-    ------
-    _maximum : ndarray
-        Padded vector
+    Returns
+    -------
+    padarr : ndarray
+        Output array, extended by `pad_amt` edge values prepended along
+        `axis`.
 
     """
-    sbvec, savec = _create_stat_vectors(vector, pad_tuple, iaxis, kwargs)
-    return _create_vector(vector, pad_tuple, max(sbvec), max(savec))
+    if pad_amt == 0:
+        return arr
 
+    edge_slice = tuple(slice(None) if i != axis else arr.shape[axis] - 1
+                       for (i, x) in enumerate(arr.shape))
 
-def _minimum(vector, pad_tuple, iaxis, kwargs):
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+    edge_arr = arr[edge_slice].reshape(pad_singleton)
+    return np.concatenate((arr, edge_arr.repeat(pad_amt, axis=axis)),
+                          axis=axis)
+
+
+def _prepend_ramp(arr, pad_amt, end, axis=-1):
     """
-    Private function to calculate the before/after vectors for pad_minimum.
+    Prepend linear ramp along `axis`.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.
-    kwargs : keyword arguments
-        Keyword arguments. Only 'stat_length' is used. 'stat_length'
-        defaults to the entire vector if not supplied.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to prepend.
+    end : scalar
+        Constal value to use. For best results should be of type `arr.dtype`;
+        if not `arr.dtype` will be cast to `arr.dtype`.
+    axis : int
+        Axis along which to pad `arr`.
 
-    Return
-    ------
-    _minimum : ndarray
-        Padded vector
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values prepended along `axis`. The
+        prepended region ramps linearly from the edge value to `end`.
 
     """
-    sbvec, savec = _create_stat_vectors(vector, pad_tuple, iaxis, kwargs)
-    return _create_vector(vector, pad_tuple, min(sbvec), min(savec))
+    if pad_amt == 0:
+        return arr
+
+    # Generate shape for final concatenated array
+    padshape = tuple(x if i != axis else pad_amt
+                     for (i, x) in enumerate(arr.shape))
 
+    # Generate an n-dimensional array incrementing along `axis`
+    ramp_arr = _arange_ndarray(arr, padshape, axis,
+                               reverse=True).astype(np.float64)
 
-def _median(vector, pad_tuple, iaxis, kwargs):
+
+    # Appropriate slicing to extract n-dimensional edge along `axis`
+    edge_slice = tuple(slice(None) if i != axis else 0
+                       for (i, x) in enumerate(arr.shape))
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+
+    # Extract edge, reshape to original rank, and extend along `axis`
+    edge_pad = arr[edge_slice].reshape(pad_singleton).repeat(pad_amt, axis)
+
+    # Linear ramp
+    slope = (end - edge_pad) / float(pad_amt)
+    ramp_arr = ramp_arr * slope
+    ramp_arr += edge_pad
+    _round_ifneeded(ramp_arr, arr.dtype)
+
+    # Ramp values will most likely be float, cast them to the same type as arr
+    return np.concatenate((ramp_arr.astype(arr.dtype), arr), axis=axis)
+
+
+def _append_ramp(arr, pad_amt, end, axis=-1):
     """
-    Private function to calculate the before/after vectors for pad_median.
+    Append linear ramp along `axis`.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.
-    kwargs : keyword arguments
-        Keyword arguments. Only 'stat_length' is used. 'stat_length'
-        defaults to the entire vector if not supplied.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to append.
+    end : scalar
+        Constal value to use. For best results should be of type `arr.dtype`;
+        if not `arr.dtype` will be cast to `arr.dtype`.
+    axis : int
+        Axis along which to pad `arr`.
 
-    Return
-    ------
-    _median : ndarray
-        Padded vector
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values appended along `axis`. The
+        appended region ramps linearly from the edge value to `end`.
 
     """
-    sbvec, savec = _create_stat_vectors(vector, pad_tuple, iaxis, kwargs)
-    return _create_vector(vector, pad_tuple, np.median(sbvec),
-                          np.median(savec))
+    if pad_amt == 0:
+        return arr
+
+    # Generate shape for final concatenated array
+    padshape = tuple(x if i != axis else pad_amt
+                     for (i, x) in enumerate(arr.shape))
+
+    # Generate an n-dimensional array incrementing along `axis`
+    ramp_arr = _arange_ndarray(arr, padshape, axis,
+                               reverse=False).astype(np.float64)
 
+    # Slice a chunk from the edge to calculate stats on
+    edge_slice = tuple(slice(None) if i != axis else -1
+                       for (i, x) in enumerate(arr.shape))
 
-def _mean(vector, pad_tuple, iaxis, kwargs):
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+
+    # Extract edge, reshape to original rank, and extend along `axis`
+    edge_pad = arr[edge_slice].reshape(pad_singleton).repeat(pad_amt, axis)
+
+    # Linear ramp
+    slope = (end - edge_pad) / float(pad_amt)
+    ramp_arr = ramp_arr * slope
+    ramp_arr += edge_pad
+    _round_ifneeded(ramp_arr, arr.dtype)
+
+    # Ramp values will most likely be float, cast them to the same type as arr
+    return np.concatenate((arr, ramp_arr.astype(arr.dtype)), axis=axis)
+
+
+def _prepend_max(arr, pad_amt, num, axis=-1):
     """
-    Private function to calculate the before/after vectors for pad_mean.
+    Prepend `pad_amt` maximum values along `axis`.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.
-    kwargs : keyword arguments
-        Keyword arguments. Only 'stat_length' is used. 'stat_length'
-        defaults to the entire vector if not supplied.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to prepend.
+    num : int
+        Depth into `arr` along `axis` to calculate maximum.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
 
-    Return
-    ------
-    _mean : ndarray
-        Padded vector
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values appended along `axis`. The
+        prepended region is the maximum of the first `num` values along
+        `axis`.
 
     """
-    sbvec, savec = _create_stat_vectors(vector, pad_tuple, iaxis, kwargs)
-    return _create_vector(vector, pad_tuple, np.average(sbvec),
-                          np.average(savec))
+    if pad_amt == 0:
+        return arr
+
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _prepend_edge(arr, pad_amt, axis)
+
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
+
+    # Slice a chunk from the edge to calculate stats on
+    max_slice = tuple(slice(None) if i != axis else slice(num)
+                      for (i, x) in enumerate(arr.shape))
 
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
 
-def _constant(vector, pad_tuple, iaxis, kwargs):
+    # Extract slice, calculate max, reshape to add singleton dimension back
+    max_chunk = arr[max_slice].max(axis=axis).reshape(pad_singleton)
+
+    # Concatenate `arr` with `max_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate((max_chunk.repeat(pad_amt, axis=axis), arr),
+                          axis=axis)
+
+
+def _append_max(arr, pad_amt, num, axis=-1):
     """
-    Private function to calculate the before/after vectors for
-    pad_constant.
+    Pad one `axis` of `arr` with the maximum of the last `num` elements.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.
-    kwargs : keyword arguments
-        Keyword arguments. Need 'constant_values' keyword argument.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to append.
+    num : int
+        Depth into `arr` along `axis` to calculate maximum.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
 
-    Return
-    ------
-    _constant : ndarray
-        Padded vector
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values appended along `axis`. The
+        appended region is the maximum of the final `num` values along `axis`.
 
     """
-    nconstant = kwargs['constant_values'][iaxis]
-    return _create_vector(vector, pad_tuple, nconstant[0], nconstant[1])
+    if pad_amt == 0:
+        return arr
+
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _append_edge(arr, pad_amt, axis)
+
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
+
+    # Slice a chunk from the edge to calculate stats on
+    end = arr.shape[axis] - 1
+    if num is not None:
+        max_slice = tuple(
+            slice(None) if i != axis else slice(end, end - num, -1)
+            for (i, x) in enumerate(arr.shape))
+    else:
+        max_slice = tuple(slice(None) for x in arr.shape)
 
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
 
-def _linear_ramp(vector, pad_tuple, iaxis, kwargs):
+    # Extract slice, calculate max, reshape to add singleton dimension back
+    max_chunk = arr[max_slice].max(axis=axis).reshape(pad_singleton)
+
+    # Concatenate `arr` with `max_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate((arr, max_chunk.repeat(pad_amt, axis=axis)),
+                          axis=axis)
+
+
+def _prepend_mean(arr, pad_amt, num, axis=-1):
     """
-    Private function to calculate the before/after vectors for
-    pad_linear_ramp.
+    Prepend `pad_amt` mean values along `axis`.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across. Not used in _linear_ramp.
-    kwargs : keyword arguments
-        Keyword arguments. Not used in _linear_ramp.
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to prepend.
+    num : int
+        Depth into `arr` along `axis` to calculate mean.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
 
-    Return
-    ------
-    _linear_ramp : ndarray
-        Padded vector
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values prepended along `axis`. The
+        prepended region is the mean of the first `num` values along `axis`.
 
     """
-    end_values = kwargs['end_values'][iaxis]
-    before_delta = ((vector[pad_tuple[0]] - end_values[0])
-                    / float(pad_tuple[0]))
-    after_delta = ((vector[-pad_tuple[1] - 1] - end_values[1])
-                   / float(pad_tuple[1]))
+    if pad_amt == 0:
+        return arr
 
-    before_vector = np.ones((pad_tuple[0], )) * end_values[0]
-    before_vector = before_vector.astype(vector.dtype)
-    for i in range(len(before_vector)):
-        before_vector[i] = before_vector[i] + i * before_delta
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _prepend_edge(arr, pad_amt, axis)
 
-    after_vector = np.ones((pad_tuple[1], )) * end_values[1]
-    after_vector = after_vector.astype(vector.dtype)
-    for i in range(len(after_vector)):
-        after_vector[i] = after_vector[i] + i * after_delta
-    after_vector = after_vector[::-1]
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
 
-    return _create_vector(vector, pad_tuple, before_vector, after_vector)
+    # Slice a chunk from the edge to calculate stats on
+    mean_slice = tuple(slice(None) if i != axis else slice(num)
+                       for (i, x) in enumerate(arr.shape))
 
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
 
-def _reflect(vector, pad_tuple, iaxis, kwargs):
+    # Extract slice, calculate mean, reshape to add singleton dimension back
+    mean_chunk = arr[mean_slice].mean(axis).reshape(pad_singleton)
+    _round_ifneeded(mean_chunk, arr.dtype)
+
+    # Concatenate `arr` with `mean_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate((mean_chunk.repeat(pad_amt, axis).astype(arr.dtype),
+                           arr), axis=axis)
+
+
+def _append_mean(arr, pad_amt, num, axis=-1):
     """
-    Private function to calculate the before/after vectors for pad_reflect.
+    Append `pad_amt` mean values along `axis`.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across. Not used in _reflect.
-    kwargs : keyword arguments
-        Keyword arguments. Not used in _reflect.
-
-    Return
-    ------
-    _reflect : ndarray
-        Padded vector
-
-    """
-    # Can't have pad_tuple[1] be used in the slice if == to 0.
-    if pad_tuple[1] == 0:
-        after_vector = vector[pad_tuple[0]:None]
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to append.
+    num : int
+        Depth into `arr` along `axis` to calculate mean.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values appended along `axis`. The
+        appended region is the maximum of the final `num` values along `axis`.
+
+    """
+    if pad_amt == 0:
+        return arr
+
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _append_edge(arr, pad_amt, axis)
+
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
+
+    # Slice a chunk from the edge to calculate stats on
+    end = arr.shape[axis] - 1
+    if num is not None:
+        mean_slice = tuple(
+            slice(None) if i != axis else slice(end, end - num, -1)
+            for (i, x) in enumerate(arr.shape))
     else:
-        after_vector = vector[pad_tuple[0]:-pad_tuple[1]]
+        mean_slice = tuple(slice(None) for x in arr.shape)
 
-    reverse = after_vector[::-1]
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
 
-    before_vector = np.resize(
-        np.concatenate((after_vector[1:-1], reverse)), pad_tuple[0])[::-1]
-    after_vector = np.resize(
-        np.concatenate((reverse[1:-1], after_vector)), pad_tuple[1])
+    # Extract slice, calculate mean, reshape to add singleton dimension back
+    mean_chunk = arr[mean_slice].mean(axis=axis).reshape(pad_singleton)
+    _round_ifneeded(mean_chunk, arr.dtype)
 
-    if kwargs['reflect_type'] == 'even':
-        pass
-    elif kwargs['reflect_type'] == 'odd':
-        before_vector = 2 * vector[pad_tuple[0]] - before_vector
-        after_vector = 2 * vector[-pad_tuple[-1] - 1] - after_vector
-    else:
-        raise ValueError("The keyword '%s' cannot have the value '%s'."
-                         % ('reflect_type', kwargs['reflect_type']))
-    return _create_vector(vector, pad_tuple, before_vector, after_vector)
+    # Concatenate `arr` with `mean_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate(
+        (arr, mean_chunk.repeat(pad_amt, axis).astype(arr.dtype)), axis=axis)
 
 
-def _symmetric(vector, pad_tuple, iaxis, kwargs):
+def _prepend_med(arr, pad_amt, num, axis=-1):
     """
-    Private function to calculate the before/after vectors for
-    pad_symmetric.
+    Prepend `pad_amt` median values along `axis`.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across. Not used in _symmetric.
-    kwargs : keyword arguments
-        Keyword arguments. Not used in _symmetric.
-
-    Return
-    ------
-    _symmetric : ndarray
-        Padded vector
-
-    """
-    if pad_tuple[1] == 0:
-        after_vector = vector[pad_tuple[0]:None]
-    else:
-        after_vector = vector[pad_tuple[0]:-pad_tuple[1]]
-
-    before_vector = np.resize(
-        np.concatenate((after_vector, after_vector[::-1])),
-        pad_tuple[0])[::-1]
-    after_vector = np.resize(
-        np.concatenate((after_vector[::-1], after_vector)),
-        pad_tuple[1])
-
-    if kwargs['reflect_type'] == 'even':
-        pass
-    elif kwargs['reflect_type'] == 'odd':
-        before_vector = 2 * vector[pad_tuple[0]] - before_vector
-        after_vector = 2 * vector[-pad_tuple[1] - 1] - after_vector
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to prepend.
+    num : int
+        Depth into `arr` along `axis` to calculate median.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values prepended along `axis`. The
+        prepended region is the median of the first `num` values along `axis`.
+
+    """
+    if pad_amt == 0:
+        return arr
+
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _prepend_edge(arr, pad_amt, axis)
+
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
+
+    # Slice a chunk from the edge to calculate stats on
+    med_slice = tuple(slice(None) if i != axis else slice(num)
+                      for (i, x) in enumerate(arr.shape))
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+
+    # Extract slice, calculate median, reshape to add singleton dimension back
+    med_chunk = np.median(arr[med_slice], axis=axis).reshape(pad_singleton)
+    _round_ifneeded(med_chunk, arr.dtype)
+
+    # Concatenate `arr` with `med_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate(
+        (med_chunk.repeat(pad_amt, axis).astype(arr.dtype), arr), axis=axis)
+
+
+def _append_med(arr, pad_amt, num, axis=-1):
+    """
+    Append `pad_amt` median values along `axis`.
+
+    Parameters
+    ----------
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to append.
+    num : int
+        Depth into `arr` along `axis` to calculate median.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values appended along `axis`. The
+        appended region is the median of the final `num` values along `axis`.
+
+    """
+    if pad_amt == 0:
+        return arr
+
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _append_edge(arr, pad_amt, axis)
+
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
+
+    # Slice a chunk from the edge to calculate stats on
+    end = arr.shape[axis] - 1
+    if num is not None:
+        med_slice = tuple(
+            slice(None) if i != axis else slice(end, end - num, -1)
+            for (i, x) in enumerate(arr.shape))
     else:
-        raise ValueError("The keyword '%s' cannot have the value '%s'."
-                         % ('reflect_type', kwargs['reflect_type']))
-    return _create_vector(vector, pad_tuple, before_vector, after_vector)
+        med_slice = tuple(slice(None) for x in arr.shape)
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+
+    # Extract slice, calculate median, reshape to add singleton dimension back
+    med_chunk = np.median(arr[med_slice], axis=axis).reshape(pad_singleton)
+    _round_ifneeded(med_chunk, arr.dtype)
+
+    # Concatenate `arr` with `med_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate(
+        (arr, med_chunk.repeat(pad_amt, axis).astype(arr.dtype)), axis=axis)
+
+
+def _prepend_min(arr, pad_amt, num, axis=-1):
+    """
+    Prepend `pad_amt` minimum values along `axis`.
+
+    Parameters
+    ----------
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to prepend.
+    num : int
+        Depth into `arr` along `axis` to calculate minimum.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values prepended along `axis`. The
+        prepended region is the minimum of the first `num` values along
+        `axis`.
+
+    """
+    if pad_amt == 0:
+        return arr
+
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _prepend_edge(arr, pad_amt, axis)
+
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
+
+    # Slice a chunk from the edge to calculate stats on
+    min_slice = tuple(slice(None) if i != axis else slice(num)
+                      for (i, x) in enumerate(arr.shape))
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+
+    # Extract slice, calculate min, reshape to add singleton dimension back
+    min_chunk = arr[min_slice].min(axis=axis).reshape(pad_singleton)
 
+    # Concatenate `arr` with `min_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate((min_chunk.repeat(pad_amt, axis=axis), arr),
+                          axis=axis)
 
-def _wrap(vector, pad_tuple, iaxis, kwargs):
+
+def _append_min(arr, pad_amt, num, axis=-1):
     """
-    Private function to calculate the before/after vectors for pad_wrap.
+    Append `pad_amt` median values along `axis`.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.  Not used in _wrap.
-    kwargs : keyword arguments
-        Keyword arguments. Not used in _wrap.
-
-    Return
-    ------
-    _wrap : ndarray
-        Padded vector
-
-    """
-    if pad_tuple[1] == 0:
-        after_vector = vector[pad_tuple[0]:None]
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : int
+        Amount of padding to append.
+    num : int
+        Depth into `arr` along `axis` to calculate minimum.
+        Range: [1, `arr.shape[axis]`] or None (entire axis)
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt` values appended along `axis`. The
+        appended region is the minimum of the final `num` values along `axis`.
+
+    """
+    if pad_amt == 0:
+        return arr
+
+    # Equivalent to edge padding for single value, so do that instead
+    if num == 1:
+        return _append_edge(arr, pad_amt, axis)
+
+    # Use entire array if `num` is too large
+    if num is not None:
+        if num >= arr.shape[axis]:
+            num = None
+
+    # Slice a chunk from the edge to calculate stats on
+    end = arr.shape[axis] - 1
+    if num is not None:
+        min_slice = tuple(
+            slice(None) if i != axis else slice(end, end - num, -1)
+            for (i, x) in enumerate(arr.shape))
     else:
-        after_vector = vector[pad_tuple[0]:-pad_tuple[1]]
+        min_slice = tuple(slice(None) for x in arr.shape)
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+
+    # Extract slice, calculate min, reshape to add singleton dimension back
+    min_chunk = arr[min_slice].min(axis=axis).reshape(pad_singleton)
+
+    # Concatenate `arr` with `min_chunk`, extended along `axis` by `pad_amt`
+    return np.concatenate((arr, min_chunk.repeat(pad_amt, axis=axis)),
+                          axis=axis)
+
+
+def _pad_ref(arr, pad_amt, method, axis=-1):
+    """
+    Pad `axis` of `arr` by reflection.
+
+    Parameters
+    ----------
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : tuple of ints, length 2
+        Padding to (prepend, append) along `axis`.
+    method : str
+        Controls method of reflection; options are 'even' or 'odd'.
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt[0]` values prepended and `pad_amt[1]`
+        values appended along `axis`. Both regions are padded with reflected
+        values from the original array.
+
+    Notes
+    -----
+    This algorithm does not pad with repetition, i.e. the edges are not
+    repeated in the reflection. For that behavior, use `method='symmetric'`.
+
+    The modes 'reflect', 'symmetric', and 'wrap' must be padded with a
+    single function, lest the indexing tricks in non-integer multiples of the
+    original shape would violate repetition in the final iteration.
+
+    """
+    # Implicit booleanness to test for zero (or None) in any scalar type
+    if pad_amt[0] == 0 and pad_amt[1] == 0:
+        return arr
+
+    ##########################################################################
+    # Prepended region
+
+    # Slice off a reverse indexed chunk from near edge to pad `arr` before
+    ref_slice = tuple(slice(None) if i != axis else slice(pad_amt[0], 0, -1)
+                      for (i, x) in enumerate(arr.shape))
+
+    ref_chunk1 = arr[ref_slice]
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+    if pad_amt[0] == 1:
+        ref_chunk1 = ref_chunk1.reshape(pad_singleton)
+
+    # Memory/computationally more expensive, only do this if `method='odd'`
+    if 'odd' in method and pad_amt[0] > 0:
+        edge_slice1 = tuple(slice(None) if i != axis else 0
+                            for (i, x) in enumerate(arr.shape))
+        edge_chunk = arr[edge_slice1].reshape(pad_singleton)
+        ref_chunk1 = 2 * edge_chunk - ref_chunk1
+        del edge_chunk
+
+    ##########################################################################
+    # Appended region
+
+    # Slice off a reverse indexed chunk from far edge to pad `arr` after
+    start = arr.shape[axis] - pad_amt[1] - 1
+    end = arr.shape[axis] - 1
+    ref_slice = tuple(slice(None) if i != axis else slice(start, end)
+                      for (i, x) in enumerate(arr.shape))
+    rev_idx = tuple(slice(None) if i != axis else slice(None, None, -1)
+                    for (i, x) in enumerate(arr.shape))
+    ref_chunk2 = arr[ref_slice][rev_idx]
+
+    if pad_amt[1] == 1:
+        ref_chunk2 = ref_chunk2.reshape(pad_singleton)
+
+    if 'odd' in method:
+        edge_slice2 = tuple(slice(None) if i != axis else -1
+                            for (i, x) in enumerate(arr.shape))
+        edge_chunk = arr[edge_slice2].reshape(pad_singleton)
+        ref_chunk2 = 2 * edge_chunk - ref_chunk2
+        del edge_chunk
+
+    # Concatenate `arr` with both chunks, extending along `axis`
+    return np.concatenate((ref_chunk1, arr, ref_chunk2), axis=axis)
+
+
+def _pad_sym(arr, pad_amt, method, axis=-1):
+    """
+    Pad `axis` of `arr` by symmetry.
+
+    Parameters
+    ----------
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : tuple of ints, length 2
+        Padding to (prepend, append) along `axis`.
+    method : str
+        Controls method of symmetry; options are 'even' or 'odd'.
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt[0]` values prepended and `pad_amt[1]`
+        values appended along `axis`. Both regions are padded with symmetric
+        values from the original array.
+
+    Notes
+    -----
+    This algorithm DOES pad with repetition, i.e. the edges are repeated.
+    For a method that does not repeat edges, use `method='reflect'`.
+
+    The modes 'reflect', 'symmetric', and 'wrap' must be padded with a
+    single function, lest the indexing tricks in non-integer multiples of the
+    original shape would violate repetition in the final iteration.
+
+    """
+    # Implicit booleanness to test for zero (or None) in any scalar type
+    if pad_amt[0] == 0 and pad_amt[1] == 0:
+        return arr
+
+    ##########################################################################
+    # Prepended region
+
+    # Slice off a reverse indexed chunk from near edge to pad `arr` before
+    sym_slice = tuple(slice(None) if i != axis else slice(0, pad_amt[0])
+                      for (i, x) in enumerate(arr.shape))
+    rev_idx = tuple(slice(None) if i != axis else slice(None, None, -1)
+                    for (i, x) in enumerate(arr.shape))
+    sym_chunk1 = arr[sym_slice][rev_idx]
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+    if pad_amt[0] == 1:
+        sym_chunk1 = sym_chunk1.reshape(pad_singleton)
+
+    # Memory/computationally more expensive, only do this if `method='odd'`
+    if 'odd' in method and pad_amt[0] > 0:
+        edge_slice1 = tuple(slice(None) if i != axis else 0
+                            for (i, x) in enumerate(arr.shape))
+        edge_chunk = arr[edge_slice1].reshape(pad_singleton)
+        sym_chunk1 = 2 * edge_chunk - sym_chunk1
+        del edge_chunk
+
+    ##########################################################################
+    # Appended region
+
+    # Slice off a reverse indexed chunk from far edge to pad `arr` after
+    start = arr.shape[axis] - pad_amt[1]
+    end = arr.shape[axis]
+    sym_slice = tuple(slice(None) if i != axis else slice(start, end)
+                      for (i, x) in enumerate(arr.shape))
+    sym_chunk2 = arr[sym_slice][rev_idx]
+
+    if pad_amt[1] == 1:
+        sym_chunk2 = sym_chunk2.reshape(pad_singleton)
+
+    if 'odd' in method:
+        edge_slice2 = tuple(slice(None) if i != axis else -1
+                            for (i, x) in enumerate(arr.shape))
+        edge_chunk = arr[edge_slice2].reshape(pad_singleton)
+        sym_chunk2 = 2 * edge_chunk - sym_chunk2
+        del edge_chunk
+
+    # Concatenate `arr` with both chunks, extending along `axis`
+    return np.concatenate((sym_chunk1, arr, sym_chunk2), axis=axis)
+
+
+def _pad_wrap(arr, pad_amt, axis=-1):
+    """
+    Pad `axis` of `arr` via wrapping.
+
+    Parameters
+    ----------
+    arr : ndarray
+        Input array of arbitrary shape.
+    pad_amt : tuple of ints, length 2
+        Padding to (prepend, append) along `axis`.
+    axis : int
+        Axis along which to pad `arr`.
+
+    Returns
+    -------
+    padarr : ndarray
+        Output array, with `pad_amt[0]` values prepended and `pad_amt[1]`
+        values appended along `axis`. Both regions are padded wrapped values
+        from the opposite end of `axis`.
+
+    Notes
+    -----
+    This method of padding is also known as 'tile' or 'tiling'.
+
+    The modes 'reflect', 'symmetric', and 'wrap' must be padded with a
+    single function, lest the indexing tricks in non-integer multiples of the
+    original shape would violate repetition in the final iteration.
+
+    """
+    # Implicit booleanness to test for zero (or None) in any scalar type
+    if pad_amt[0] == 0 and pad_amt[1] == 0:
+        return arr
+
+    ##########################################################################
+    # Prepended region
+
+    # Slice off a reverse indexed chunk from near edge to pad `arr` before
+    start = arr.shape[axis] - pad_amt[0]
+    end = arr.shape[axis]
+    wrap_slice = tuple(slice(None) if i != axis else slice(start, end)
+                       for (i, x) in enumerate(arr.shape))
+    wrap_chunk1 = arr[wrap_slice]
+
+    # Shape to restore singleton dimension after slicing
+    pad_singleton = tuple(x if i != axis else 1
+                          for (i, x) in enumerate(arr.shape))
+    if pad_amt[0] == 1:
+        wrap_chunk1 = wrap_chunk1.reshape(pad_singleton)
+
+    ##########################################################################
+    # Appended region
+
+    # Slice off a reverse indexed chunk from far edge to pad `arr` after
+    wrap_slice = tuple(slice(None) if i != axis else slice(0, pad_amt[1])
+                       for (i, x) in enumerate(arr.shape))
+    wrap_chunk2 = arr[wrap_slice]
+
+    if pad_amt[1] == 1:
+        wrap_chunk2 = wrap_chunk2.reshape(pad_singleton)
+
+    # Concatenate `arr` with both chunks, extending along `axis`
+    return np.concatenate((wrap_chunk1, arr, wrap_chunk2), axis=axis)
+
+
+def _normalize_shape(narray, shape):
+    """
+    Private function which does some checks and normalizes the possibly
+    much simpler representations of 'pad_width', 'stat_length',
+    'constant_values', 'end_values'.
 
-    before_vector = np.resize(after_vector[::-1], pad_tuple[0])[::-1]
-    after_vector = np.resize(after_vector, pad_tuple[1])
+    Parameters
+    ----------
+    narray : ndarray
+        Input ndarray
+    shape : {sequence, int}, optional
+        The width of padding (pad_width) or the number of elements on the
+        edge of the narray used for statistics (stat_length).
+        ((before_1, after_1), ... (before_N, after_N)) unique number of
+        elements for each axis where `N` is rank of `narray`.
+        ((before, after),) yields same before and after constants for each
+        axis.
+        (constant,) or int is a shortcut for before = after = constant for
+        all axes.
 
-    return _create_vector(vector, pad_tuple, before_vector, after_vector)
+    Returns
+    -------
+    _normalize_shape : tuple of tuples
+        int                               => ((int, int), (int, int), ...)
+        [[int1, int2], [int3, int4], ...] => ((int1, int2), (int3, int4), ...)
+        ((int1, int2), (int3, int4), ...) => no change
+        [[int1, int2], ]                  => ((int1, int2), (int1, int2), ...)
+        ((int1, int2), )                  => ((int1, int2), (int1, int2), ...)
+        [[int ,     ], ]                  => ((int, int), (int, int), ...)
+        ((int ,     ), )                  => ((int, int), (int, int), ...)
+
+    """
+    normshp = None
+    shapelen = len(np.shape(narray))
+    if (isinstance(shape, int)) or shape is None:
+        normshp = ((shape, shape), ) * shapelen
+    elif (isinstance(shape, (tuple, list))
+            and isinstance(shape[0], (tuple, list))
+            and len(shape) == shapelen):
+        normshp = shape
+        for i in normshp:
+            if len(i) != 2:
+                fmt = "Unable to create correctly shaped tuple from %s"
+                raise ValueError(fmt % (normshp,))
+    elif (isinstance(shape, (tuple, list))
+            and isinstance(shape[0], (int, float, long))
+            and len(shape) == 1):
+        normshp = ((shape[0], shape[0]), ) * shapelen
+    elif (isinstance(shape, (tuple, list))
+            and isinstance(shape[0], (int, float, long))
+            and len(shape) == 2):
+        normshp = (shape, ) * shapelen
+    if normshp is None:
+        fmt = "Unable to create correctly shaped tuple from %s"
+        raise ValueError(fmt % (shape,))
+    return normshp
 
 
-def _edge(vector, pad_tuple, iaxis, kwargs):
+def _validate_lengths(narray, number_elements):
     """
-    Private function to calculate the before/after vectors for pad_edge.
+    Private function which does some checks and reformats pad_width and
+    stat_length using _normalize_shape.
 
     Parameters
     ----------
-    vector : ndarray
-        Input vector that already includes empty padded values.
-    pad_tuple : tuple
-        This tuple represents the (before, after) width of the padding
-        along this particular iaxis.
-    iaxis : int
-        The axis currently being looped across.  Not used in _edge.
-    kwargs : keyword arguments
-        Keyword arguments. Not used in _edge.
+    narray : ndarray
+        Input ndarray
+    number_elements : {sequence, int}, optional
+        The width of padding (pad_width) or the number of elements on the edge
+        of the narray used for statistics (stat_length).
+        ((before_1, after_1), ... (before_N, after_N)) unique number of
+        elements for each axis.
+        ((before, after),) yields same before and after constants for each
+        axis.
+        (constant,) or int is a shortcut for before = after = constant for all
+        axes.
 
-    Return
-    ------
-    _edge : ndarray
-        Padded vector
+    Returns
+    -------
+    _validate_lengths : tuple of tuples
+        int                               => ((int, int), (int, int), ...)
+        [[int1, int2], [int3, int4], ...] => ((int1, int2), (int3, int4), ...)
+        ((int1, int2), (int3, int4), ...) => no change
+        [[int1, int2], ]                  => ((int1, int2), (int1, int2), ...)
+        ((int1, int2), )                  => ((int1, int2), (int1, int2), ...)
+        [[int ,     ], ]                  => ((int, int), (int, int), ...)
+        ((int ,     ), )                  => ((int, int), (int, int), ...)
 
     """
-    return _create_vector(vector, pad_tuple, vector[pad_tuple[0]],
-                          vector[-pad_tuple[1] - 1])
+    normshp = _normalize_shape(narray, number_elements)
+    for i in normshp:
+        chk = [1 if x is None else x for x in i]
+        chk = [1 if x > 0 else -1 for x in chk]
+        if (chk[0] < 0) or (chk[1] < 0):
+            fmt = "%s cannot contain negative values."
+            raise ValueError(fmt % (number_elements,))
+    return normshp
 
 
 ###############################################################################
@@ -714,19 +1281,6 @@ def pad(array, pad_width, mode=None, **kwargs):
     narray = np.array(array)
     pad_width = _validate_lengths(narray, pad_width)
 
-    modefunc = {
-        'constant': _constant,
-        'edge': _edge,
-        'linear_ramp': _linear_ramp,
-        'maximum': _maximum,
-        'mean': _mean,
-        'median': _median,
-        'minimum': _minimum,
-        'reflect': _reflect,
-        'symmetric': _symmetric,
-        'wrap': _wrap,
-        }
-
     allowedkwargs = {
         'constant': ['constant_values'],
         'edge': [],
@@ -748,8 +1302,6 @@ def pad(array, pad_width, mode=None, **kwargs):
         }
 
     if isinstance(mode, str):
-        function = modefunc[mode]
-
         # Make sure have allowed kwargs appropriate for mode
         for key in kwargs:
             if key not in allowedkwargs[mode]:
@@ -762,35 +1314,156 @@ def pad(array, pad_width, mode=None, **kwargs):
 
         # Need to only normalize particular keywords.
         for i in kwargs:
-            if i == 'stat_length' and kwargs[i]:
+            if i == 'stat_length':
                 kwargs[i] = _validate_lengths(narray, kwargs[i])
             if i in ['end_values', 'constant_values']:
                 kwargs[i] = _normalize_shape(narray, kwargs[i])
     elif mode is None:
         raise ValueError('Keyword "mode" must be a function or one of %s.' %
-                         (list(modefunc.keys()),))
+                         (list(allowedkwargs.keys()),))
     else:
-        # User supplied function, I hope
+        # Drop back to old, slower np.apply_along_axis mode for user-supplied
+        # vector function
         function = mode
 
-    # Create a new padded array
-    rank = list(range(len(narray.shape)))
-    total_dim_increase = [np.sum(pad_width[i]) for i in rank]
-    offset_slices = [slice(pad_width[i][0],
-                           pad_width[i][0] + narray.shape[i])
-                     for i in rank]
-    new_shape = np.array(narray.shape) + total_dim_increase
-    newmat = np.zeros(new_shape).astype(narray.dtype)
-
-    # Insert the original array into the padded array
-    newmat[offset_slices] = narray
-
-    # This is the core of pad ...
-    for iaxis in rank:
-        np.apply_along_axis(function,
-                            iaxis,
-                            newmat,
-                            pad_width[iaxis],
-                            iaxis,
-                            kwargs)
+        # Create a new padded array
+        rank = list(range(len(narray.shape)))
+        total_dim_increase = [np.sum(pad_width[i]) for i in rank]
+        offset_slices = [slice(pad_width[i][0],
+                               pad_width[i][0] + narray.shape[i])
+                         for i in rank]
+        new_shape = np.array(narray.shape) + total_dim_increase
+        newmat = np.zeros(new_shape).astype(narray.dtype)
+
+        # Insert the original array into the padded array
+        newmat[offset_slices] = narray
+
+        # This is the core of pad ...
+        for iaxis in rank:
+            np.apply_along_axis(function,
+                                iaxis,
+                                newmat,
+                                pad_width[iaxis],
+                                iaxis,
+                                kwargs)
+        return newmat
+
+    # If we get here, use new padding method
+    newmat = narray.copy()
+
+    # API preserved, but completely new algorithm which pads by building the
+    # entire block to pad before/after `arr` with in one step, for each axis.
+    if mode == 'constant':
+        for axis, ((pad_before, pad_after), (before_val, after_val)) \
+                in enumerate(zip(pad_width, kwargs['constant_values'])):
+            newmat = _prepend_const(newmat, pad_before, before_val, axis)
+            newmat = _append_const(newmat, pad_after, after_val, axis)
+
+    elif mode == 'edge':
+        for axis, (pad_before, pad_after) in enumerate(pad_width):
+            newmat = _prepend_edge(newmat, pad_before, axis)
+            newmat = _append_edge(newmat, pad_after, axis)
+
+    elif mode == 'linear_ramp':
+        for axis, ((pad_before, pad_after), (before_val, after_val)) \
+                in enumerate(zip(pad_width, kwargs['end_values'])):
+            newmat = _prepend_ramp(newmat, pad_before, before_val, axis)
+            newmat = _append_ramp(newmat, pad_after, after_val, axis)
+
+    elif mode == 'maximum':
+        for axis, ((pad_before, pad_after), (chunk_before, chunk_after)) \
+                in enumerate(zip(pad_width, kwargs['stat_length'])):
+            newmat = _prepend_max(newmat, pad_before, chunk_before, axis)
+            newmat = _append_max(newmat, pad_after, chunk_after, axis)
+
+    elif mode == 'mean':
+        for axis, ((pad_before, pad_after), (chunk_before, chunk_after)) \
+                in enumerate(zip(pad_width, kwargs['stat_length'])):
+            newmat = _prepend_mean(newmat, pad_before, chunk_before, axis)
+            newmat = _append_mean(newmat, pad_after, chunk_after, axis)
+
+    elif mode == 'median':
+        for axis, ((pad_before, pad_after), (chunk_before, chunk_after)) \
+                in enumerate(zip(pad_width, kwargs['stat_length'])):
+            newmat = _prepend_med(newmat, pad_before, chunk_before, axis)
+            newmat = _append_med(newmat, pad_after, chunk_after, axis)
+
+    elif mode == 'minimum':
+        for axis, ((pad_before, pad_after), (chunk_before, chunk_after)) \
+                in enumerate(zip(pad_width, kwargs['stat_length'])):
+            newmat = _prepend_min(newmat, pad_before, chunk_before, axis)
+            newmat = _append_min(newmat, pad_after, chunk_after, axis)
+
+    elif mode == 'reflect':
+        for axis, (pad_before, pad_after) in enumerate(pad_width):
+            # Recursive padding along any axis where `pad_amt` is too large
+            # for indexing tricks. We can only safely pad the original axis
+            # length, to keep the period of the reflections consistent.
+            if ((pad_before > 0) or
+                    (pad_after > 0)) and newmat.shape[axis] == 1:
+                # Extending singleton dimension for 'reflect' is legacy
+                # behavior; it really should raise an error.
+                newmat = _prepend_edge(newmat, pad_before, axis)
+                newmat = _append_edge(newmat, pad_after, axis)
+                continue
+
+            method = kwargs['reflect_type']
+            safe_pad = newmat.shape[axis] - 1
+            repeat = safe_pad
+            while ((pad_before > safe_pad) or (pad_after > safe_pad)):
+                offset = 0
+                pad_iter_b = min(safe_pad,
+                                 safe_pad * (pad_before // safe_pad))
+                pad_iter_a = min(safe_pad, safe_pad * (pad_after // safe_pad))
+                newmat = _pad_ref(newmat, (pad_iter_b,
+                                           pad_iter_a), method, axis)
+                pad_before -= pad_iter_b
+                pad_after -= pad_iter_a
+                if pad_iter_b > 0:
+                    offset += 1
+                if pad_iter_a > 0:
+                    offset += 1
+                safe_pad += pad_iter_b + pad_iter_a
+            newmat = _pad_ref(newmat, (pad_before, pad_after), method, axis)
+
+    elif mode == 'symmetric':
+        for axis, (pad_before, pad_after) in enumerate(pad_width):
+            # Recursive padding along any axis where `pad_amt` is too large
+            # for indexing tricks. We can only safely pad the original axis
+            # length, to keep the period of the reflections consistent.
+            method = kwargs['reflect_type']
+            safe_pad = newmat.shape[axis]
+            repeat = safe_pad
+            while ((pad_before > safe_pad) or
+                   (pad_after > safe_pad)):
+                pad_iter_b = min(safe_pad,
+                                 safe_pad * (pad_before // safe_pad))
+                pad_iter_a = min(safe_pad, safe_pad * (pad_after // safe_pad))
+                newmat = _pad_sym(newmat, (pad_iter_b,
+                                           pad_iter_a), method, axis)
+                pad_before -= pad_iter_b
+                pad_after -= pad_iter_a
+                safe_pad += pad_iter_b + pad_iter_a
+            newmat = _pad_sym(newmat, (pad_before, pad_after), method, axis)
+
+    elif mode == 'wrap':
+        for axis, (pad_before, pad_after) in enumerate(pad_width):
+            # Recursive padding along any axis where `pad_amt` is too large
+            # for indexing tricks. We can only safely pad the original axis
+            # length, to keep the period of the reflections consistent.
+            safe_pad = newmat.shape[axis]
+            repeat = safe_pad
+            while ((pad_before > safe_pad) or
+                   (pad_after > safe_pad)):
+                pad_iter_b = min(safe_pad,
+                                 safe_pad * (pad_before // safe_pad))
+                pad_iter_a = min(safe_pad, safe_pad * (pad_after // safe_pad))
+                newmat = _pad_wrap(newmat, (pad_iter_b, pad_iter_a), axis)
+
+                pad_before -= pad_iter_b
+                pad_after -= pad_iter_a
+                safe_pad += pad_iter_b + pad_iter_a
+
+            newmat = _pad_wrap(newmat, (pad_before, pad_after), axis)
+
     return newmat