search

hgrecco / numbakit-ode

Leveraging numba to speed up ODE integration
Other
68 stars 3 forks source link

Support for events? #13

Closed astrojuanlu closed 3 years ago

astrojuanlu commented 3 years ago

(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 https://github.com/hgrecco/numbakit-ode/issues/11#issuecomment-752248440

hgrecco commented 3 years ago

I have a work in progress that is promising. The proposed API is the following:

  1. Keep events as defined in the SciPy API (i.e. a function or an iterable of functions with optional direction and is_terminal attributes, events occur when the function cross zero) Reason: this API has proven to work and keeping it will simplify porting code.
  2. 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.
  3. Provide a new method (run_events?) which is like run but with an extra parameters (events) Reason: do not use run as the output will be different. (see point 4)
  4. Output: I see three possible outputs: a. t, y, t_events, y_events: where t, y are the normal result of run and t_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 where t, y are the normal result of run and t_events, y_events, n_events are the times, states and event source) c. t_all, y_all, n_all where t_all, y_all, n_all is interleaving t, y from run and t_events, y_events, and adding a fake 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.

maurosilber commented 3 years ago

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.

hgrecco commented 3 years ago

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.

astrojuanlu commented 3 years ago

That was fast, thanks a lot! I'll take a look... next year 😄 Feliz 2021 desde el otro lado del charco 🥂

hgrecco commented 3 years ago

A few comments:

  1. Version 4a is implemented
  2. 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 when run (and not run_events) is called Proposal: defer this discussion until the method events API is stabilized and in use
  3. Not sure if implementing skip_events is necessary as run_events(np.inf, events) will do exactly that. What will be the use of step_events?
hgrecco commented 3 years ago

Merged. Docs are still missing.