Maxlag2023

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

@@ -2309,6 +2309,27 @@ 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)
 

numpy/core/code_generators/cversions.txt

@@ -71,3 +71,6 @@
 
 # Version 17 (NumPy 1.25) No actual change.
 0x00000011 = ca1aebdad799358149567d9d93cbca09
+
+# Version 18 (NumPy 1.24) Added PyArray_CorrelateLags.
+0x00000012 = 923ebe245d30b2d9c98027ddbb7789f4

numpy/core/code_generators/numpy_api.py

@@ -370,6 +370,8 @@ def get_annotations():
     'PyDataMem_SetHandler':                 (304, MinVersion("1.22")),
     'PyDataMem_GetHandler':                 (305, MinVersion("1.22")),
     # End 1.22 API
+    'PyArray_CorrelateLags':    (307, MinVersion("1.25")),
+    # End 1.26 API
 }
 
 ufunc_types_api = {

numpy/core/include/numpy/ndarraytypes.h

@@ -215,7 +215,8 @@ typedef enum {
 typedef enum {
         NPY_VALID=0,
         NPY_SAME=1,
-        NPY_FULL=2
+        NPY_FULL=2,
+        NPY_LAGS=3
 } NPY_CORRELATEMODE;
 
 /* The special not-a-time (NaT) value */

numpy/core/numeric.py

@@ -25,7 +25,7 @@
 from . import umath
 from . import shape_base
 from .overrides import set_array_function_like_doc, set_module
-from .umath import (multiply, invert, sin, PINF, NAN)
+from .umath import (multiply, invert, sin, PINF, NAN, ceil)
 from . import numerictypes
 from .numerictypes import longlong, intc, int_, float_, complex_, bool_
 from ..exceptions import ComplexWarning, TooHardError, AxisError
@@ -653,12 +653,73 @@ def flatnonzero(a):
     return np.nonzero(np.ravel(a))[0]
 
 
-def _correlate_dispatcher(a, v, mode=None):
+_mode_from_name_dict = {'v': 0,
+                        's': 1,
+                        'f': 2,
+                        'l': 3,
+                        }
+keys = [key for key in _mode_from_name_dict]
+for key in keys:
+    _mode_from_name_dict[key.upper()] = _mode_from_name_dict[key]
+_mode_from_name_dict_values = _mode_from_name_dict.values()
+
+
+def _mode_from_name(mode):
+    # guarantees that output is a value in _mode_from_name_dict
+    if mode in _mode_from_name_dict_values:
+        return mode
+    try:
+        mode = _mode_from_name_dict[mode[0]]
+    except KeyError:
+        raise ValueError("correlate/convolve mode argument must be unused or" +
+                         " one of {'valid', 'same', 'full', 'lags'}")
+    return mode
+
+
+def _lags_from_lags(lag):
+    if type(lag) is int:          # maxlag
+        lags = (-lag+1, lag, 1)
+    elif type(lag) is tuple:      # minlag and maxlag
+        if len(lag) > 2:
+            lags = (int(lag[0]), int(lag[1]), int(lag[2]))
+        else:
+            lags = (int(lag[0]), int(lag[1]), 1)
+    else:
+        raise ValueError("correlate/convolve lags argument must be " +
+                         "int or int tuple.")
+    return lags
+
+
+def _lags_from_mode(alen, vlen, mode):
+    inverted = 0
+    if alen < vlen:
+        alen, vlen = vlen, alen
+        inverted = 1
+
+    if mode == 0:
+        mode_lags = (0, alen-vlen+1, 1)
+    elif mode == 1:
+        mode_lags = (-int(vlen/2), alen - int(vlen/2), 1)
+    elif mode == 2:
+        mode_lags = (-vlen+1, alen, 1)
+    else:
+        raise ValueError("correlate/convolve mode argument must be unused or" +
+                         " one of {'valid', 'same', 'full', 'lags'}")
+
+    if inverted:
+        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])
+
+    return mode_lags
+
+
+def _correlate_dispatcher(a, v, mode=None, lags=None, returns_lagvector=None):
     return (a, v)
 
 
 @array_function_dispatch(_correlate_dispatcher)
-def correlate(a, v, mode='valid'):
+def correlate(a, v, mode='default', lags=None, returns_lagvector=False):
     r"""
     Cross-correlation of two 1-dimensional sequences.
 
@@ -674,9 +735,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`.
@@ -685,6 +755,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
     --------
@@ -710,10 +783,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")
-    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="same")
+    array([ 2. ,  3.5,  3. ])
+    >>> 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:
 
@@ -728,15 +805,36 @@ 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])
 
     """
-    return multiarray.correlate2(a, v, mode)
+    if mode == 'default':
+        if lags:
+            mode = 'lags'
+        else:
+            mode = 'valid'
+    mode = _mode_from_name(mode)  # guaranteed a value in _mode_from_name_dict
+    if mode in (0, 1, 2):
+        if lags:
+            raise ValueError("correlate mode keyword argument must be 'lags'" +
+                             " or unused if the lags keyword argument is used.")
+        result = multiarray.correlate2(a, v, mode)
+        if returns_lagvector:
+            alen, vlen = len(a), len(v)
+            lags = _lags_from_mode(alen, vlen, mode)
+    elif mode == 3:
+        lags = _lags_from_lags(lags)
+        result = multiarray.correlate2(a, v, 3, lags[0], lags[1], lags[2])
+
+    if returns_lagvector:
+        return result, arange(lags[0], lags[1], lags[2])
+    else:
+        return result
 
 
-def _convolve_dispatcher(a, v, mode=None):
+def _convolve_dispatcher(a, v, mode=None, lags=None, returns_lagvector=False):
     return (a, v)
 
 
 @array_function_dispatch(_convolve_dispatcher)
-def convolve(a, v, mode='full'):
+def convolve(a, v, mode='full', lags=None, returns_lagvector=False):
     """
     Returns the discrete, linear convolution of two one-dimensional sequences.
 
@@ -754,27 +852,60 @@ 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
     --------
@@ -814,24 +945,63 @@ 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:
-
-    >>> np.convolve([1,2,3],[0,1,0.5], 'valid')
-    array([2.5])
+    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]))
+
+    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:
+    alen, vlen = len(a), len(v)
+    if alen == 0:
         raise ValueError('a cannot be empty')
-    if len(v) == 0:
+    if vlen == 0:
         raise ValueError('v cannot be empty')
-    return multiarray.correlate(a, v[::-1], mode)
+
+    if mode == 'default':
+        if lags:
+            mode = 'lags'
+        else:
+            mode = 'full'
+    mode = _mode_from_name(mode)  # guaranteed a value in _mode_from_name_dict
+    if mode in (0, 1, 2):
+        if lags:
+            raise ValueError("convolve mode keyword argument must be 'lags' " +
+                             "or unused if the lags keyword argument is used.")
+        result = multiarray.correlate2(a, v[::-1], mode)
+        if returns_lagvector:
+            lags = _lags_from_mode(alen, vlen, mode)
+    elif mode == 3:
+        lags = _lags_from_lags(lags)
+        result = multiarray.correlate2(a, v[::-1], 3, lags[0], lags[1], lags[2])
+
+    if returns_lagvector:
+        return result, arange(lags[0], lags[1], lags[2])
+    else:
+        return result
 
 
 def _outer_dispatcher(a, b, out=None):

numpy/core/setup_common.py

@@ -49,7 +49,8 @@
 # 0x00000010 - 1.23.x
 # 0x00000010 - 1.24.x
 # 0x00000011 - 1.25.x
-C_API_VERSION = 0x00000011
+# 0x00000012 - 1.26.x
+C_API_VERSION = 0x00000012
 
 # By default, when compiling downstream libraries against NumPy,```
 # pick an older feature version.  For example, for 1.25.x we default to the