Repository: scipy/scipy
Pull request: DEP: integrate: Deprecate `initial` values other than None or 0 in cumulative_trapezoid
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Related to [#125](https://github.com/scipy/scipy/pull/125`)

#### What does this implement/fix?
<!--Please explain your changes.-->
The `initial` kw in `integrate.cumulative_trapezoid` was only intended as a way to make the output array have the same length as the input array. Ie: if initial=None, len(output) = len(input) - 1. If initial=0, len(output) = len(input), and the first value is 0.
But, `initial` can currently be any other scalar, and `cumulative_trapezoid` will prepend the value to the result. This is confusing because `initial` may be assumed to be an initial condition to the integration which is not true.

To remove the confusion this PR updates `cumulative_trapezoid` to raise a ValueError if `initial` is not None or 0.

#### Additional information
<!--Any additional information you think is important.-->
This issue has been discussed in #18151 and a Scipy mailing list [thread](https://mail.python.org/archives/list/scipy-dev@python.org/thread/5UYCLPVODQJRMMXCUQA5DNPBOIBKAJPK/).


Problem statement:
DOC: scipy.integrate.cumulative_trapezoid initial clarification
### Issue with current documentation:

Currently it is stated that the `initial` argument of `scipy.integrate.cumulative_trapezoid` is inserted as the resulting array of the cumulative integration.

However, I would have expected a cumulative function to add the `initial` to *all* values (since it is part of the cumulative function).

However, this is not clear from the documentation.

### Idea or request for content:

I can see that this may really break things, to I would suggest to add a small example clarifying the `initial` argument:

```python
from scipy.integrate import cumulative_trapezoid

print(cumulative_trapezoid([0, 1], dx=1, initial=2))
```
which yields `[2.   0.5]`

To get the initial in a cumulative fashion one should do:
```python
print(cumulative_trapezoid([0, 1], dx=1, initial=0)+2)
```
which yields `[2.   2.5]`

### Additional context (e.g. screenshots, GIFs)

_No response_

Hints from issue discussion:
Thanks @zerothi this will hopefully be fixed in #18219
Great! Thanks for this!

Commit message: STY: Remove trailing whitespace

Patch to review:
diff --git a/scipy/integrate/_quadrature.py b/scipy/integrate/_quadrature.py
index 1048e352675e..b54a69927fb9 100644
--- a/scipy/integrate/_quadrature.py
+++ b/scipy/integrate/_quadrature.py
@@ -411,6 +411,11 @@ def cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None):
         value at ``x[0]`` is returned and `res` has one element less than `y`
         along the axis of integration.
 
+        .. deprecated:: 1.12.0
+            The option for non-zero inputs for `initial` will be deprecated in
+            SciPy 1.14.0. After this time, a ValueError will be raised if
+            `initial` is not None or 0.
+
     Returns
     -------
     res : ndarray
@@ -472,6 +477,13 @@ def cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None):
     res = np.cumsum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)
 
     if initial is not None:
+        if initial != 0:
+            warnings.warn(
+                "The option for values for `initial` other than None or 0 is "
+                "deprecated as of SciPy 1.12.0 and will be removed in SciPy "
+                "1.14.0.",
+                DeprecationWarning, stacklevel=2
+            )
         if not np.isscalar(initial):
             raise ValueError("`initial` parameter should be a scalar.")
 
diff --git a/scipy/integrate/tests/test_quadrature.py b/scipy/integrate/tests/test_quadrature.py
index c9b713339b68..b2869eac152c 100644
--- a/scipy/integrate/tests/test_quadrature.py
+++ b/scipy/integrate/tests/test_quadrature.py
@@ -286,7 +286,7 @@ def test_y_nd_x_nd(self):
         # Try with all axes
         shapes = [(2, 2, 4), (3, 1, 4), (3, 2, 3)]
         for axis, shape in zip([0, 1, 2], shapes):
-            y_int = cumulative_trapezoid(y, x, initial=3.45, axis=axis)
+            y_int = cumulative_trapezoid(y, x, initial=0, axis=axis)
             assert_equal(y_int.shape, (3, 2, 4))
             y_int = cumulative_trapezoid(y, x, initial=None, axis=axis)
             assert_equal(y_int.shape, shape)
@@ -322,18 +322,28 @@ def test_x_none(self):
         y_expected = [-1.5, -2., -1.5, 0.]
         assert_allclose(y_int, y_expected)
 
-        y_int = cumulative_trapezoid(y, initial=1.23)
-        y_expected = [1.23, -1.5, -2., -1.5, 0.]
+        y_int = cumulative_trapezoid(y, initial=0)
+        y_expected = [0, -1.5, -2., -1.5, 0.]
         assert_allclose(y_int, y_expected)
 
         y_int = cumulative_trapezoid(y, dx=3)
         y_expected = [-4.5, -6., -4.5, 0.]
         assert_allclose(y_int, y_expected)
 
-        y_int = cumulative_trapezoid(y, dx=3, initial=1.23)
-        y_expected = [1.23, -4.5, -6., -4.5, 0.]
+        y_int = cumulative_trapezoid(y, dx=3, initial=0)
+        y_expected = [0, -4.5, -6., -4.5, 0.]
         assert_allclose(y_int, y_expected)
 
+    @pytest.mark.parametrize(
+        "initial", [1, 0.5]
+    )
+    def test_initial_warning(self, initial):
+        """If initial is not None or 0, a ValueError is raised."""
+        y = np.linspace(0, 10, num=10)
+        with pytest.deprecated_call():
+            res = cumulative_trapezoid(y, initial=initial)
+        assert_allclose(res, [initial, *np.cumsum(y[1:] + y[:-1])/2]) 
+
     def test_cumtrapz(self):
         # Basic coverage test for the alias
         x = np.arange(3 * 2 * 4).reshape(3, 2, 4)


Code review comments:
diff_hunk: @@ -322,18 +322,28 @@ def test_x_none(self):
         y_expected = [-1.5, -2., -1.5, 0.]
         assert_allclose(y_int, y_expected)
 
-        y_int = cumulative_trapezoid(y, initial=1.23)
-        y_expected = [1.23, -1.5, -2., -1.5, 0.]
+        y_int = cumulative_trapezoid(y, initial=0)
+        y_expected = [0, -1.5, -2., -1.5, 0.]
         assert_allclose(y_int, y_expected)
 
         y_int = cumulative_trapezoid(y, dx=3)
         y_expected = [-4.5, -6., -4.5, 0.]
         assert_allclose(y_int, y_expected)
 
-        y_int = cumulative_trapezoid(y, dx=3, initial=1.23)
-        y_expected = [1.23, -4.5, -6., -4.5, 0.]
+        y_int = cumulative_trapezoid(y, dx=3, initial=0)
+        y_expected = [0, -4.5, -6., -4.5, 0.]
         assert_allclose(y_int, y_expected)
 
+    @pytest.mark.parametrize(
+        "initial", [1, 0.5]
+    )
+    def test_initial_warning(self, initial):
+        """If initial is not None or 0, a ValueError is raised."""
+        y = np.linspace(0, 10, num=10)
+        with pytest.deprecated_call():
path: scipy/integrate/tests/test_quadrature.py
text: @user1:
```suggestion
        with pytest.deprecated_call(match="`initial`"):
```

@author:
Thanks @user2-bowhay , added in [aba400b](https://github.com/scipy/scipy/pull/18219/commits/aba400b702b2af0fec52a6607e6c547688cab6cd)