From c5f43db25e37c6b48bde5bb7b9a1b6989597bbc7 Mon Sep 17 00:00:00 2001
From: Honi Sanders <sandersh6000@gmail.com>
Date: Mon, 22 May 2023 09:38:02 -0400
Subject: [PATCH 1/9] starting off work with files from maxlag. still need to:
 - fix correlatemode_parser and PyArray_CorrelatemodeConverter in
 core/src/multiarray/conversion_utils.c - fix array_correlate and
 array_correlate2 in core/src/multiarray/multiarraymodule.c based on those
 fixes and based on maxlag PR

---
 doc/source/reference/c-api/array.rst         |  21 ++
 numpy/core/code_generators/numpy_api.py      |   2 +
 numpy/core/src/multiarray/multiarraymodule.c | 307 +++++++++++++++----
 numpy/core/tests/test_numeric.py             | 122 +++++++-
 4 files changed, 382 insertions(+), 70 deletions(-)

diff --git a/doc/source/reference/c-api/array.rst b/doc/source/reference/c-api/array.rst
index baebfdb0217d..1cb97f94afde 100644
--- a/doc/source/reference/c-api/array.rst
+++ b/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)
 
diff --git a/numpy/core/code_generators/numpy_api.py b/numpy/core/code_generators/numpy_api.py
index bb3b15806a55..a9ee9f23cd35 100644
--- a/numpy/core/code_generators/numpy_api.py
+++ b/numpy/core/code_generators/numpy_api.py
@@ -424,6 +424,8 @@ def get_annotations():
     # End 1.8 API
     'PyUFunc_FromFuncAndDataAndSignatureAndIdentity': (42, MinVersion("1.16")),
     # End 1.16 API
+    'PyArray_CorrelateLags':    (43,),
+    # End 1.17 API
 }
 
 # List of all the dicts which define the C API
diff --git a/numpy/core/src/multiarray/multiarraymodule.c b/numpy/core/src/multiarray/multiarraymodule.c
index d7d19493b754..3a37a1db04b2 100644
--- a/numpy/core/src/multiarray/multiarraymodule.c
+++ b/numpy/core/src/multiarray/multiarraymodule.c
@@ -1210,18 +1210,25 @@ PyArray_CopyAndTranspose(PyObject *op)
  * inverted is set to 1 if computed correlate(ap2, ap1), 0 otherwise
  */
 static PyArrayObject*
-_pyarray_correlate(PyArrayObject *ap1, PyArrayObject *ap2, int typenum,
-                   int mode, int *inverted)
+_pyarray_correlate(PyArrayObject *ap1, PyArrayObject *ap2,
+                   int typenum, int *inverted,
+                   npy_intp minlag, npy_intp maxlag, npy_intp lagstep)
 {
     PyArrayObject *ret;
     npy_intp length;
-    npy_intp i, n1, n2, n, n_left, n_right;
+    npy_intp i, i1, n1, n2, n, n11;
+    npy_intp lag, tmplag, maxleft, maxright;
     npy_intp is1, is2, os;
     char *ip1, *ip2, *op;
     PyArray_DotFunc *dot;
 
     NPY_BEGIN_THREADS_DEF;
 
+    if (lagstep == 0) {
+        lagstep = 1;
+    }
+
+    /* size of x (n1) and y (n2) */
     n1 = PyArray_DIMS(ap1)[0];
     n2 = PyArray_DIMS(ap2)[0];
     if (n1 == 0) {
@@ -1239,41 +1246,42 @@ _pyarray_correlate(PyArrayObject *ap1, PyArrayObject *ap2, int typenum,
         ret = NULL;
         i = n1;
         n1 = n2;
-        n2 = i;
+        n2 = i;        
+        minlag = -minlag;
+        maxlag = -maxlag;
+        lagstep = -lagstep;
+    }
+
+    if (lagstep < 0) {
         *inverted = 1;
-    } else {
+        i = minlag;
+        i1 = (npy_intp)(npy_ceil((maxlag - minlag)/(float)lagstep))*lagstep;
+        minlag =  i1 + minlag - lagstep;
+        maxlag = i - lagstep;
+        lagstep = -lagstep;
+    }
+    else {
         *inverted = 0;
     }
-
-    length = n1;
-    n = n2;
-    switch(mode) {
-    case 0:
-        length = length - n + 1;
-        n_left = n_right = 0;
-        break;
-    case 1:
-        n_left = (npy_intp)(n/2);
-        n_right = n - n_left - 1;
-        break;
-    case 2:
-        n_right = n - 1;
-        n_left = n - 1;
-        length = length + n - 1;
-        break;
-    default:
-        PyErr_SetString(PyExc_ValueError, "mode must be 0, 1, or 2");
-        return NULL;
+    if (maxlag <= minlag) {
+        length = 0;
+    }
+    else {
+        length = (maxlag - minlag + lagstep - 1)/lagstep;
     }
 
     /*
      * Need to choose an output array that can hold a sum
      * -- use priority to determine which subtype.
+     * ret is the array that will be returned as the answer
      */
     ret = new_array_for_sum(ap1, ap2, NULL, 1, &length, typenum, NULL);
     if (ret == NULL) {
         return NULL;
     }
+    if (length > 0) {
+        PyArray_FILLWBYTE(ret, 0);
+    }
     dot = PyArray_DESCR(ret)->f->dotfunc;
     if (dot == NULL) {
         PyErr_SetString(PyExc_ValueError,
@@ -1282,36 +1290,67 @@ _pyarray_correlate(PyArrayObject *ap1, PyArrayObject *ap2, int typenum,
     }
 
     NPY_BEGIN_THREADS_DESCR(PyArray_DESCR(ret));
+    /* p is pointer to array start and s is stride */
+    /* i1 is x, i2, is y, o is answer */
+    ip1 = PyArray_DATA(ap1);
     is1 = PyArray_STRIDES(ap1)[0];
+    ip2 = PyArray_DATA(ap2);
     is2 = PyArray_STRIDES(ap2)[0];
     op = PyArray_DATA(ret);
     os = PyArray_DESCR(ret)->elsize;
-    ip1 = PyArray_DATA(ap1);
-    ip2 = PyArray_BYTES(ap2) + n_left*is2;
-    n = n - n_left;
-    for (i = 0; i < n_left; i++) {
-        dot(ip1, is1, ip2, is2, op, n, ret);
-        n++;
-        ip2 -= is2;
+
+    lag = minlag;
+    if (lag < -n2+1) {
+        /* if minlag is before any overlap between the vectors,
+         * then skip to first relevant lag
+         */
+        op += os*((-n2+1) - lag);
+        lag = -n2+1;
+    }
+    maxleft = (0 < maxlag ? 0 : maxlag);
+    tmplag = lag;
+    for (lag = tmplag; lag < maxleft; lag+=lagstep) {
+        /* for lags where y is left of x,
+         * overlap is length of y - lag
+         */
+        n = n2 + lag;
+        dot(ip1, is1, ip2 - lag*is2, is2, op, n, ret);
+        /* iterate over answer vector */
         op += os;
     }
-    if (small_correlate(ip1, is1, n1 - n2 + 1, PyArray_TYPE(ap1),
-                        ip2, is2, n, PyArray_TYPE(ap2),
-                        op, os)) {
-        ip1 += is1 * (n1 - n2 + 1);
-        op += os * (n1 - n2 + 1);
+    if (maxlag < n1 - n2) {
+        /* if maxlag doesn't take y all the way to the end of x,
+         * relevant length of x is smaller
+         */
+        n11 = maxlag + n2;
     }
     else {
-        for (i = 0; i < (n1 - n2 + 1); i++) {
-            dot(ip1, is1, ip2, is2, op, n, ret);
-            ip1 += is1;
+        n11 = n1;
+    }
+    /* starts at lag=minlag if minlag>0.
+     * Does lags where y entirely overlaps with x.
+     */
+    if (lagstep == 1 && lag < maxlag &&
+            small_correlate(ip1 + lag*is1, is1,
+                            n11 - n2 + 1 - lag, PyArray_TYPE(ap1),
+                            ip2, is2, n2, PyArray_TYPE(ap2),
+                            op, os)) {
+        lag = n11 - n2 + 1;
+        op += os * (n11 - n2 + 1);
+    }
+    else if (lag < maxlag) {
+        tmplag = lag;
+        for (lag = tmplag; lag < (n11 - n2 + 1); lag+=lagstep) {
+            dot(ip1 + lag*is1, is1, ip2, is2, op, n2, ret);
             op += os;
         }
     }
-    for (i = 0; i < n_right; i++) {
-        n--;
-        dot(ip1, is1, ip2, is2, op, n, ret);
-        ip1 += is1;
+    maxright = (maxlag < n1  ? maxlag : n1 );
+    tmplag = lag;
+    for (lag = tmplag; lag < maxright; lag+=lagstep) {
+        /* for lags where y is right of x */
+        n = n1 - lag;
+        dot(ip1 + lag*is1, is1, ip2, is2, op, n, ret);
         op += os;
     }
 
@@ -1374,14 +1413,70 @@ _pyarray_revert(PyArrayObject *ret)
     return 0;
 }
 
+/*
+ * Generate the lags corresponding to each mode
+ * n1 and n2 are sizes of two vectors
+ * minlag, maxlag, and lagstep are edited in-place
+ */
+void
+_lags_from_mode(int mode, npy_intp n1, npy_intp n2,
+                npy_intp *minlag, npy_intp *maxlag, npy_intp *lagstep) {
+
+    int i, i1, inverted = 0;
+    /* place holders for minlag, maxlag, and lagstep respectively */
+    npy_intp m0, m1, s;
+
+    if (n1 < n2) {
+        i = n1;
+        n1 = n2;
+        n2 = i;
+        inverted = 1;
+    }
+
+    s = 1;
+    switch(mode) {
+    case 0:
+        /* mode = 'valid' */
+        m0 = 0;
+        m1 = n1 - n2 + 1;
+        break;
+    case 1:
+        /* mode = 'same' */
+        m0 = -(npy_intp)(n2/2);
+        m1 = n1 + m0;
+        break;
+    case 2:
+        /* mode = 'full' */
+        m0 = -n2 + 1;
+        m1 = n1;
+        break;
+    default:
+        PyErr_SetString(PyExc_ValueError, "mode must be 0, 1, or 2");
+    }
+
+    if (inverted) {
+        i = m0;
+        i1 = (npy_intp)(npy_ceil((m1 - m0)/(float)s))*s;
+        m0 =  -(i1 + m0 - s);
+        m1 = -(i - s);
+    }
+
+    /* set minlag, maxlag, and lagstep values */
+    *minlag = m0;
+    *maxlag = m1;
+    *lagstep = s;
+}
+
 /*NUMPY_API
- * correlate(a1,a2,mode)
+ * correlate(a1,a2,minlag,maxlag,lagstep)
  *
- * This function computes the usual correlation (correlate(a1, a2) !=
- * correlate(a2, a1), and conjugate the second argument for complex inputs
+ * This function computes the cross-correlation of a1 with a2.
+ * This function does not accept modes.
+ * See _lags_from_mode() to convert modes to lags.
  */
 NPY_NO_EXPORT PyObject *
-PyArray_Correlate2(PyObject *op1, PyObject *op2, int mode)
+PyArray_CorrelateLags(PyObject *op1, PyObject *op2,
+                      npy_intp minlag, npy_intp maxlag, npy_intp lagstep)
 {
     PyArrayObject *ap1, *ap2, *ret = NULL;
     int typenum;
@@ -1389,14 +1484,8 @@ PyArray_Correlate2(PyObject *op1, PyObject *op2, int mode)
     int inverted;
     int st;
 
-    typenum = PyArray_ObjectType(op1, NPY_NOTYPE);
-    if (typenum == NPY_NOTYPE) {
-        return NULL;
-    }
+    typenum = PyArray_ObjectType(op1, 0);
     typenum = PyArray_ObjectType(op2, typenum);
-    if (typenum == NPY_NOTYPE) {
-        return NULL;
-    }
 
     typec = PyArray_DescrFromType(typenum);
     Py_INCREF(typec);
@@ -1422,18 +1511,19 @@ PyArray_Correlate2(PyObject *op1, PyObject *op2, int mode)
         ap2 = cap2;
     }
 
-    ret = _pyarray_correlate(ap1, ap2, typenum, mode, &inverted);
+    ret = _pyarray_correlate(ap1, ap2, typenum, &inverted,
+                             minlag, maxlag, lagstep);
     if (ret == NULL) {
         goto clean_ap2;
     }
 
     /*
-     * If we inverted input orders, we need to reverse the output array (i.e.
-     * ret = ret[::-1])
+     * If we inverted input orders, we need to reverse the output array
+     * (i.e. ret = ret[::-1])
      */
     if (inverted) {
         st = _pyarray_revert(ret);
-        if (st) {
+        if(st) {
             goto clean_ret;
         }
     }
@@ -1451,6 +1541,86 @@ PyArray_Correlate2(PyObject *op1, PyObject *op2, int mode)
     return NULL;
 }
 
+/*NUMPY_API
+ * correlate(a1,a2,mode)
+ *
+ * This function computes the usual correlation (correlate(a1, a2) !=
+ * correlate(a2, a1), and conjugate the second argument for complex inputs
+ *
+ * This function requires a mode.
+ */
+NPY_NO_EXPORT PyObject *
+PyArray_Correlate2(PyObject *op1, PyObject *op2, int mode)
+{
+    PyArrayObject *ap1, *ap2, *ret = NULL;
+    int typenum;
+    PyArray_Descr *typec;
+    npy_intp minlag, maxlag, lagstep;
+    npy_intp n1, n2;
+    int inverted, st;
+
+    typenum = PyArray_ObjectType(op1, 0);
+    typenum = PyArray_ObjectType(op2, typenum);
+
+    typec = PyArray_DescrFromType(typenum);
+    Py_INCREF(typec);
+    ap1 = (PyArrayObject *)PyArray_FromAny(op1, typec, 1, 1,
+                                           NPY_ARRAY_DEFAULT, NULL);
+    if (ap1 == NULL) {
+        Py_DECREF(typec);
+        return NULL;
+    }
+    ap2 = (PyArrayObject *)PyArray_FromAny(op2, typec, 1, 1,
+                                           NPY_ARRAY_DEFAULT, NULL);
+    if (ap2 == NULL) {
+        goto clean_ap1;
+    }
+
+    if (PyArray_ISCOMPLEX(ap2)) {
+        PyArrayObject *cap2;
+        cap2 = (PyArrayObject *)PyArray_Conjugate(ap2, NULL);
+        if (cap2 == NULL) {
+            goto clean_ap2;
+        }
+        Py_DECREF(ap2);
+        ap2 = cap2;
+    }
+
+    /* size of x (n1) and y (n2) */
+    n1 = PyArray_DIMS(ap1)[0];
+    n2 = PyArray_DIMS(ap2)[0];
+    _lags_from_mode(mode, n1, n2, &minlag, &maxlag, &lagstep);
+    ret = _pyarray_correlate(ap1, ap2, typenum, &inverted,
+                             minlag, maxlag, lagstep);
+    if (ret == NULL) {
+        goto clean_ap2;
+    }
+
+    /*
+     * If we inverted input orders, we need to reverse the output array
+     * (i.e. ret = ret[::-1])
+     */
+    if (inverted) {
+        st = _pyarray_revert(ret);
+        if(st) {
+            goto clean_ret;
+        }
+    }
+
+    Py_DECREF(ap1);
+    Py_DECREF(ap2);
+    return (PyObject *)ret;
+
+clean_ret:
+    Py_DECREF(ret);
+clean_ap2:
+    Py_DECREF(ap2);
+clean_ap1:
+    Py_DECREF(ap1);
+    return NULL;
+}
+
+
 /*NUMPY_API
  * Numeric.correlate(a1,a2,mode)
  */
@@ -1461,15 +1631,11 @@ PyArray_Correlate(PyObject *op1, PyObject *op2, int mode)
     int typenum;
     int unused;
     PyArray_Descr *typec;
+    npy_intp minlag, maxlag, lagstep;
+    npy_intp n1, n2;
 
-    typenum = PyArray_ObjectType(op1, NPY_NOTYPE);
-    if (typenum == NPY_NOTYPE) {
-        return NULL;
-    }
+    typenum = PyArray_ObjectType(op1, 0);
     typenum = PyArray_ObjectType(op2, typenum);
-    if (typenum == NPY_NOTYPE) {
-        return NULL;
-    }
 
     typec = PyArray_DescrFromType(typenum);
     Py_INCREF(typec);
@@ -1485,8 +1651,13 @@ PyArray_Correlate(PyObject *op1, PyObject *op2, int mode)
         goto fail;
     }
 
-    ret = _pyarray_correlate(ap1, ap2, typenum, mode, &unused);
-    if (ret == NULL) {
+    /* size of x (n1) and y (n2) */
+    n1 = PyArray_DIMS(ap1)[0];
+    n2 = PyArray_DIMS(ap2)[0];
+    _lags_from_mode(mode, n1, n2, &minlag, &maxlag, &lagstep);
+    ret = _pyarray_correlate(ap1, ap2, typenum, &unused,
+                             minlag, maxlag, lagstep);
+    if(ret == NULL) {
         goto fail;
     }
     Py_DECREF(ap1);
diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py
index 832a47c926cf..3309237a2696 100644
--- a/numpy/core/tests/test_numeric.py
+++ b/numpy/core/tests/test_numeric.py
@@ -2942,21 +2942,66 @@ def _setup(self, dt):
         self.xs = np.arange(1, 20)[::3]
         self.y = np.array([-1, -2, -3], dtype=dt)
         self.z1 = np.array([-3., -8., -14., -20., -26., -14., -5.], dtype=dt)
+        self.z1s = np.array([-8., -14., -20., -26., -14.], dtype=dt)
+        self.z1v = np.array([-14., -20., -26.], dtype=dt)
         self.z1_4 = np.array([-2., -5., -8., -11., -14., -5.], dtype=dt)
         self.z1r = np.array([-15., -22., -22., -16., -10., -4., -1.], dtype=dt)
         self.z2 = np.array([-5., -14., -26., -20., -14., -8., -3.], dtype=dt)
+        self.z2s = np.array([-14, -26., -20., -14., -8.], dtype=dt)
+        self.z2v = np.array([-26., -20., -14.], dtype=dt)
         self.z2r = np.array([-1., -4., -10., -16., -22., -22., -15.], dtype=dt)
         self.zs = np.array([-3., -14., -30., -48., -66., -84.,
                            -102., -54., -19.], dtype=dt)
 
     def test_float(self):
         self._setup(float)
-        z = np.correlate(self.x, self.y, 'full')
+        z, lagvec = np.correlate(self.x, self.y, 'full', returns_lagvector=True)
         assert_array_almost_equal(z, self.z1)
+        assert_equal(lagvec.size, z.size)
+        assert_array_equal(lagvec, np.arange(-self.y.size + 1, self.x.size))
+        z, lagvec = np.correlate(self.x, self.y, 'same', returns_lagvector=True)
+        assert_array_almost_equal(z, self.z1s)
+        assert_array_almost_equal(z, self.z1[lagvec+self.y.size-1])
+        assert_equal(lagvec.size, z.size)
+        if self.x.size > self.y.size:
+            assert_array_equal(lagvec,
+                np.arange(-int(self.y.size/2), self.x.size - int(self.y.size/2)))
+        else:
+            assert_array_equal(lagvec,
+                np.arange(-self.y.size + int(self.x.size/2) + 1, int(self.x.size/2) + 1))
+        z, lagvec = np.correlate(self.x, self.y, 'valid', returns_lagvector=True)
+        assert_array_almost_equal(z, self.z1v)
+        assert_array_almost_equal(z, self.z1[lagvec+self.y.size-1])
+        assert_equal(lagvec.size, z.size)
+        if self.x.size > self.y.size:
+            assert_array_equal(lagvec, np.arange(0, self.x.size - self.y.size + 1))
+        else:
+            assert_array_equal(lagvec, np.arange(self.x.size - self.y.size, 1))
+        
         z = np.correlate(self.x, self.y[:-1], 'full')
         assert_array_almost_equal(z, self.z1_4)
-        z = np.correlate(self.y, self.x, 'full')
+        z, lagvec = np.correlate(self.y, self.x, 'full', returns_lagvector=True)
         assert_array_almost_equal(z, self.z2)
+        assert_equal(lagvec.size, z.size)
+        assert_array_equal(lagvec, np.arange(-self.x.size + 1, self.y.size))
+        z, lagvec = np.correlate(self.y, self.x, 'same', returns_lagvector=True)
+        assert_array_almost_equal(z, self.z2s)
+        assert_array_almost_equal(z, self.z2[lagvec+self.x.size-1])
+        assert_equal(lagvec.size, z.size)
+        if self.y.size > self.x.size:
+            assert_array_equal(lagvec,
+                np.arange(-int(self.x.size/2), self.y.size - int(self.x.size/2)))
+        else:
+            assert_array_equal(lagvec,
+                np.arange(-self.x.size + int(self.y.size/2) + 1, int(self.y.size/2) + 1))
+        z, lagvec = np.correlate(self.y, self.x, 'valid', returns_lagvector=True)
+        assert_array_almost_equal(z, self.z2v)
+        assert_array_almost_equal(z, self.z2[lagvec+self.x.size-1])
+        assert_equal(lagvec.size, z.size)
+        if self.y.size > self.x.size:
+            assert_array_equal(lagvec, np.arange(0, self.y.size - self.x.size + 1))
+        else:
+            assert_array_equal(lagvec, np.arange(self.y.size - self.x.size, 1))        z = np.correlate(self.x[::-1], self.y, 'full')
         z = np.correlate(self.x[::-1], self.y, 'full')
         assert_array_almost_equal(z, self.z1r)
         z = np.correlate(self.y, self.x[::-1], 'full')
@@ -2964,6 +3009,79 @@ def test_float(self):
         z = np.correlate(self.xs, self.y, 'full')
         assert_array_almost_equal(z, self.zs)
 
+     def test_lags(self):
+         self._setup(np.float)
+         longlen = max(self.x.size, self.y.size)
+         maxlag = 3 + longlen
+         minlag = -2
+         lagstep = 2
+         tmp = np.correlate(self.x, self.y, 'full')
+         z_full = np.zeros(tmp.size + 100)
+         z_full[:tmp.size] = tmp
+         z, lagvec = np.correlate(self.x, self.y, lags=maxlag,
+                                  returns_lagvector=True)
+         assert_array_equal(z[0:3], np.zeros(3))
+         assert_array_equal(z[-3:], np.zeros(3))
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(-maxlag+1, maxlag))
+         z, lagvec = np.correlate(self.x, self.y, mode='lags',
+                                  lags=(minlag, maxlag), returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(minlag, maxlag))
+         z, lagvec = np.correlate(self.x, self.y, lags=(minlag, maxlag, lagstep),
+                                  returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(minlag, maxlag, lagstep))
+         lagstep = 3
+         z, lagvec = np.correlate(self.x, self.y, mode='lags',
+                                  lags=(minlag, maxlag, lagstep),
+                                  returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         minlag = 10
+         maxlag = -2
+         lagstep = -2
+         z, lagvec = np.correlate(self.x, self.y, mode='lags',
+                                  lags=(minlag, maxlag, lagstep),
+                                  returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(int(minlag), int(maxlag), lagstep))
+         maxlag = 10.2
+         minlag = -2.2
+         z, lagvec = np.correlate(self.x, self.y, mode='lags',
+                                  lags=(minlag, maxlag),
+                                  returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(int(minlag), int(maxlag)))
+         z, lagvec = np.correlate(self.x, self.y, mode='lags',
+                                  lags=(maxlag, minlag),
+                                  returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(int(maxlag), int(minlag)))
+         maxlag = -1
+         minlag = -2
+         z, lagvec = np.correlate(self.x, self.y, mode='lags',
+                                  lags=(minlag, maxlag),
+                                  returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(int(minlag), int(maxlag)))
+         maxlag = 2
+         minlag = 4
+         z, lagvec = np.correlate(self.x, self.y, mode='lags',
+                                  lags=(minlag, maxlag),
+                                  returns_lagvector=True)
+         assert_array_almost_equal(z, z_full[lagvec+self.y.size-1])
+         assert_equal(lagvec.size, z.size)
+         assert_array_equal(lagvec, np.arange(int(minlag), int(maxlag)))
+         # check non-integer lagsteps?
+         # i want to be able to switch x and y and test
+
     def test_object(self):
         self._setup(Decimal)
         z = np.correlate(self.x, self.y, 'full')

From 3ec6e8c86bd7e14cef6e16aad5d4679932117cd4 Mon Sep 17 00:00:00 2001
From: Honi Sanders <sandersh6000@gmail.com>
Date: Mon, 22 May 2023 11:13:43 -0400
Subject: [PATCH 2/9] some changes to numeric.py from maxlag branch

---
 numpy/core/numeric.py | 217 +++++++++++++++++++++++++++++++++++++-----
 1 file changed, 192 insertions(+), 25 deletions(-)

diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index 91ac3f8606fe..85ffe8b21ff1 100644
--- a/numpy/core/numeric.py
+++ b/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,70 @@ def flatnonzero(a):
     return np.nonzero(np.ravel(a))[0]
 
 
+_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(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)
+    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 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)
+    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):
     return (a, v)
 
 
 @array_function_dispatch(_correlate_dispatcher)
-def correlate(a, v, mode='valid'):
+def correlate(a, v, mode='default', lags=(), returns_lagvector=False):
     r"""
     Cross-correlation of two 1-dimensional sequences.
 
@@ -674,9 +732,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 +752,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 +780,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,7 +802,28 @@ 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):
@@ -736,7 +831,7 @@ def _convolve_dispatcher(a, v, mode=None):
 
 
 @array_function_dispatch(_convolve_dispatcher)
-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.
 
@@ -754,27 +849,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 +942,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):

From bdfadbba5390067aafa4e223696ea697317e9727 Mon Sep 17 00:00:00 2001
From: Honi Sanders <sandersh6000@gmail.com>
Date: Mon, 22 May 2023 11:17:45 -0400
Subject: [PATCH 3/9] formatting  in numeric.py

---
 numpy/core/numeric.py | 37 ++++++++++++++++++++-----------------
 1 file changed, 20 insertions(+), 17 deletions(-)

diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index 85ffe8b21ff1..aa1f3cd93aad 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -663,30 +663,33 @@ def flatnonzero(a):
     _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(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]))
+    # 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(l[0]), int(l[1]), 1)
+            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:

From 9e52bd9d03b857303c76fe7eb750c8f0980ebba9 Mon Sep 17 00:00:00 2001
From: Honi Sanders <sandersh6000@gmail.com>
Date: Mon, 22 May 2023 11:53:16 -0400
Subject: [PATCH 4/9] more formatting in numeric.py

---
 numpy/core/numeric.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index aa1f3cd93aad..9a817a1640f8 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -696,16 +696,16 @@ def _lags_from_mode(alen, vlen, mode):
         alen, vlen = vlen, alen
         inverted = 1
 
-    if mode is 0:
+    if mode == 0:
         mode_lags = (0, alen-vlen+1, 1)
-    elif mode is 1:
+    elif mode == 1:
         mode_lags = (-int(vlen/2), alen - int(vlen/2), 1)
-    elif mode is 2:
+    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],
@@ -814,7 +814,7 @@ def correlate(a, v, mode='default', lags=(), returns_lagvector=False):
     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.")
+                             " or unused if the lags keyword argument is used.")
         result = multiarray.correlate2(a, v, mode)
         if returns_lagvector:
             alen, vlen = len(a), len(v)
@@ -876,10 +876,10 @@ def convolve(a, v, mode='full', lags=(), returns_lagvector=False):
         '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.

From 955f18c18b2659017548b70d260795c61633aa48 Mon Sep 17 00:00:00 2001
From: Honi Sanders <sandersh6000@gmail.com>
Date: Mon, 22 May 2023 12:29:47 -0400
Subject: [PATCH 5/9] formatting in numeric.py

---
 numpy/core/numeric.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index 9a817a1640f8..c1fbb2dcb0e7 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -989,8 +989,8 @@ def convolve(a, v, mode='full', lags=(), returns_lagvector=False):
     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.")
+            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)

From 42e5549c4abedd9f398fc407ad11465fd8b368a9 Mon Sep 17 00:00:00 2001
From: Honi Sanders <sandersh6000@gmail.com>
Date: Mon, 22 May 2023 13:30:27 -0400
Subject: [PATCH 6/9] getting numpy api to work, finishing incorporation of
 multiarraymodule.c changes from maxlag branch, requiring changes to
 conversion_utils.c and ndarratypes.h

---
 numpy/core/code_generators/cversions.txt     |  3 ++
 numpy/core/code_generators/numpy_api.py      |  4 +--
 numpy/core/include/numpy/ndarraytypes.h      |  3 +-
 numpy/core/setup_common.py                   |  3 +-
 numpy/core/src/multiarray/conversion_utils.c | 13 +++++---
 numpy/core/src/multiarray/multiarraymodule.c | 32 ++++++++++++++++++--
 6 files changed, 47 insertions(+), 11 deletions(-)

diff --git a/numpy/core/code_generators/cversions.txt b/numpy/core/code_generators/cversions.txt
index e52193d7a462..1eea0dfe198b 100644
--- a/numpy/core/code_generators/cversions.txt
+++ b/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
diff --git a/numpy/core/code_generators/numpy_api.py b/numpy/core/code_generators/numpy_api.py
index a9ee9f23cd35..3adf26e3b8a2 100644
--- a/numpy/core/code_generators/numpy_api.py
+++ b/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 = {
@@ -424,8 +426,6 @@ def get_annotations():
     # End 1.8 API
     'PyUFunc_FromFuncAndDataAndSignatureAndIdentity': (42, MinVersion("1.16")),
     # End 1.16 API
-    'PyArray_CorrelateLags':    (43,),
-    # End 1.17 API
 }
 
 # List of all the dicts which define the C API
diff --git a/numpy/core/include/numpy/ndarraytypes.h b/numpy/core/include/numpy/ndarraytypes.h
index 742ba5261225..98d5fc16d5d3 100644
--- a/numpy/core/include/numpy/ndarraytypes.h
+++ b/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 */
diff --git a/numpy/core/setup_common.py b/numpy/core/setup_common.py
index b5bc0dec1f89..d66aeb47fba5 100644
--- a/numpy/core/setup_common.py
+++ b/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
diff --git a/numpy/core/src/multiarray/conversion_utils.c b/numpy/core/src/multiarray/conversion_utils.c
index 84bf08206f51..ac214ce107d7 100644
--- a/numpy/core/src/multiarray/conversion_utils.c
+++ b/numpy/core/src/multiarray/conversion_utils.c
@@ -842,6 +842,10 @@ static int correlatemode_parser(char const *str, Py_ssize_t length, void *data)
         *val = NPY_FULL;
         is_exact = (length == 4 && strcmp(str, "full") == 0);
     }
+    else if (str[0] == 'L' || str[0] == 'l') {
+        *val = NPY_LAGS;
+        is_exact = (length == 4 && strcmp(str, "lags") == 0);
+    }
     else {
         return -1;
     }
@@ -853,7 +857,8 @@ static int correlatemode_parser(char const *str, Py_ssize_t length, void *data)
         /* Numpy 1.21, 2021-01-19 */
         if (DEPRECATE("inexact matches and case insensitive matches for "
                       "convolve/correlate mode are deprecated, please "
-                      "use one of 'valid', 'same', or 'full' instead.") < 0) {
+                      "use one of 'valid', 'same', 'full', or 'lags' "
+                      "instead.") < 0) {
             return -1;
         }
     }
@@ -870,7 +875,7 @@ PyArray_CorrelatemodeConverter(PyObject *object, NPY_CORRELATEMODE *val)
     if (PyUnicode_Check(object)) {
         return string_converter_helper(
             object, (void *)val, correlatemode_parser, "mode",
-            "must be one of 'valid', 'same', or 'full'");
+            "must be one of 'valid', 'same', 'full', or 'lags'");
     }
 
     else {
@@ -881,14 +886,14 @@ PyArray_CorrelatemodeConverter(PyObject *object, NPY_CORRELATEMODE *val)
                         "convolve/correlate mode not understood");
             return NPY_FAIL;
         }
-        if (number <= (int) NPY_FULL
+        if (number <= (int) NPY_LAGS
                 && number >= (int) NPY_VALID) {
             *val = (NPY_CORRELATEMODE) number;
             return NPY_SUCCEED;
         }
         else {
             PyErr_Format(PyExc_ValueError,
-                    "integer convolve/correlate mode must be 0, 1, or 2");
+                    "integer convolve/correlate mode must be 0, 1, 2, or 3");
             return NPY_FAIL;
         }
     }
diff --git a/numpy/core/src/multiarray/multiarraymodule.c b/numpy/core/src/multiarray/multiarraymodule.c
index 3a37a1db04b2..41fab7315072 100644
--- a/numpy/core/src/multiarray/multiarraymodule.c
+++ b/numpy/core/src/multiarray/multiarraymodule.c
@@ -3294,25 +3294,51 @@ array_correlate(PyObject *NPY_UNUSED(dummy),
             NULL, NULL, NULL) < 0) {
         return NULL;
     }
+    if (mode == 3) {
+        PyErr_SetString(PyExc_ValueError,
+                        "correlate() accepts only modes 0, 1, and 2 "
+                        "(valid, same and full). "
+                        "Use correlate2() for mode 3 (lags).");
+        return NULL;
+    }
     return PyArray_Correlate(a0, shape, mode);
 }
 
-static PyObject*
+static PyObject *
 array_correlate2(PyObject *NPY_UNUSED(dummy),
         PyObject *const *args, Py_ssize_t len_args, PyObject *kwnames)
 {
     PyObject *shape, *a0;
-    int mode = 0;
+    int mode = -1;
+    npy_intp minlag = 0, maxlag = 0, lagstep = 0;
+    static char *kwlist[] = {"a", "v", "mode", "minlag", "maxlag", "lagstep",
+                             NULL};
     NPY_PREPARE_ARGPARSER;
 
     if (npy_parse_arguments("correlate2", args, len_args, kwnames,
             "a", NULL, &a0,
             "v", NULL, &shape,
             "|mode", &PyArray_CorrelatemodeConverter, &mode,
+            "|minlag", NULL, &minlag,
+            "|maxlag", NULL, &maxlag,
+            "|lagstep", NULL, &lagstep,
             NULL, NULL, NULL) < 0) {
         return NULL;
     }
-    return PyArray_Correlate2(a0, shape, mode);
+    if (mode == -1) {
+        if (minlag == 0 && maxlag == 0 && lagstep == 0) {
+            /* if no lag parameters passed, use default: mode = 'valid' */
+            mode = 0;
+        }
+        else {
+            /* if lag parameters were passed, use them */
+            mode = 3;
+        }
+    }
+    if (mode != 3) {
+        return PyArray_Correlate2(a0, shape, mode);
+    }
+    return PyArray_CorrelateLags(a0, shape, minlag, maxlag, lagstep);
 }
 
 static PyObject *

From 9cfb40424fbe1b620badefffb3347c98627d078f Mon Sep 17 00:00:00 2001
From: Honi Sanders <sandersh6000@gmail.com>
Date: Mon, 22 May 2023 13:33:24 -0400
Subject: [PATCH 7/9] formatting of docstring for correlate

---
 doc/scipy-sphinx-theme | 1 +
 doc/sphinxext          | 1 +
 numpy/core/numeric.py  | 4 ++--
 3 files changed, 4 insertions(+), 2 deletions(-)
 create mode 160000 doc/scipy-sphinx-theme
 create mode 160000 doc/sphinxext

diff --git a/doc/scipy-sphinx-theme b/doc/scipy-sphinx-theme
new file mode 160000
index 000000000000..c466764e2231
--- /dev/null
+++ b/doc/scipy-sphinx-theme
@@ -0,0 +1 @@
+Subproject commit c466764e2231ba132c09826b5b138fffa1cfcec3
diff --git a/doc/sphinxext b/doc/sphinxext
new file mode 160000
index 000000000000..ef988a4a4658
--- /dev/null
+++ b/doc/sphinxext
@@ -0,0 +1 @@
+Subproject commit ef988a4a4658c991f4445f6241ab02d74710c6e3
diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index c1fbb2dcb0e7..e2416edbf478 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -877,8 +877,8 @@ def convolve(a, v, mode='full', lags=(), returns_lagvector=False):
           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.
+    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

From c1eed58c5b5c36b50b7f73165b8d46826108fc55 Mon Sep 17 00:00:00 2001
From: Yinnon Sanders <yinnons@gmail.com>
Date: Mon, 26 Jun 2023 11:08:38 -0400
Subject: [PATCH 8/9] Add lag arguments to dispatchers

---
 numpy/core/numeric.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index e2416edbf478..3d495e55b3d7 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -714,7 +714,7 @@ def _lags_from_mode(alen, vlen, mode):
     return mode_lags
 
 
-def _correlate_dispatcher(a, v, mode=None):
+def _correlate_dispatcher(a, v, mode=None, lags=None, returns_lagvector=None):
     return (a, v)
 
 
@@ -829,7 +829,7 @@ def correlate(a, v, mode='default', lags=(), returns_lagvector=False):
         return result
 
 
-def _convolve_dispatcher(a, v, mode=None):
+def _convolve_dispatcher(a, v, mode=None, lags=None, returns_lagvector=False):
     return (a, v)
 
 

From 05bce7de8b887d2c6df53c69f66dc8a1ca4849e7 Mon Sep 17 00:00:00 2001
From: Yinnon Sanders <yinnons@gmail.com>
Date: Mon, 26 Jun 2023 11:14:25 -0400
Subject: [PATCH 9/9] Change default lags argument to None

---
 numpy/core/numeric.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/numpy/core/numeric.py b/numpy/core/numeric.py
index 3d495e55b3d7..59115a57e0d8 100644
--- a/numpy/core/numeric.py
+++ b/numpy/core/numeric.py
@@ -719,7 +719,7 @@ def _correlate_dispatcher(a, v, mode=None, lags=None, returns_lagvector=None):
 
 
 @array_function_dispatch(_correlate_dispatcher)
-def correlate(a, v, mode='default', lags=(), returns_lagvector=False):
+def correlate(a, v, mode='default', lags=None, returns_lagvector=False):
     r"""
     Cross-correlation of two 1-dimensional sequences.
 
@@ -834,7 +834,7 @@ def _convolve_dispatcher(a, v, mode=None, lags=None, returns_lagvector=False):
 
 
 @array_function_dispatch(_convolve_dispatcher)
-def convolve(a, v, mode='full', lags=(), returns_lagvector=False):
+def convolve(a, v, mode='full', lags=None, returns_lagvector=False):
     """
     Returns the discrete, linear convolution of two one-dimensional sequences.