Research/document desired functionality of Python (scipy) integrators

[nicklafarge]

We need to understand what functionality from [scipy.integrate.solve_ivp] (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) we want to see incorporated into the python API

[rjpower4]

"New" API

• solve_ivp
• RK23
• RK45
• DOP853
• Radau
• BDF
• LSODA
• OdeSolver
• DenseOutput
• OdeSolution

"Old" API

• odeint
• ode
• complex_ode

[DhruvJ22]

We should incorporate most of the features into our Python API. I have strikethrough the ones that are overkill/unnecessary to include for our purpose and added some comments for things that we can add.

Parameters:

1. t_span: integration interval (to, tf)
2. y0 {only states, no Identity even if STM is to be computed}
3. Integration method {we can discuss which methods to include for v 1.0 in Select which Boost ODE methods to initially support #3}
4. t_eval: specify @t_i at which to store solutions
5. dense_output
6. events:

• predefine common events functions like y=0, x = 1-mu and make them accessible using a string name
• option to define a "custom" events function in python and pass that to the underneath c++
• option for terminal and direction
• option to get the n^th, or an array of i (e.g. [1,2,5,6]) events and terminate after the last event specified.

7. rtol, atol
8. vectorized
9. firsrt_step
10. max_step
11. jac
12. jac_sparsity
13. lband, uband
14. min_step

Returns: {we should try to keep the "Returns" API similar to solve_ivp to make it easier for users to integrate forest in their existing work}

1. t (ndarray)
2. y (values of states at t)
3. sol
4. t_events (list of ndarray)
5. y_events (list, values of states at t)
6. n_events (new parameter, list of # event numbers corresponding to the value of t_events; e.g. [1,2,5,6])
7. nfev
8. njev
9. nlu
10. status (int, reason for algorithm termination)
11. message (string)
12. success (bool)

[noahsadaka]

I think that rather than have several specific axis-crossing event functions, it would be more flexible to have a single event function where you can specify an axis and a value along that axis for the hyperplane to exist at. So setting axis y and value 0 would give our typical x-axis crossing function, but users wouldn't be constrained to using only that event function or whatever other specific axis-crossing event function that is included in forest.

[rjpower4]

Axis-crossing events are a special type of the more general set of affine events of the form

$$ \vec{f}\left(\vec{x}\right) = A\vec{x} + \vec{b} $$

with $\vec{x} \in \mathbb{R}^n$, $\vec{f}, \vec{b} \in \mathbb{R}^m$, $A \in \mathbb{R}^{m \times n}$. So I think that implementation should encompass most of these along. We can provide convenience methods for plane crossing events which are just cases where $A$ is a one-hot vector.

[nicklafarge]

Agreed on all points above.

Also, I think we also should think carefully before we expose too much Python functionality (@rjpower4 and I have discussed this previously). For example, Python-defined EOMs sound convenient and fun in theory, but ultimately may defeat the point of a rapid integration package if it suddenly requires going back-and-forth to Python each step of the integrator. The convenience of Python-defined functionality may also make people less likely to implement their model/event/whatever directly in C++ (which may then benefit the user community in general). The same could be potentially said for event functions.

Not saying we shouldn't do any of those things, but it is a perspective to consider when deciding where to draw that line. For events specifically, RJ's affine event suggestion sounds great.

[DhruvJ22]

I agree that we should not expose EOMs, that is why I did not explicitly include 'fun' as an input.

I think it will be a good option to allow users to pass A and $\bar{b}$ on the Python side to define an events function. With this, the users can retain the ease of defining straightforward event functions, and I don't think it will significantly affect the C++ performance. I see it as constant variables like abs_tol that are set at the beginning of the integration loop. We can push the users to define more complex event functions, like periapsis and altitude, directly in C++, which can then be used by other community members.

[nicklafarge]

Sounds good! Thanks for documenting all this @DhruvJ22 !

[nicklafarge]

This discussion provides some good insight for #13 as well