ENH: add cumulative_simpson integration to scipy.integrate

[nprav]

Reference issue

N/A

What does this implement/fix?

The cumulative integration and differentiation of numerical samples is a common operation in the engineering domains.
Scipy currently only provides the integrate.cumulative_trapezoid function, but this is insufficiently accurate for engineering purposes. The integrate.simpson method is available but only as a complete integral, not a cumulative integral. This PR adds an implementation for the cumulative integration of an array of samples using the Simpson's 1/3 rule (wiki link).

CAE softwares such as ANSYS/LS-DYNA provide simpson's integration as a default processing option. Adding this to Scipy.integrate will help expand usage in the engineering domains.

Additional information

Refer to this jupyter notebook: Cumulative Simpson derivation for further information.

[tupui]

Thank you for the PR @nprav!

Before we start any review, we need to get the approval for inclusion in SciPy. The process involves sending an email to the dev mailing list. Could you do that?

https://mail.python.org/mailman3/lists/scipy-dev.python.org/

[nprav]

Hi @tupui , I have emailed scipy-dev-owner@python.org and scipy-dev@python.org.
My email to scipy-dev@python.org got rejected since I still has a subscription request pending to be a list member.
Thanks!

[tupui]

@nprav thank you. The first time you send an email to the list, an admin needs to validate. Once this is done your email should go through.

[tupui]

@MatteoRaso Thank you for starting a review but please wait until there is a decision from a maintainer. We don't want people to spend time on something that could in the end be rejected.

EDIT: we now got an approval 😃

[lorentzenchr]

Some review comments. IMHO, it is already looking quite good and a very useful addition to scipy.

[andyfaff]

A comment from the sidelines. I would expect such a function to have the same functionality as the existing integrate.simpson function, and for its behaviour to eventually reduce to the existing function.

i.e. it should handle end points similarly, and the behaviour of overlapping keywords should be the same, etc.

[nprav]

Thank you for the comment @andyfaff ! This is an important discussion point because there are currently some differences in the function signatures. Please find the differences and justifications below and let me know what you think.

Current function signatures:

• simpson(y, x=None, dx=1.0, axis=-1, even='avg')
• cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None)
• cumulative_simpson(y, x=None, dx=1.0, axis=-1, initial=None)

Comparisons

Integrate.simpson(y, x=None, dx=1.0, axis=-1, even="avg")

• y, x, dx, axis work the same way
• even is not used in cumulative_simpson

The 'even' keyword relates to what happens if there are an odd number of intervals. Integrate.simpson will necessarily apply a trapezoidal rule for the odd interval. However, this approach violates the assumption of a quadratic relationship. In reality, the final interval can be approximated with a modified formula while maintaining the quadratic assumption. This is described in the wikipedia article.

It is also depicted with figures here: Cumulative Simpson Derivation.

Because of the above, I had implemented cumulative_simpson without the trapezoidal rule logic. This also ensures that when applied to samples with a quadratic relationship, cumulative_simpson will always give the exact solution. On the other hand, integrate.simpson will be inaccurate if there are an odd number of intervals.

I also compared accuracy when integrating oscillatory functions and found cumulative_simpson to be exhibit equal or better accuracy than integrate.simpson (above notebook link).

cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None)

• y, x, dx, axis work the same way.
• Initial is slightly different
  • In cumulative_trapezoid, initial is simply prepended to the result array.
  • In cumulative_simpson, initial is prepended to the result array AND summed to every value in the result

The second approach is more 'correct' because the 'initial' value is analogous to a constant of integration that must be added to every point in the cumulative integral. Consider the following practical example:

• Let y be sampled velocity data of an oscillating object. We wish to cumulatively integrate to attain the displacement of the object. We know the starting position (initial value) is 10m. Without the initial value input, the displacement data will show oscillations around 0m. With the initial value, the displacement should show oscillations around 10m.
• But, with the cumulative_trapezoid implementation, the first value would be 10m, and the object would then be oscillating around 0m, which is incorrect. The summing approach would fix this.

[andyfaff]

I don't have time to go through this PR myself (plus I'm not an expert in this area of the codebase). However, I should reiterate that behaviour for the proposed PR should tally exactly with existing behaviour for simpson and cumulative_trapezoid.

For example:

-Initial is slightly different
-In cumulative_trapezoid, initial is simply prepended to the result array.
-In cumulative_simpson, initial is prepended to the result array AND summed to every value in the result

The cumulative_simpson would not keep with existing behaviour of cumulative_trapezoid because it adds initial to every value in the result. This would create confusion with users because of the differing behaviour. See #125 for why cumulative_trapezoid has that behaviour.

Similarly I would expect the exact same behaviour for simpson + cumulative_simpson, down to the behaviour of the even keyword being included, which is not what you outlined here. Different behaviour leads to confusion.

However, I do see the benefit of using the Cartwright2017 approach, as you linked to on Wikipedia. It may be that there's benefit to changing the existing code in simpson to take that into account, but that would have to be discussed more thoroughly as it has to take into account back-compatibility.

[mdhaber]

Crossref gh-5618 and gh-9871. Please consider linking gh-5618.

[andyfaff]

@mdhaber I'ce prepped a fix for #5618 in https://github.com/andyfaff/scipy/tree/simpson. The mechanics are there but the API needs a bit of discussion. E.g.. how do we deal with back compat, do we ignore even kwd, add an extra option to even, what's that option called, is that option default, etc.

[mdhaber]

Yup, your recent email reminded me of those issues.

[nprav]

Thanks @nprav. I think we'll finish this before we branch. I just want to make sure it's all expressed in a way that can easily be maintained. It looks like you've been working on this for a long time and have gone through several reviewers, so I don't want to push your patience. If there are other things (like the input validation and corresponding tests) that are technically correct but can be expressed more clearly (it seems that you agreed that the changes were improvements), are you willing to help with that?

Yes @mdhaber, I'm willing to help with any specific comments.

[mdhaber]

Here are some thoughts on the actual implementation to make the code easier to verify and improve performance.

[mdhaber]

FAILED scipy/integrate/tests/test_quadrature.py::TestCumulativeSimpson::test_simpson_exceptions[`axis=4` is not valid for `y` with `y.ndim=1`-kwarg_update4] - AttributeError: module 'numpy' has no attribute 'AxisError'

Not sure why this is failing on just one platform. Maybe it does need to be numpy.exceptions.AxisError, not just numpy.AxisError?

Maybe except IndexError (which it subclasses) if we can't use AxisError.

[mdhaber]

After #18151 (review) and nprav#2 are considered, I think this is good to go. Only merge them if you agree; LMK either way. Thanks for your patience @nprav!

[mdhaber]

There are some details to improve, but let's get this in. Thanks @nprav!

[nprav]

Thank you all for helping to add this to SciPy! @mdhaber @tupui @andyfaff @lucascolley @lorentzenchr @MatteoRaso @j-bowhay

[lucascolley]

thanks @nprav and congratulations on your first contribution to SciPy! Keep them coming :)

[mdhaber]

@nprav Here are some things to consider for follow-up work. If you are willing to do this soon, we can get it in the version of the function that is released in SciPy 1.12.

There are other things about documentation that I'd change in all the related function, but I've avoided commenting about those sorts of things here if they were just copied over.

[nprav]

@mdhaber , should i add the follow up changes in a separate PR?

[mdhaber]

Yes, thanks!

[nprav]

@mdhaber
Working on these comments in: https://github.com/nprav/scipy/tree/cumulative_simpson-clean-up
Will add a PR to the maintenance branch soon.

Regarding the case of "... elements of x identical"
If all elements of x are identical, it is easy: integral is 0.

But if only some elements are identical, it is tricky. Eg: x = [0, 0, 1, 1, 2, 2]
This causes ZeroDivisionErrors in the internal functions.
In Simpson, this is treated by just skipping the ZeroDivisionError points and just taking them as 0: https://github.com/scipy/scipy/blob/main/scipy/integrate/_quadrature.py#L560

We can do the same in cumulative_simpson, but the results may not match with simpson because cumulative_simpson explicitly calculates the intermediate points. More importantly, this does not make much mathematical sense. I think it should be up to the user to determine how to clean/handle their inputs. At the same time, I understand we want to keep the API consistent.

Do you have any thoughts on this?

[lucascolley]

Will add a PR to the maintenance branch soon.

All PRs should be to main, then the release manager can backport.

[nprav]

@mdhaber
Added #19709 for the follow-up comments.