Support for events?
Closed this issue · 6 comments
(Comes from #11)
Event are trickier but it would be nice to have also for the original purpose of the numbakit-ode (biological simulations). It would be nice to define the API first.
Originally posted by @hgrecco in #11 (comment)
I have a work in progress that is promising. The proposed API is the following:
- Keep events as defined in the SciPy API (i.e. a function or an iterable of functions with optional
direction
andis_terminal
attributes, events occur when the function cross zero)
Reason: this API has proven to work and keeping it will simplify porting code. - Provide events at method calling, not at class initialization.
Reason: flexibly as it will allow running the same system with different set of events in different parts. - Provide a new method (
run_events
?) which is likerun
but with an extra parameters (events
)
Reason: do not userun
as the output will be different. (see point 4) - Output: I see three possible outputs:
a.t, y, t_events, y_events
: wheret, y
are the normal result ofrun
andt_events, y_events
mimic the SciPy API output (i.e. a pair list of ndarrays, the nth element of the each list corresponds to then nth event function)
b.t, y, t_events, y_events, n_events
wheret, y
are the normal result ofrun
andt_events, y_events, n_events
are the times, states and event source)
c.t_all, y_all, n_all
wheret_all, y_all, n_all
is interleavingt, y
from run andt_events, y_events
, and adding afake
source for n_events (value = -1) that corresponds to non event (i.e times added because the user requested when calling the method)
I see advantages and disadvantages for all options of 4.
2- How about providing them at both places? It would support continuing to run the same solver instance without having to explicitly pass events at each run
call. But we'd have to decide what to do if there are instance and call events simultaneously. Should them override or extend the instance events?
4- I think I prefer version a
.
Maybe it helps to think of events in the other methods (skip
and step
) too. For example, when computing a Poincaré map, we might be only interested in the events, without saving the full trajectory. Would that be a skip_events
that returns t_events
and y_events
?
I would exclude c
, so as not to mix the values given by the steps, with the values interpolated for events. Thinking in terms of dense output, we would need to either get the t
and y
where n_all == -1
to perform the interpolation, or keep a separate copy of the t
and y
given by the steps.
I have added a PR to develop with an initial implementation. I have chosen a jitclass based implementation as I think it makes things easier to read and understand. I have not seen a performance hit as compared to using just functions and passing data around.
I would love to hear your thoughts.
That was fast, thanks a lot! I'll take a look... next year 😄 Feliz 2021 desde el otro lado del charco 🥂
A few comments:
- Version 4a is implemented
- My problem with providing events at the instance level are:
a. How to deal when they are also provided at the method level (as mentioned)
b. What to do whenrun
(and notrun_events
) is called
Proposal: defer this discussion until the method events API is stabilized and in use - Not sure if implementing
skip_events
is necessary asrun_events(np.inf, events)
will do exactly that. What will be the use ofstep_events
?
Merged. Docs are still missing.