Maxlag2016b

doc/source/reference/c-api.array.rst

@@ -2161,6 +2161,32 @@ Array Functions
 
       z[k] = sum_n op1[n] * conj(op2[n+k])
 
+.. c:function:: PyObject* PyArray_CorrelateLags(PyObject* op1, PyObject* op2, npy_intp minlag, npy_intp maxlag, npy_intp lagstep)
+
+    Updated version of PyArray_Correlate2, which allows for custom lags in calculation.
+    The correlation is computed at each output point by multiplying *op1* by
+    a shifted version of *op2* and summing the result.
+    As a result of the shift, needed values outside of the defined range of
+    *op1* and *op2* are interpreted as zero.
+
+    This version differs from the other PyArray_Correlate* functions in that
+    PyArray_CorrelateLags calculates the correlation for lags starting at minlag
+    and ending at (maxlag - 1), with steps of size lagstep.
+    The lags for which the correlation will be calculated correspond with the values in
+    the vector formed by numpy.arange((minlag, maxlag, lagstep)).
+
+    .. rubric:: Notes
+
+    The following lag values should be used to reproduce the mode parameter used by
+    PyArray_Correlate2 (n1 and n2 are the vector lengths):
+       mode     |   minlag   |   maxlag   |   lagstep
+    0 ("valid") |     0      |  n1-n2+1   |      1
+     1 ("same") |   -n2/2    |  n1-n2/2   |      1
+     2 ("full") |   -n2+1    |     n1     |      1
+
+    PyArray_CorrelateLags uses the same definition of correlation as PyArray_Correlate2
+    in that the conjugate is taken and the argument order matters.
+
 .. c:function:: PyObject* PyArray_Where(PyObject* condition, PyObject* x, PyObject* y)
 
     If both ``x`` and ``y`` are ``NULL``, then return

numpy/core/code_generators/numpy_api.py

@@ -344,6 +344,8 @@
     # End 1.9 API
     'PyArray_CheckAnyScalarExact':          (300, NonNull(1)),
     # End 1.10 API
+    'PyArray_CorrelateLags':                (301,),
+    # End 1.11 API
 }
 
 ufunc_types_api = {

numpy/core/numeric.py

@@ -23,7 +23,7 @@
 from . import umath
 from .umath import (invert, sin, UFUNC_BUFSIZE_DEFAULT, ERR_IGNORE,
                     ERR_WARN, ERR_RAISE, ERR_CALL, ERR_PRINT, ERR_LOG,
-                    ERR_DEFAULT, PINF, NAN)
+                    ERR_DEFAULT, PINF, NAN, ceil)
 from . import numerictypes
 from .numerictypes import longlong, intc, int_, float_, complex_, bool_
 from ._internal import TooHardError
@@ -377,7 +377,8 @@ def extend_all(module):
 
 
 def asarray(a, dtype=None, order=None):
-    """Convert the input to an array.
+    """
+    Convert the input to an array.
 
     Parameters
     ----------
@@ -388,9 +389,8 @@ def asarray(a, dtype=None, order=None):
     dtype : data-type, optional
         By default, the data-type is inferred from the input data.
     order : {'C', 'F'}, optional
-        Whether to use row-major (C-style) or
-        column-major (Fortran-style) memory representation.
-        Defaults to 'C'.
+        Whether to use row-major ('C') or column-major ('F' for FORTRAN)
+        memory representation.  Defaults to 'C'.
 
     Returns
     -------
@@ -448,7 +448,8 @@ class ndarray is returned.
 
 
 def asanyarray(a, dtype=None, order=None):
-    """Convert the input to an ndarray, but pass ndarray subclasses through.
+    """
+    Convert the input to an ndarray, but pass ndarray subclasses through.
 
     Parameters
     ----------
@@ -459,8 +460,8 @@ def asanyarray(a, dtype=None, order=None):
     dtype : data-type, optional
         By default, the data-type is inferred from the input data.
     order : {'C', 'F'}, optional
-        Whether to use row-major (C-style) or column-major
-        (Fortran-style) memory representation.  Defaults to 'C'.
+        Whether to use row-major ('C') or column-major ('F') memory
+        representation.  Defaults to 'C'.
 
     Returns
     -------
@@ -814,16 +815,55 @@ def flatnonzero(a):
 
 _mode_from_name_dict = {'v': 0,
                         's': 1,
-                        'f': 2}
-
+                        'f': 2,
+                        'l': 3,
+                        }
 
 def _mode_from_name(mode):
-    if isinstance(mode, basestring):
+    if type(mode) is int:
+        assert(mode in _mode_from_name_dict.values())
+        return mode
+    elif isinstance(mode, basestring):
         return _mode_from_name_dict[mode.lower()[0]]
     return mode
 
-
-def correlate(a, v, mode='valid'):
+def _lags_from_lags(l):
+    if type(l) is int:          # maxlag
+        lags = (-l+1, l, 1)
+    elif type(l) is tuple:      # minlag and maxlag
+        if len(l) > 2:
+            lags = (int(l[0]), int(l[1]), int(l[2]))
+        else:
+            lags = (int(l[0]), int(l[1]), 1)
+    return lags
+
+def _lags_from_mode(alen, vlen, mode, lags=()):
+    inverted = 0
+    if alen < vlen:
+        alen, vlen = vlen, alen
+        inverted = 1
+
+    if mode is 0:
+        mode_lags = (0, alen-vlen+1, 1)
+    elif mode is 1:
+        mode_lags = (-int(vlen/2), alen - int(vlen/2), 1)
+    elif mode is 2:
+        mode_lags = (-vlen+1, alen, 1)
+    elif mode is 3:
+        mode_lags = _lags_from_lags(lags)
+
+    if inverted and mode != 3:
+        mode_lags = (-int(ceil((mode_lags[1]-mode_lags[0])/float(mode_lags[2])))
+                     *mode_lags[2]-mode_lags[0]+mode_lags[2],
+                     -mode_lags[0]+mode_lags[2], mode_lags[2])
+
+    if lags:
+        lags = _lags_from_lags(lags)
+        assert lags == mode_lags, \
+               'Custom lags should be passed to correlate() with mode="lags".'
+    return mode_lags
+
+def correlate(a, v, mode='default', lags=(), returns_lagvector=False):
     """
     Cross-correlation of two 1-dimensional sequences.
 
@@ -839,9 +879,18 @@ def correlate(a, v, mode='valid'):
     ----------
     a, v : array_like
         Input sequences.
-    mode : {'valid', 'same', 'full'}, optional
+    mode : {'valid', 'same', 'full', 'lags'}, optional
         Refer to the `convolve` docstring.  Note that the default
-        is 'valid', unlike `convolve`, which uses 'full'.
+        is `valid`, unlike `convolve`, which uses `full`.
+    lags : int or int tuple, optional
+        Refer to the `convolve` docstring.
+        mode should be unset or set to 'lags' to use the lags argument
+    returns_lagvector : bool, optional
+        If True, the function returns a lagvector array in addition to the
+        cross-correlation result.  The lagvector contains the indices of
+        the lags for which the cross-correlation was calculated.  It is
+        the same length as the return array, and corresponds one-to-one.
+        False is default.
     old_behavior : bool
         `old_behavior` was removed in NumPy 1.10. If you need the old
         behavior, use `multiarray.correlate`.
@@ -850,6 +899,9 @@ def correlate(a, v, mode='valid'):
     -------
     out : ndarray
         Discrete cross-correlation of `a` and `v`.
+    lagvector : ndarray, optional
+        The indices of the lags for which the cross-correlation was calculated.
+        It is the same length as out, and corresponds one-to-one.
 
     See Also
     --------
@@ -869,10 +921,14 @@ def correlate(a, v, mode='valid'):
     --------
     >>> np.correlate([1, 2, 3], [0, 1, 0.5])
     array([ 3.5])
-    >>> np.correlate([1, 2, 3], [0, 1, 0.5], "same")
+    >>> np.correlate([1, 2, 3], [0, 1, 0.5], mode="same")
     array([ 2. ,  3.5,  3. ])
-    >>> np.correlate([1, 2, 3], [0, 1, 0.5], "full")
-    array([ 0.5,  2. ,  3.5,  3. ,  0. ])
+    >>> np.correlate([1, 2, 3], [0, 1, 0.5], mode="full", returns_lagvector=True)
+    (array([ 0.5,  2. ,  3.5,  3. ,  0. ]), array([-2, -1,  0,  1,  2]))
+    >>> np.correlate([1, 2, 3], [0, 1, 0.5], mode="lags", lags=2)
+    array([ 2. ,  3.5,  3. ])
+    >>> np.correlate([1, 2, 3], [0, 1, 0.5], mode="lags", lags=(-1,2,2), returns_lagvector=True)
+    (array([ 2.,  3.]), array([-1,  1]))
 
     Using complex sequences:
 
@@ -887,11 +943,24 @@ def correlate(a, v, mode='valid'):
     array([ 0.0+0.j ,  3.0+1.j ,  1.5+1.5j,  1.0+0.j ,  0.5+0.5j])
 
     """
+    if mode == 'default':
+        if lags:
+            mode = 'lags'
+        else:
+            mode = 'valid'
     mode = _mode_from_name(mode)
-    return multiarray.correlate2(a, v, mode)
+    if mode in (0, 1, 2) and not returns_lagvector:
+        return multiarray.correlate2(a, v, mode)
 
+    lags = _lags_from_mode(len(a), len(v), mode, lags)
+    if returns_lagvector:
+        return multiarray.correlate2(a, v, 3, lags[0], lags[1], lags[2]), \
+                                     arange(lags[0], lags[1], lags[2])
+    else:
+        return multiarray.correlate2(a, v, 3, lags[0], lags[1], lags[2])
 
-def convolve(a,v,mode='full'):
+
+def convolve(a, v, mode='full', lags=(), returns_lagvector=False):
     """
     Returns the discrete, linear convolution of two one-dimensional sequences.
 
@@ -909,27 +978,59 @@ def convolve(a,v,mode='full'):
         First one-dimensional input array.
     v : (M,) array_like
         Second one-dimensional input array.
-    mode : {'full', 'valid', 'same'}, optional
+    mode : int, int tuple, or {'full', 'valid', 'same'}, optional
         'full':
           By default, mode is 'full'.  This returns the convolution
           at each point of overlap, with an output shape of (N+M-1,). At
           the end-points of the convolution, the signals do not overlap
-          completely, and boundary effects may be seen.
+          completely, and boundary effects may be seen.  This corresponds
+          with a lag tuple of (-M+1, N, 1) for N>M or (-N+1, M, 1)
+          for M>N.
 
         'same':
-          Mode 'same' returns output of length ``max(M, N)``.  Boundary
-          effects are still visible.
+          Mode `same` returns output of length ``max(M, N)``.  Boundary
+          effects are still visible. This corresponds with a lag tuple of
+          (-M/2, N-M/2, 1) for N>M or (-M+N/2+1, N/2+1, 1) for M>N.
 
         'valid':
           Mode 'valid' returns output of length
           ``max(M, N) - min(M, N) + 1``.  The convolution product is only given
           for points where the signals overlap completely.  Values outside
-          the signal boundary have no effect.
+          the signal boundary have no effect. This corresponds with a lag tuple
+          of (0, N-M+1, 1) for N>M or (-M+N, 1, 1) for M>N.
+
+        'lags':
+          Mode 'lags' uses the lags argument to define the lags for which
+          to perform the convolution.
+    lags : int or int tuple, optional
+        mode should be unset or set to 'lags' to use the lags argument
+
+        int (maxlag):
+          This calculates the convolution for all lags starting at
+          (-maxlag + 1) and ending at (maxlag - 1), with steps of size 1.
+          See the optional returns_lagvec argument to get an array containing
+          lags corresponding to the convolution values in the return array.
+
+        tuple (minlag, maxlag) or (minlag, maxlag, lagstep):
+          This calculates the convolution for all lags starting at
+          minlag and ending at (maxlag - 1), with steps of size lagstep.
+          The lags for which the convolution will be calculated correspond
+          with the values in the vector formed by numpy.arange() with the
+          same tuple argument.
+    returns_lagvector : bool, optional
+        If True, the function returns a lagvector array in addition to the
+        convolution result.  The lagvector contains the indices of
+        the lags for which the convolution was calculated.  It is
+        the same length as the return array, and corresponds one-to-one.
+        False is default.
 
     Returns
     -------
     out : ndarray
         Discrete, linear convolution of `a` and `v`.
+    lagvector : ndarray, optional
+        The indices of the lags for which the convolution was calculated.
+        It is the same length as out, and corresponds one-to-one.
 
     See Also
     --------
@@ -968,25 +1069,58 @@ def convolve(a,v,mode='full'):
     Contains boundary effects, where zeros are taken
     into account:
 
-    >>> np.convolve([1,2,3],[0,1,0.5], 'same')
+    >>> np.convolve([1,2,3],[0,1,0.5], mode='same')
     array([ 1. ,  2.5,  4. ])
 
     The two arrays are of the same length, so there
-    is only one position where they completely overlap:
+    is only one position where they completely overlap,
+    corresponding to a lag of 0.  lagvector=True causes
+    the function to return the lagvector corresponding
+    to the convolution in addition to the convolution
+    itself:
+
+    >>> np.convolve([1,2,3],[0,1,0.5], mode='valid', returns_lagvector=True)
+    (array([ 2.5]), array([0]))
 
-    >>> np.convolve([1,2,3],[0,1,0.5], 'valid')
-    array([ 2.5])
+    Find the convolution for lags ranging from -1 to 1
+    (0 is the lag for which the left sides of the arrays
+    are aligned, -1 has the second vector to the left of
+    the first, and +1 has the second vector to the right
+    of the first):
+
+    >>> np.convolve([1,2,3],[0,1,0.5], mode='lags', lags=2, returns_lagvector=True)
+    (array([ 1. ,  2.5,  4. ]), array([-1,  0,  1]))
+
+    Find the convolution for lags ranging from -2 to 4
+    with steps of length 2 (the maxlag member of the
+    lag range tuple is non-inclusive, similar to np.arange()):
+
+    >>> np.convolve([1,2,3,4,5],[0,1,0.5], mode='lags', lags=(-2,6,2), returns_lagvector=True)
+    (array([ 0. ,  2.5,  5.5,  2.5]), array([-2,  0,  2,  4]))
 
     """
     a, v = array(a, copy=False, ndmin=1), array(v, copy=False, ndmin=1)
-    if (len(v) > len(a)):
-        a, v = v, a
     if len(a) == 0:
         raise ValueError('a cannot be empty')
     if len(v) == 0:
         raise ValueError('v cannot be empty')
+    if (len(v) > len(a)):
+        a, v = v, a
+    if mode == 'default':
+        if lags:
+            mode = 'lags'
+        else:
+            mode = 'full'
     mode = _mode_from_name(mode)
-    return multiarray.correlate(a, v[::-1], mode)
+    if mode in (0, 1, 2) and not returns_lagvector:
+        return multiarray.correlate2(a, v[::-1], mode)
+
+    lags = _lags_from_mode(len(a), len(v), mode, lags)
+    if returns_lagvector:
+        return multiarray.correlate2(a, v[::-1], 3, lags[0], lags[1], lags[2]),\
+                                     arange(lags[0], lags[1], lags[2])
+    else:
+        return multiarray.correlate2(a, v[::-1], 3, lags[0], lags[1], lags[2])
 
 
 def outer(a, b, out=None):
@@ -1096,7 +1230,6 @@ def alterdot():
     restoredot : `restoredot` undoes the effects of `alterdot`.
 
     """
-    # 2014-08-13, 1.10
     warnings.warn("alterdot no longer does anything.", DeprecationWarning)
 
 
@@ -1120,7 +1253,6 @@ def restoredot():
     alterdot : `restoredot` undoes the effects of `alterdot`.
 
     """
-    # 2014-08-13, 1.10
     warnings.warn("restoredot no longer does anything.", DeprecationWarning)