ENH: integrate.solve_ivp: allow event `terminal` attribute to specify maximum number of events #14771

Patol75 · 2021-09-27T00:01:33Z

Reference issue

What does this implement/fix?

The terminal attribute of an event function as part of solve_ivp is deprecated in favour of the new max_events attribute. Instead of asserting if integration should terminate on an event occurrence, max_events allows specifying the number of occurrences of a given event after which integration must terminate. Through setting event.max_events = 1, the behaviour of event.terminal = True is recovered. Additionally, using event.max_events = np.inf or, more simply, not setting max_events will ensure that integration never terminates on the occurrence of a given event, similar to writing event.terminal = False. Finally, the enhancement lies in the fact that event.max_events can be given any value in ]1, np.inf[ and integration will only terminate if the number of occurrences of the specified event becomes greater or equal to the latter value.

Additional information

Original concept by @drhagen in gh-14176.

This sounds like a good feature to me. One thing to note is that the current terminal behavior is equivalent to terminating after 1 event, so the most appropriate API would be deprecating terminal as an attribute on the event and replacing it with max_events or some similar name with a default of inf.

This would be a somewhat disruptive deprecation, but the deprecation period could be very long because it is not particularly bad to have these features side-by-side. In the interim, I would propose that setting both terminal and max_events would be an error. Setting terminal=True would be equivalent to max_events=1.

scipy/integrate/_ivp/tests/test_ivp.py

scipy/integrate/_ivp/ivp.py

tylerjereddy · 2021-09-28T02:30:42Z

Probably SciPy 1.10 at the earliest for removal if this functionality is public/heavily used, but maybe wait for an integrate regular to confirm. assert_allclose() is defined to work on array_like, which in the NumPy glossary is defined as:

Any scalar or sequence that can be interpreted as an ndarray. In addition to ndarrays and scalars this category includes lists (possibly nested and with different element types) and tuples. Any argument accepted by numpy.array is array_like.

Patol75 · 2021-09-28T05:51:11Z

Thanks, I have updated the test file, happy to have your feedback about it.

drhagen · 2021-10-02T19:12:53Z

This implementation looks good to me. I'll defer to Tyler on the deprecation cycle.

Patol75 · 2021-10-19T23:26:34Z

Hi @tylerjereddy,

Considering @drhagen is supportive of the new implementation and that CI is green except for a failure on CircleCI that I cannot manage to reproduce locally and I believe is unrelated to the changes proposed here, do you think we can move towards merging this PR if I update the SciPy target version to 1.10? Thanks for your input.

Patol75 · 2021-11-16T03:27:42Z

It looks like the DeprecationWarning is now considered as an error by the CI; is that the expected behaviour?

mdhaber · 2022-04-16T15:02:31Z

Yes, please use pytest.warns to catch the deprecation warning.

…r to allow for integration termination to occur only after the event was triggered a given amount of times.

…the 'terminate_after' attribute.

…onal check within test_events to ensure multiple events can be captured before termination.

…d with deprecation; Update tests.

mdhaber · 2022-04-26T15:02:34Z

@patol what do you think about just adding a callback function, since that has also been requested, and that could take care of verbose and max_events at the same time without complicating the interface or needing deprecation?
Also, I'd be happy to review and merge a callback PR, but I (personally) would prefer not to review three separate PRs for these when we could get all three features with one addition to the API.

Patol75 · 2022-04-26T17:22:33Z

I see the idea. You would have access to the solution object after every integration step within the callback function and be able to interact with it as you wish. However, regarding the max_events functionality, you would still need most of the changes included in the present PR. Otherwise, there would not be any counter for event occurrences. Another thing to consider is how user-friendly such a "powerful" callback would be. It may not be obvious to many users how to deal with it, and it would definitely be more difficult than, for example, setting verbose=True.

mdhaber · 2022-04-26T17:33:15Z

The idea was to give examples of how to use callback to do the same job as verbose and max_events in the documentation. That way, they have turnkey solutions for those simple things, but they also learn a bit so they can customize them for their own needs.

Otherwise, there would not be any counter for event occurrence

For that, I thought callback would be a method of an object that remembered the counts, or it could be a function that increments a counter declared outside its scope.

I agree that, say, toggling verbose is less work and thought for the user than implementing a callback. My concern is that solve_ivp already has ten arguments, one of them being an options structure with several other components. It's getting quite complicated, so I thought we should be more conservative with additions. (And I realized this only after I did the work to finish that verbose PR gh-11815!)

Patol75 · 2023-03-02T00:56:36Z

@mdhaber Sorry, I never got to properly consider implementing a callback feature. Do you still believe it should be the way forward? Or should this PR actually go ahead?

j-bowhay

Can you add a test that checks the deprecation warning is produced is terminal is used (if you decide to go ahead)

j-bowhay · 2023-03-02T08:43:59Z

scipy/integrate/_ivp/ivp.py

+            terminal = event.terminal
+            warnings.warn(
+                "The 'terminal' event function attribute is deprecated and "
+                "will be removed in SciPy release 1.11.0. Please use the new "


Suggested change

"will be removed in SciPy release 1.11.0. Please use the new "

"will be removed in SciPy release 1.13.0. Please use the new "

mdhaber · 2023-03-02T22:37:36Z

Do you still believe it should be the way forward?

I'm open to other opinions, but I've seen a lot of specialized feature requests in integrate and optimize that could be resolved if the relevant functions had a callback interface, so I think that a callback interface would be more efficient in terms of developer/reviewer time and API.

Patol75 · 2023-03-04T04:08:51Z

Do you still believe it should be the way forward?

I'm open to other opinions, but I've seen a lot of specialized feature requests in integrate and optimize that could be resolved if the relevant functions had a callback interface, so I think that a callback interface would be more efficient in terms of developer/reviewer time and API.

I think I will be willing to look into implementing this callback functionality. However, it would have to wait a little while, as I am quite busy at the moment. What I am realising is that, with such a functionality, one would have as much control over the solver evolution as one currently has by manually calling, for example, RK45 or LSODA. Would that make the solver classes redundant? If yes, should they become part of the private API?

Continuing on the callback feature, I perfectly see how one would easily replicate the proposed verbose addition to solve_ivp, but I still struggle to see how the same can be said about the changes proposed here. To stop the solver execution, solve_ivp relies on the terminal attribute. On the one hand, if we do not change that part of the API, then how would you allow for multiple occurrences of the same event? On the other hand, if we retire the direction and terminal attributes and it becomes up to the user to keep track of the number of event occurrences, would it then also be up to the user to interrupt the solver execution through the callback function? I am not sure this is super desirable...

mdhaber · 2023-03-04T16:49:57Z

Would that make the solver classes redundant?

Maybe. (But this would not be a reason to withhold such a capability from the function interface.)

If yes, should they become part of the private API?

Sometimes I would prefer that we have only one interface instead of two (e.g. root_scalar instead of root_scalar and brent), especially when there is only enough difference between the two to make them incompatible but too little to make each situationally useful. ATM, I think solve_ivp and ODESolver are different enough and each has advantages. Perhaps it could have been argued to release only one to begin with, but I'm not motivated to go through the trouble of deprecating one that already exists. I suspect this would not be popular.

but I still struggle to see how the same can be said about the changes proposed here.

When I posted that, I thought that the callback function would be informed of events, and it could keep track of the number of events. It could set the terminal attribute to True before the final event or (probably preferably) be capable of halting integration by raising a StopIteration exception or similar. This is something we're now adding throughout optimize.

I am not sure this is super desirable..

Why not? In optimize, for example, users have asked for exotic termination criteria. Rather than extending the API ad infinitum, we can refer them to callback.

zachjweiner · 2023-03-04T19:39:03Z

I might be able to help with a callback feature (although I am also pretty swamped in the near term). I once hacked such an interface by overriding _step_impl to allow manual control over step acceptance. Here's a gist of the implementation - obviously this is specific to my use case (which was to trigger adjusting approximation schemes in the system and then redo the step) and not generic. Just sharing in case it's a useful starting point.

drhagen · 2023-03-18T15:42:13Z

I'm open to other opinions, but I've seen a lot of specialized feature requests in integrate and optimize that could be resolved if the relevant functions had a callback interface, so I think that a callback interface would be more efficient in terms of developer/reviewer time and API.

The main purpose of exposing the object-oriented interface solver interface (the step solvers) was to make it easy for users to implement whatever specialized features they wanted. solve_ivp was a convenience function to expose the most common features in way that was familiar to users of Matlab's ODE solvers and the solvers in other languages inspired by them (like SciPy's own odeint). I thought that having the step solvers would lessen the need to add more and more specialized features to solve_ivp, but that does not seem to be true. I don't know if people prefer the callback style because of unfamiliarity or because they'd like it even if they knew how to use the step solvers.

mdhaber · 2023-03-19T19:43:42Z

I don't know if people prefer the callback style because of unfamiliarity or because they'd like it even if they knew how to use the step solvers.

Familiarity is probably a big part of it. In addition to what you wrote about Matlab's ODE solvers and odeint:

I have personal experience being taught and teaching function-based ODE solvers in schools. I haven't seen OO options taught in academia, personally.
There are many other functions in SciPy (and beyond) that accept a callable and use it to solve a problem iteratively (e.g. in scipy.optimize). There could be OO options for these (and we're discussing an OO version of scipy.optimize.linprog) but users are more familiar with the more common function-based versions. After experiences with these, it would be natural to choose a function-based ODE solver.

Even if one knows about the OO option, solve_ivp is easier to get started with. (It's a convenience function, after all.) Unlike the OO options, it has an example in the API documentation, a section in the scipy.integrate tutorial, and the code is more compact, which could all make beginners lean toward it instead of the OO interface. When their needs become more advanced, why learn something new or restructure existing code when you can just add a callback?

That said, I wouldn't object to a consensus that solve_ivp is only for the most common use-cases and the OO interface is for the rest. Any unified solution to the many separate enhancement requests will probably sound good to me.

jkitchin · 2023-06-08T12:34:57Z

Would a callback allow something like this to work?

A recent problem I have been working on is detecting steady state cyclic behavior in an ODE solution. One way to do this is to keep track of some solution property in different time windows, e.g. we can integrate over different time windows, and if those are constant within a tolerance, terminate. Right now the only solution seems to be solve over some interval, do these analyses, and if they are not satisfied then restart the integrator at the last position.

Its not currently possible to accumulate events from an event function because it is called many times in the root finding code to locate the event.

j-bowhay · 2023-06-08T12:50:56Z

Would a callback allow something like this to work?

A recent problem I have been working on is detecting steady state cyclic behavior in an ODE solution. One way to do this is to keep track of some solution property in different time windows, e.g. we can integrate over different time windows, and if those are constant within a tolerance, terminate. Right now the only solution seems to be solve over some interval, do these analyses, and if they are not satisfied then restart the integrator at the last position.

Its not currently possible to accumulate events from an event function because it is called many times in the root finding code to locate the event.

As much as it pains me to say this as a scipy maintainer, I would recommend checking out BifurcationKit.jl. It has implementations of the shooting, colocation and trapezoid methods for detecting periodic orbits/limit cycles in dynamical systems and is really quite nice to use. https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/periodicOrbit/

tylerjereddy · 2023-12-04T22:43:43Z

It looks like @mdhaber is suggesting consideration of an alternative approach that uses a callback design, and I don't see a completely unified agreement about things in the discussion just yet, so I'd be hesitant to merge this.

I'll bump the milestone to 1.13.0; one potentially useful thing if there is demand for a solution here is that 1.13.0 should also come out fairly soon as it aims to support NumPy 2.0.0 rather than just being six months after 1.12.0.

tylerjereddy · 2024-03-12T17:07:16Z

Still unresolved, bumping milestone again ahead of 1.13.0 branching.

mdhaber · 2024-03-14T20:42:49Z

@j-bowhay Since it doesn't look like callback will happen soon, I'm thinking about resolving the conflicts and merging this since the feature has been needed for a while now. To make it as simple as possible, what do you think about 1) just accepting an integer rather than just a boolean for terminal and 2) not refactoring all the existing tests? Then, 1) renaming terminal to max_events or similar and 2) refactoring the tests can be separate PRs if desired.

j-bowhay · 2024-03-15T14:36:55Z

@j-bowhay Since it doesn't look like callback will happen soon, I'm thinking about resolving the conflicts and merging this since the feature has been needed for a while now. To make it as simple as possible, what do you think about 1) just accepting an integer rather than just a boolean for terminal and 2) not refactoring all the existing tests? Then, 1) renaming terminal to max_events or similar and 2) refactoring the tests can be separate PRs if desired.

👍 Sound much easier to review!

mdhaber · 2024-03-15T18:37:15Z

scipy/integrate/_ivp/ivp.py

@@ -1,5 +1,6 @@
 import inspect
 import numpy as np
+import warnings


Suggested change

import warnings

mdhaber · 2024-03-15T18:41:45Z

@j-bowhay would you be willing to take a look at this before 1.13 branches? I simplified prepare_events, left the logic in handle_events alone, updated the documentation, and added a simple test; others can be added (e.g. multiple events, use of direction) if desired, but I didn't think they were needed. Regarding input validation, I opted for simplicity above all else: non-positive integers are treated as 0, non-integral floats are just rounded down, and invalid values will get whatever message int(event.terminal) provides.

j-bowhay

Core logic lgtm, one query and one nonblocking suggestion

j-bowhay · 2024-03-15T19:16:24Z

scipy/integrate/_ivp/ivp.py

+    direction = np.empty(len(events))
+    for i, event in enumerate(events):
+        terminal = int(event.terminal) if hasattr(event, 'terminal') else np.inf
+        max_events[i] = terminal if terminal > 0 else np.inf


Any reason for not raising an error here for negative values? If I passed nonsense I think I'd rather know about it rather than silently receive some other behaviour

When I write new functions, I am pretty careful about this sort of thing. When I work on old functions - especially for minor work like this - I consider how it will affect the overall quality of the function's input validation and try to balance that against simplicity of the PR. Currently, for example:

If event.direction=None, no events are detected.

If event.terminal = 'no', then the function terminates after the first event.

No errors are raised and these behavior are surprising, so I didn't think treating negative values like 0 made the function's edges substantially rougher.

However, seeing now that the current behavior for event.terminal = -1 is also to cause the function to terminate after the first event, I agree that we should raise the error instead of changing that behavior. I'll add the error.

Oops, I'll also replace the ternary and hasattr with getattr.

scipy/integrate/_ivp/ivp.py

mdhaber · 2024-03-16T04:02:56Z

Thanks @j-bowhay. Fixed that and added a test.

j-bowhay

Thanks @Patol75 @mdhaber

… maximum number of events (scipy#14771) * ENH: The terminal attribute of an event function can now be an integer to allow for integration termination to occur only after the event was triggered a given amount of times. * ENH: Deprecate the 'terminal' event function attribute and introduce the 'terminate_after' attribute. * ENH: Change attribute name terminate_after to max_events; TST: Additional check within test_events to ensure multiple events can be captured before termination. * MAINT: Rebase on top of current main branch; Update release associated with deprecation; Update tests. * MAINT: integrate.solve_ivp: simplify prepare_events * MAINT: integreate.solve_ivp: update documentation, add test * MAINT: integrate.solve_ivp: add input validation to attribute --------- Co-authored-by: Matt Haberland <mhaberla@calpoly.edu>

github-actions bot added the scipy.integrate label Sep 27, 2021

tylerjereddy added the enhancement A new feature or improvement label Sep 27, 2021

tylerjereddy reviewed Sep 27, 2021

View reviewed changes

scipy/integrate/_ivp/tests/test_ivp.py Outdated Show resolved Hide resolved

scipy/integrate/_ivp/ivp.py Outdated Show resolved Hide resolved

This was referenced Apr 16, 2022

ENH: Add max_events parameter for flexible termination of scipy.solve_ivp #14395

Closed

ENH: integrate.solve_ivp: add verbose as optional argument to print time step info #11815

Open

Patol75 added 4 commits April 26, 2022 12:55

ENH: The terminal attribute of an event function can now be an intege…

db54492

…r to allow for integration termination to occur only after the event was triggered a given amount of times.

ENH: Deprecate the 'terminal' event function attribute and introduce …

92d0a65

…the 'terminate_after' attribute.

ENH: Change attribute name terminate_after to max_events; TST: Additi…

12d9610

…onal check within test_events to ensure multiple events can be captured before termination.

MAINT: Rebase on top of current main branch; Update release associate…

c5e63ce

…d with deprecation; Update tests.

Patol75 force-pushed the Max_Events_Solve_IVP branch from b3b80fd to c5e63ce Compare April 26, 2022 14:09

j-bowhay requested changes Mar 2, 2023

View reviewed changes

j-bowhay added this to the 1.11.0 milestone Mar 2, 2023

h-vetinari modified the milestones: 1.11.0, 1.12.0 May 20, 2023

tylerjereddy modified the milestones: 1.12.0, 1.13.0 Dec 4, 2023

tylerjereddy modified the milestones: 1.13.0, 1.14.0 Mar 12, 2024

mdhaber added 3 commits March 15, 2024 10:50

Merge remote-tracking branch 'upstream/main' into Max_Events_Solve_IVP

0e42e73

MAINT: integrate.solve_ivp: simplify prepare_events

224c301

MAINT: integreate.solve_ivp: update documentation, add test

d8bd33b

Patol75 requested a review from steppi as a code owner March 15, 2024 18:36

mdhaber reviewed Mar 15, 2024

View reviewed changes

mdhaber changed the title ~~ENH: Replace terminal attribute by max_events in solve_ivp~~ ENH: integrate.solve_ivp: allow event terminal attribute to specify maximum number of events Mar 15, 2024

mdhaber requested review from j-bowhay and removed request for steppi March 15, 2024 18:51

j-bowhay reviewed Mar 15, 2024

View reviewed changes

MAINT: integrate.solve_ivp: add input validation to attribute

c58d0b7

mdhaber modified the milestones: 1.14.0, 1.13.0 Mar 16, 2024

mdhaber approved these changes Mar 16, 2024

View reviewed changes

j-bowhay approved these changes Mar 16, 2024

View reviewed changes

j-bowhay merged commit 8e5b96a into scipy:main Mar 16, 2024
29 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: integrate.solve_ivp: allow event `terminal` attribute to specify maximum number of events #14771

ENH: integrate.solve_ivp: allow event `terminal` attribute to specify maximum number of events #14771

Patol75 commented Sep 27, 2021 •

edited

tylerjereddy commented Sep 28, 2021

Patol75 commented Sep 28, 2021

drhagen commented Oct 2, 2021

Patol75 commented Oct 19, 2021

Patol75 commented Nov 16, 2021

mdhaber commented Apr 16, 2022

mdhaber commented Apr 26, 2022 •

edited

Patol75 commented Apr 26, 2022

mdhaber commented Apr 26, 2022 •

edited

Patol75 commented Mar 2, 2023

j-bowhay left a comment •

edited

j-bowhay Mar 2, 2023

mdhaber commented Mar 2, 2023

Patol75 commented Mar 4, 2023

mdhaber commented Mar 4, 2023 •

edited

zachjweiner commented Mar 4, 2023

drhagen commented Mar 18, 2023

mdhaber commented Mar 19, 2023

jkitchin commented Jun 8, 2023

j-bowhay commented Jun 8, 2023

tylerjereddy commented Dec 4, 2023

tylerjereddy commented Mar 12, 2024

mdhaber commented Mar 14, 2024 •

edited

j-bowhay commented Mar 15, 2024

mdhaber Mar 15, 2024

mdhaber commented Mar 15, 2024 •

edited

j-bowhay left a comment

j-bowhay Mar 15, 2024 •

edited

mdhaber Mar 15, 2024 •

edited

mdhaber Mar 15, 2024

mdhaber commented Mar 16, 2024

j-bowhay left a comment

	"will be removed in SciPy release 1.11.0. Please use the new "
	"will be removed in SciPy release 1.13.0. Please use the new "

ENH: integrate.solve_ivp: allow event terminal attribute to specify maximum number of events #14771

ENH: integrate.solve_ivp: allow event terminal attribute to specify maximum number of events #14771

Conversation

Patol75 commented Sep 27, 2021 • edited

Reference issue

What does this implement/fix?

Additional information

tylerjereddy commented Sep 28, 2021

Patol75 commented Sep 28, 2021

drhagen commented Oct 2, 2021

Patol75 commented Oct 19, 2021

Patol75 commented Nov 16, 2021

mdhaber commented Apr 16, 2022

mdhaber commented Apr 26, 2022 • edited

Patol75 commented Apr 26, 2022

mdhaber commented Apr 26, 2022 • edited

Patol75 commented Mar 2, 2023

j-bowhay left a comment • edited

Choose a reason for hiding this comment

j-bowhay Mar 2, 2023

Choose a reason for hiding this comment

mdhaber commented Mar 2, 2023

Patol75 commented Mar 4, 2023

mdhaber commented Mar 4, 2023 • edited

zachjweiner commented Mar 4, 2023

drhagen commented Mar 18, 2023

mdhaber commented Mar 19, 2023

jkitchin commented Jun 8, 2023

j-bowhay commented Jun 8, 2023

tylerjereddy commented Dec 4, 2023

tylerjereddy commented Mar 12, 2024

mdhaber commented Mar 14, 2024 • edited

j-bowhay commented Mar 15, 2024

mdhaber Mar 15, 2024

Choose a reason for hiding this comment

mdhaber commented Mar 15, 2024 • edited

j-bowhay left a comment

Choose a reason for hiding this comment

j-bowhay Mar 15, 2024 • edited

Choose a reason for hiding this comment

mdhaber Mar 15, 2024 • edited

Choose a reason for hiding this comment

mdhaber Mar 15, 2024

Choose a reason for hiding this comment

mdhaber commented Mar 16, 2024

j-bowhay left a comment

Choose a reason for hiding this comment

ENH: integrate.solve_ivp: allow event `terminal` attribute to specify maximum number of events #14771

ENH: integrate.solve_ivp: allow event `terminal` attribute to specify maximum number of events #14771

Patol75 commented Sep 27, 2021 •

edited

mdhaber commented Apr 26, 2022 •

edited

mdhaber commented Apr 26, 2022 •

edited

j-bowhay left a comment •

edited

mdhaber commented Mar 4, 2023 •

edited

mdhaber commented Mar 14, 2024 •

edited

mdhaber commented Mar 15, 2024 •

edited

j-bowhay Mar 15, 2024 •

edited

mdhaber Mar 15, 2024 •

edited