New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: New ODE solvers #6326

Merged
merged 94 commits into from Apr 21, 2017

Conversation

Projects
None yet
@nmayorov
Contributor

nmayorov commented Jun 28, 2016

Hi!

First some considerations about the current ode class.

  1. It allows to use several good quality solvers.
  2. It feels like a giant black box, the usability and docs are not the best. It's hard to understand if you get your setup right or not, seeing output from fortran code is awkward, etc.
  3. It's not clear how to efficiently compute the solution at many points ("dense output"). For Runge-Kutta methods (dopri5 and dop853) just calling integrate by small steps would be very wrong, but that's what a user would probably do (a "proper" way is discussed here) . For other solvers the situation should be better, but still it is not efficient to make small steps (but this is limitation of solvers and not on scipy side).
  4. No continuous solution output, i.e. which can be evaluated at any point with the same accuracy as the values computed at discrete points.
  5. No event detection capabilities.

So I want to provide a new monolithic implementation, which is easy to understand and use and also provides missing features 4 and 5. Yes, it will be in Python and not as fast. But the same situation is in MATLAB and people are totally fine with it, as the solvers are well designed and easy to use. In fact I was largely motivated by MATLAB suite. I realize that this is not a super clear situation, and there are pros and cons. But comparing the current state of scipy and MATLAB in ODE integration, I come to conclusion that scipy could do better.

Currently I only implemented two explicit Runge-Kutta methods (admittedly the easiest methods, but no doubt useful). And getting right implicit solvers for stiff problems will be a challenge.

I prepared an example with two main features, continuous output and event detection: http://nbviewer.jupyter.org/gist/nmayorov/ca99381bb20a87762758f7dad4fc2220

The code is in preliminary stage. The current TODO list:

  • Implement at least one implicit method for stiff problems.
  • Refine API: add necessary parameters, rename current parameters, etc.
  • [ x ] Develop a benchmark suite and assess quality of the algorithms.

Any feedback and suggestions are very welcome. A concrete technical question: which implicit method you believe would be the most useful (if implement only one)? The BDF (and its modifications) seems to be the most widely used, so I want to start looking into it.

@charris

This comment has been minimized.

Show comment
Hide comment
@charris

charris Jun 28, 2016

Member

@aarchiba IIRC, you had an interest in some of the topics, e.g., event detection.

Member

charris commented Jun 28, 2016

@aarchiba IIRC, you had an interest in some of the topics, e.g., event detection.

@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Jun 28, 2016

Member

Didn't look at this yet, but cross-linking some previous (inconclusive) discussion on the ivp api redesign:

Regarding pure-python implementation --- I think the Python overhead is more significant for ivps that bvps as the latter can be vectorized. For ivps you can, and often want to, write the evaluation function in a compiled language.

Member

pv commented Jun 28, 2016

Didn't look at this yet, but cross-linking some previous (inconclusive) discussion on the ivp api redesign:

Regarding pure-python implementation --- I think the Python overhead is more significant for ivps that bvps as the latter can be vectorized. For ivps you can, and often want to, write the evaluation function in a compiled language.

@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Jun 28, 2016

Member

What I would like to avoid here is a situation where we end up with three different interfaces to different solvers, written by different people, each providing 80% solutions with differing tradeoffs. So I would like to see that the new interface covers all the use cases of the previous ones.

Member

pv commented Jun 28, 2016

What I would like to avoid here is a situation where we end up with three different interfaces to different solvers, written by different people, each providing 80% solutions with differing tradeoffs. So I would like to see that the new interface covers all the use cases of the previous ones.

@aarchiba

This comment has been minimized.

Show comment
Hide comment
@aarchiba

aarchiba Jun 28, 2016

Contributor

Returning a PPoly object turns out not to be the best way to handle dense output. I thought they might be helpful for this, which is why I implemented them, but they actually produce rather poor results (perhaps because of the Runge phenomenon). Regardless, most modern ODE solvers implement dense output as part of the algorithm: when you advance the solver by one step (possibly adaptively-chosen), the solver retains some state. This state is used to compute the position at the end of the step, but a dense output algorithm also allows evaluation within that last step, using that state, and provides a solution of guaranteed order and accuracy (usually the same order as the method, sometimes one lower). These algorithms usually also permit evaluation of the derivative (without RHS calls) albeit at a lower order of accuracy.

The interface to most solvers, including VODE in scipy itself, do not make it easy to keep this state around. This is partly because they're in FORTRAN, but partly also because in many applications the full freely-evaluable solution object is not needed. For example, event detection can be carried out by root-finding within the last step. Nevertheless, under most circumstances, the state usually has size something like the method order (maybe squared) times the number of dimensions; this need not be expensive to keep around.

Contributor

aarchiba commented Jun 28, 2016

Returning a PPoly object turns out not to be the best way to handle dense output. I thought they might be helpful for this, which is why I implemented them, but they actually produce rather poor results (perhaps because of the Runge phenomenon). Regardless, most modern ODE solvers implement dense output as part of the algorithm: when you advance the solver by one step (possibly adaptively-chosen), the solver retains some state. This state is used to compute the position at the end of the step, but a dense output algorithm also allows evaluation within that last step, using that state, and provides a solution of guaranteed order and accuracy (usually the same order as the method, sometimes one lower). These algorithms usually also permit evaluation of the derivative (without RHS calls) albeit at a lower order of accuracy.

The interface to most solvers, including VODE in scipy itself, do not make it easy to keep this state around. This is partly because they're in FORTRAN, but partly also because in many applications the full freely-evaluable solution object is not needed. For example, event detection can be carried out by root-finding within the last step. Nevertheless, under most circumstances, the state usually has size something like the method order (maybe squared) times the number of dimensions; this need not be expensive to keep around.

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Jun 28, 2016

Contributor

@aarchiba I want to clarify that I don't compute PPoly like a Lagrange polynomial passing through all computed points or anything like that (this would be very wrong). But actually I use well thought mathematical procedures for computing the polynomial of the same order of accuracy (well, 1 order less for RK45 but it's OK because of the global error) as the integration method ("as part of the algorithm").You can check a reference from the docstring or "Solving Ordinary Differential Equation I, p. 191". For example, a cubic Hermit polynomial will be good for any method with order <= 3 and it is used for 'RK23'.

Contributor

nmayorov commented Jun 28, 2016

@aarchiba I want to clarify that I don't compute PPoly like a Lagrange polynomial passing through all computed points or anything like that (this would be very wrong). But actually I use well thought mathematical procedures for computing the polynomial of the same order of accuracy (well, 1 order less for RK45 but it's OK because of the global error) as the integration method ("as part of the algorithm").You can check a reference from the docstring or "Solving Ordinary Differential Equation I, p. 191". For example, a cubic Hermit polynomial will be good for any method with order <= 3 and it is used for 'RK23'.

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Jun 28, 2016

Contributor

@pv

  1. Totally agree that IVP will suffer more from overhead. My guess is that it will be really bad only if a solver makes too much steps (like 10000+), but this really should not happen if a problem is properly formulated and an appropriate (and well implemented) solver is used. If the problem dimension is big, then numpy's linear algebra will kick in and save the situation. It's only theory though, but I think it will be fine for interactive usage (sure not for some real time processing). Rewriting in Cython is also an option, right?
  2. I agree that it becomes too many different interfaces (and it is one of disadvantageous I thought of). I think one of the old interfaces should be removed and one retained. So a user can select if he wants efficiency or convenience. But it's too early to decide.
Contributor

nmayorov commented Jun 28, 2016

@pv

  1. Totally agree that IVP will suffer more from overhead. My guess is that it will be really bad only if a solver makes too much steps (like 10000+), but this really should not happen if a problem is properly formulated and an appropriate (and well implemented) solver is used. If the problem dimension is big, then numpy's linear algebra will kick in and save the situation. It's only theory though, but I think it will be fine for interactive usage (sure not for some real time processing). Rewriting in Cython is also an option, right?
  2. I agree that it becomes too many different interfaces (and it is one of disadvantageous I thought of). I think one of the old interfaces should be removed and one retained. So a user can select if he wants efficiency or convenience. But it's too early to decide.

@rgommers rgommers referenced this pull request Jun 28, 2016

Closed

Odeint rewrite #2818

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jun 28, 2016

Member

There's so much history here and several alternative proposals (this one, the one from @Juanlu001 in gh-2818, one from @bmcage in the odes scikit here. Plus a bunch of discussion/ideas on the mailing list, moving to Sundials, etc. I have to admit I can't quite fit that in my head - a PEP-style summary of the desired features and implementation alternatives may be very useful.

Member

rgommers commented Jun 28, 2016

There's so much history here and several alternative proposals (this one, the one from @Juanlu001 in gh-2818, one from @bmcage in the odes scikit here. Plus a bunch of discussion/ideas on the mailing list, moving to Sundials, etc. I have to admit I can't quite fit that in my head - a PEP-style summary of the desired features and implementation alternatives may be very useful.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jun 28, 2016

Member

I prepared an example with two main features, continuous output and event detection: http://nbviewer.jupyter.org/gist/nmayorov/ca99381bb20a87762758f7dad4fc2220

@nmayorov very nice examples!

Member

rgommers commented Jun 28, 2016

I prepared an example with two main features, continuous output and event detection: http://nbviewer.jupyter.org/gist/nmayorov/ca99381bb20a87762758f7dad4fc2220

@nmayorov very nice examples!

@aarchiba

This comment has been minimized.

Show comment
Hide comment
@aarchiba

aarchiba Jun 28, 2016

Contributor

@nmayorov yes the examples are lovely! I don't have easy access to the references right now, but I think you are confusing the degree of the polynomial with the order of the approximation. In particular, for (say) Dormand-Prince 853, the value at the end of a step of length h converges to the right value like h^9, but with your quadric polynomials it's not even guaranteed to converge as the degree of the polynomial. You really need to use the interior points (some of which are not on the solution curve) to construct an approximation of the solution curve, and you need to know something about the guts of the algorithm to do this. The dopri FORTRAN implementations do just this.

Contributor

aarchiba commented Jun 28, 2016

@nmayorov yes the examples are lovely! I don't have easy access to the references right now, but I think you are confusing the degree of the polynomial with the order of the approximation. In particular, for (say) Dormand-Prince 853, the value at the end of a step of length h converges to the right value like h^9, but with your quadric polynomials it's not even guaranteed to converge as the degree of the polynomial. You really need to use the interior points (some of which are not on the solution curve) to construct an approximation of the solution curve, and you need to know something about the guts of the algorithm to do this. The dopri FORTRAN implementations do just this.

@aarchiba

This comment has been minimized.

Show comment
Hide comment
@aarchiba

aarchiba Jun 28, 2016

Contributor

@rgommers I put forward just such a design on the mailing list ages ago: https://mail.scipy.org/pipermail/scipy-dev/2015-June/020739.html

It doesn't cover this use case - anywhere-evaluable solution objects - but this functionality is easily implemented on top of it.

Contributor

aarchiba commented Jun 28, 2016

@rgommers I put forward just such a design on the mailing list ages ago: https://mail.scipy.org/pipermail/scipy-dev/2015-June/020739.html

It doesn't cover this use case - anywhere-evaluable solution objects - but this functionality is easily implemented on top of it.

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Jun 28, 2016

Contributor

Thank you for complimenting my examples :) They are simple really.

@aarchiba

I probably should have explained everything straight away rather than pointing to references.

You really need to use the interior point.

And I do that for RK45 (which we also can call dopri45), I compute y(x + 0.5 * h) using a special linear combination of 7 Runge-Kutta stages (including f(x + h)) and then construct a quartic polynomial using 3 values and 2 derivatives. This approach is suggested as an alternative to construction of a 5-th order polynomial (which nevertheless only 4-th order accurate) as done in DOPRI5. I decided to stick with a 4-th order polynomial mainly because 4 < 5. Also I could verify it in sympy as there is a linear system provided to compute these coefficients (in fact I used slightly non-optimal value of c6 to avoid 8 digit numbers). And then construction of a quartic polynomial is a trivial matter.

As for 853 method --- I didn't consider it. I assume that 45 will be enough (I even considering removing 23 method).

Thank you for the feedback, I hope to hear more.

@rgommers

The basic idea is to make a new transparent implementation in Python (maybe later Cythonize some parts) with a convenient and unified interface. I want it to be user friendly and feel very seamless, which is hard to achieve by wrapping different existing codes. I think it will benefit scipy greatly as a research/education environment, provided the solvers will be algorithmically correct and efficient. Run time in pure Python is an issue, but I'm not very pessimistic about that.

It's far from being ready. I guess the main criterion is if it will be complete and efficient enough. For now I will continue working on it, it should get more clear is it worth including.

Contributor

nmayorov commented Jun 28, 2016

Thank you for complimenting my examples :) They are simple really.

@aarchiba

I probably should have explained everything straight away rather than pointing to references.

You really need to use the interior point.

And I do that for RK45 (which we also can call dopri45), I compute y(x + 0.5 * h) using a special linear combination of 7 Runge-Kutta stages (including f(x + h)) and then construct a quartic polynomial using 3 values and 2 derivatives. This approach is suggested as an alternative to construction of a 5-th order polynomial (which nevertheless only 4-th order accurate) as done in DOPRI5. I decided to stick with a 4-th order polynomial mainly because 4 < 5. Also I could verify it in sympy as there is a linear system provided to compute these coefficients (in fact I used slightly non-optimal value of c6 to avoid 8 digit numbers). And then construction of a quartic polynomial is a trivial matter.

As for 853 method --- I didn't consider it. I assume that 45 will be enough (I even considering removing 23 method).

Thank you for the feedback, I hope to hear more.

@rgommers

The basic idea is to make a new transparent implementation in Python (maybe later Cythonize some parts) with a convenient and unified interface. I want it to be user friendly and feel very seamless, which is hard to achieve by wrapping different existing codes. I think it will benefit scipy greatly as a research/education environment, provided the solvers will be algorithmically correct and efficient. Run time in pure Python is an issue, but I'm not very pessimistic about that.

It's far from being ready. I guess the main criterion is if it will be complete and efficient enough. For now I will continue working on it, it should get more clear is it worth including.

@aarchiba

This comment has been minimized.

Show comment
Hide comment
@aarchiba

aarchiba Jun 28, 2016

Contributor

@nmayorov They say in HNW that using a Hermite-type interpolant is enough (and easier) for low orders, but for more elaborate methods (Bulirsch-Stoer, dopri853, or even this one) you have to poke around inside the guts to get out the additional points and derivatives you need. Further, we are going to want to use established codes like CVODE that do more-careful error analysis and stiff detection, so we're going to need to use their dense evaluation functions. Fortunately all the codes I'd consider using provide dense evaluation! So while the PPoly code may be enough for this particular implementation, we're going to need something that can take advantage of this feature of underlying solvers.

I think that we can in fact provide one or more convenient seamless interfaces atop existing codes; what we need from the underlying solvers is limited and specific, so writing adapters to fit them into a well-designed user-facing layer need not be difficult. I'm sketching out what I'd want from such a user-facing layer and will post something about it in a bit.

Contributor

aarchiba commented Jun 28, 2016

@nmayorov They say in HNW that using a Hermite-type interpolant is enough (and easier) for low orders, but for more elaborate methods (Bulirsch-Stoer, dopri853, or even this one) you have to poke around inside the guts to get out the additional points and derivatives you need. Further, we are going to want to use established codes like CVODE that do more-careful error analysis and stiff detection, so we're going to need to use their dense evaluation functions. Fortunately all the codes I'd consider using provide dense evaluation! So while the PPoly code may be enough for this particular implementation, we're going to need something that can take advantage of this feature of underlying solvers.

I think that we can in fact provide one or more convenient seamless interfaces atop existing codes; what we need from the underlying solvers is limited and specific, so writing adapters to fit them into a well-designed user-facing layer need not be difficult. I'm sketching out what I'd want from such a user-facing layer and will post something about it in a bit.

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Jun 30, 2016

Member

@rgommers I put forward just such a design on the mailing list ages ago: https://mail.scipy.org/pipermail/scipy-dev/2015-June/020739.html

Thanks @aarchiba, that helps!

Member

rgommers commented Jun 30, 2016

@rgommers I put forward just such a design on the mailing list ages ago: https://mail.scipy.org/pipermail/scipy-dev/2015-June/020739.html

Thanks @aarchiba, that helps!

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Jul 11, 2016

Contributor

I implemented the Radau implicit solver. I updated my examples with demonstration that it works for stiff problems http://nbviewer.jupyter.org/gist/nmayorov/ca99381bb20a87762758f7dad4fc2220

For now I don't provide a support for the mass matrix. Here is my considerations why:

  1. If allow singular matrices (which is supported by the Radau method) then the solver will be able to solve index 1 implicit DAEs (formulated as M y' = f(x, y)). While this is good, it will complicate the status of this function. a) More explanations will be necessary. b) In this case consistent initial conditions are required and if a user will fail to provide those, the result will be to some extent incorrect or strange (I'm not sure now how exactly it will behave). Ideally the solver should correct the initial conditions in this case treating them only as a guess, but this is another non-trivial addition.
  2. If allow only constant and non-singular matrix then the value of this addition won't be very high (easy to handle on a user's side).
  3. If allow time and state dependent mass matrix then the algorithm will get more complicated. If mass matrix is state dependent then it will be necessary to estimate its derivatives by finite differences (or ask to provide another Jacobian).

So it's not easy to make a good support of the mass matrix, currently I'm not sure if I want to do it.

Next I want to focus on the quality of the algorithms in terms of number of steps and correctness of a solution on a set of test problems.

Does anyone want to run them on his/her problems and see if it works satisfactory?

Contributor

nmayorov commented Jul 11, 2016

I implemented the Radau implicit solver. I updated my examples with demonstration that it works for stiff problems http://nbviewer.jupyter.org/gist/nmayorov/ca99381bb20a87762758f7dad4fc2220

For now I don't provide a support for the mass matrix. Here is my considerations why:

  1. If allow singular matrices (which is supported by the Radau method) then the solver will be able to solve index 1 implicit DAEs (formulated as M y' = f(x, y)). While this is good, it will complicate the status of this function. a) More explanations will be necessary. b) In this case consistent initial conditions are required and if a user will fail to provide those, the result will be to some extent incorrect or strange (I'm not sure now how exactly it will behave). Ideally the solver should correct the initial conditions in this case treating them only as a guess, but this is another non-trivial addition.
  2. If allow only constant and non-singular matrix then the value of this addition won't be very high (easy to handle on a user's side).
  3. If allow time and state dependent mass matrix then the algorithm will get more complicated. If mass matrix is state dependent then it will be necessary to estimate its derivatives by finite differences (or ask to provide another Jacobian).

So it's not easy to make a good support of the mass matrix, currently I'm not sure if I want to do it.

Next I want to focus on the quality of the algorithms in terms of number of steps and correctness of a solution on a set of test problems.

Does anyone want to run them on his/her problems and see if it works satisfactory?

@drhagen

This comment has been minimized.

Show comment
Hide comment
@drhagen

drhagen Aug 12, 2016

Contributor

I found this pull request while looking for a way to do ODE event detection in scipy. I think this is a great effort. I am willing to provide help as I am able. I am quite familiar with Matlab's ODE solvers.

I noticed that the type of object returned as x_events is dependent on the number of events provided. This will be problematic for code that accepts a variable number of events. Example: Say I write a function count_events that takes a list of events, runs an ODE with those events, and returns the number of each event found. I may think that I simply need to loop over x_events and return the size of each ndarray, and this will work fine until I send a list of one event, and then it will fail because the type of output depends on the value of the input.

For now I don't provide a support for the mass matrix.

I think this is the right thing to do. While there is mathematical overlap between DAEs and ODEs, the concepts and the interfaces are quite different.

I agree that it becomes too many different interfaces (and it is one of disadvantageous I thought of). I think one of the old interfaces should be removed and one retained.

Ideally, there would be (1) an interface for stepping through the integration one step at a time, (2) a function that gets the value of an ODE system at a set of predefined times, and (3) a function that gets the dense solution until a final time.

Contributor

drhagen commented Aug 12, 2016

I found this pull request while looking for a way to do ODE event detection in scipy. I think this is a great effort. I am willing to provide help as I am able. I am quite familiar with Matlab's ODE solvers.

I noticed that the type of object returned as x_events is dependent on the number of events provided. This will be problematic for code that accepts a variable number of events. Example: Say I write a function count_events that takes a list of events, runs an ODE with those events, and returns the number of each event found. I may think that I simply need to loop over x_events and return the size of each ndarray, and this will work fine until I send a list of one event, and then it will fail because the type of output depends on the value of the input.

For now I don't provide a support for the mass matrix.

I think this is the right thing to do. While there is mathematical overlap between DAEs and ODEs, the concepts and the interfaces are quite different.

I agree that it becomes too many different interfaces (and it is one of disadvantageous I thought of). I think one of the old interfaces should be removed and one retained.

Ideally, there would be (1) an interface for stepping through the integration one step at a time, (2) a function that gets the value of an ODE system at a set of predefined times, and (3) a function that gets the dense solution until a final time.

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Aug 12, 2016

Contributor

@drhagen, it's great that you have an interest in this work. It should give me a motivation to actually continue and finish it. Any help (including discussion) will be very useful.

Noted about different types of x_events. Looks like it is seldom a good idea to use different types of return values.

Ideally, there would be (1) an interface for stepping through the integration one step at a time, (2) a function that gets the value of an ODE system at a set of predefined times, and (3) a function that gets the dense solution until a final time.

My understanding that (2) is solved in terms of (3). And (1), could you explain how it should work and what benefits it should bring?


How would you suggest to continue this PR? Here are some thoughts and questions:

  1. What integration methods you believe are essential? I'm sure that implemented explicit RK methods are enough for non-stiff problems. For stiff problems there a lot of methods besides Radau. But if implemented Radau will be proved to work well, then there are good reasons not to add anything else.
  2. After the algorithms are established, it should be beneficial to translate the code into Cython. Even though there are likely will be many Python interactions ("yellow lines" in cython -a), I think it will become considerably faster than MATLAB for large number of steps.
  3. As I wrote before, I believe now it's time to do thorough algorithm's testing. If the results are good and API / features are established we can do Cython translation and be done with it.
Contributor

nmayorov commented Aug 12, 2016

@drhagen, it's great that you have an interest in this work. It should give me a motivation to actually continue and finish it. Any help (including discussion) will be very useful.

Noted about different types of x_events. Looks like it is seldom a good idea to use different types of return values.

Ideally, there would be (1) an interface for stepping through the integration one step at a time, (2) a function that gets the value of an ODE system at a set of predefined times, and (3) a function that gets the dense solution until a final time.

My understanding that (2) is solved in terms of (3). And (1), could you explain how it should work and what benefits it should bring?


How would you suggest to continue this PR? Here are some thoughts and questions:

  1. What integration methods you believe are essential? I'm sure that implemented explicit RK methods are enough for non-stiff problems. For stiff problems there a lot of methods besides Radau. But if implemented Radau will be proved to work well, then there are good reasons not to add anything else.
  2. After the algorithms are established, it should be beneficial to translate the code into Cython. Even though there are likely will be many Python interactions ("yellow lines" in cython -a), I think it will become considerably faster than MATLAB for large number of steps.
  3. As I wrote before, I believe now it's time to do thorough algorithm's testing. If the results are good and API / features are established we can do Cython translation and be done with it.
@drhagen

This comment has been minimized.

Show comment
Hide comment
@drhagen

drhagen Aug 12, 2016

Contributor

Noted about different types of x_events. Looks like it is seldom a good idea to use different types of return values.

The Julia people call this "type stability". It gives a name to something I don't like about Matlab, which does this a lot, a vector of length 1 or 2 is treated differently from a longer vector. It would be fine if a Callable events returned an ndarray x_events and a List[Callable] events returned a List[ndarray] x_events.

My understanding that (2) is solved in terms of (3).

Matlab has both (2) and (3). In my biological modeling toolbox written in Matlab, I tried to simplify the simulation code by always integrating the full solution and down sampling when the user requested only particular points, but ultimately had to make two parallel code paths because of the impact on performance. I learned that there is a substantial impact on performance if I stored the dense solution and then sampled later versus keeping only the desired samples along the way. This becomes particularly important when the system is large and stiff. The integrator may need to take a huge number of tiny steps in rapidly changing regions that the user doesn't ultimately care about. It is computationally expensive to allocate memory for all that only to throw it away.

That being said leaving (2) for later or implementing (2) in terms of (3) is fine for now. Performance can be improved once the interface is working well.

And (1), could you explain how it should work and what benefits it should bring?

This is essentially the scipy.integrate.ode class. This is one of the features that I really like about scipy, even though the current implementation could use a little sprucing up. Both (2) and (3) can easily be solved with (1). The advantage of having an interface like this is that you can write anything else with fairly easily. Say a user wants to do something unusual with an ODE solver, like have a function called at every step, or apply a function to decide when to stop the solver, or make an object that can return the solution at time 0 to inf by extending the solution when necessary.

What integration methods you believe are essential?

This is very problem dependent. I am going test my systems on Radau. I use ode15s in Matlab, and the little bit of interaction I've had with SUNDIALS has made suspect that CVODE is the best that there is. But I don't know the source for the underlying algorithms. I only work with very stiff systems (biology), so I have no experience with good solvers for non-stiff systems.

After the algorithms are established, it should be beneficial to translate the code into Cython

I think this is a great idea. Get the interface perfect and then find the slow spots and speed them up with whatever tools are available.

As I wrote before, I believe now it's time to do thorough algorithm's testing.

I have a Windows 7 machine and have been unable to get scipy to compile today. Where is the best place to get help with this? The scipy mailing list? Stackoverflow?

Contributor

drhagen commented Aug 12, 2016

Noted about different types of x_events. Looks like it is seldom a good idea to use different types of return values.

The Julia people call this "type stability". It gives a name to something I don't like about Matlab, which does this a lot, a vector of length 1 or 2 is treated differently from a longer vector. It would be fine if a Callable events returned an ndarray x_events and a List[Callable] events returned a List[ndarray] x_events.

My understanding that (2) is solved in terms of (3).

Matlab has both (2) and (3). In my biological modeling toolbox written in Matlab, I tried to simplify the simulation code by always integrating the full solution and down sampling when the user requested only particular points, but ultimately had to make two parallel code paths because of the impact on performance. I learned that there is a substantial impact on performance if I stored the dense solution and then sampled later versus keeping only the desired samples along the way. This becomes particularly important when the system is large and stiff. The integrator may need to take a huge number of tiny steps in rapidly changing regions that the user doesn't ultimately care about. It is computationally expensive to allocate memory for all that only to throw it away.

That being said leaving (2) for later or implementing (2) in terms of (3) is fine for now. Performance can be improved once the interface is working well.

And (1), could you explain how it should work and what benefits it should bring?

This is essentially the scipy.integrate.ode class. This is one of the features that I really like about scipy, even though the current implementation could use a little sprucing up. Both (2) and (3) can easily be solved with (1). The advantage of having an interface like this is that you can write anything else with fairly easily. Say a user wants to do something unusual with an ODE solver, like have a function called at every step, or apply a function to decide when to stop the solver, or make an object that can return the solution at time 0 to inf by extending the solution when necessary.

What integration methods you believe are essential?

This is very problem dependent. I am going test my systems on Radau. I use ode15s in Matlab, and the little bit of interaction I've had with SUNDIALS has made suspect that CVODE is the best that there is. But I don't know the source for the underlying algorithms. I only work with very stiff systems (biology), so I have no experience with good solvers for non-stiff systems.

After the algorithms are established, it should be beneficial to translate the code into Cython

I think this is a great idea. Get the interface perfect and then find the slow spots and speed them up with whatever tools are available.

As I wrote before, I believe now it's time to do thorough algorithm's testing.

I have a Windows 7 machine and have been unable to get scipy to compile today. Where is the best place to get help with this? The scipy mailing list? Stackoverflow?

@rgommers

This comment has been minimized.

Show comment
Hide comment
@rgommers

rgommers Aug 13, 2016

Member

I have a Windows 7 machine and have been unable to get scipy to compile today. Where is the best place to get help with this? The scipy mailing list? Stackoverflow?

The scipy-user mailing list is a good place, or open an issue if you think it's a bug. Note that there's an existing Windows issue with the latest setuptools under Windows / Python 3.5. If you have that combination, downgrading to setuptools 23 may solve it.

Member

rgommers commented Aug 13, 2016

I have a Windows 7 machine and have been unable to get scipy to compile today. Where is the best place to get help with this? The scipy mailing list? Stackoverflow?

The scipy-user mailing list is a good place, or open an issue if you think it's a bug. Note that there's an existing Windows issue with the latest setuptools under Windows / Python 3.5. If you have that combination, downgrading to setuptools 23 may solve it.

@aarchiba

This comment has been minimized.

Show comment
Hide comment
@aarchiba

aarchiba Aug 14, 2016

Contributor

I don't have any experience with stiff ODEs, but I have spent a lot of time on solving a very smooth problem to very high accuracy. I used Bulirsch-Stoer, and it is generally agreed to be a useful method for such a situation; I have a partially-written implementation intended for inclusion in scipy, along with three levels of interface to it:

(1) A self-contained class that takes one (adaptive) step at a time and allows evaluation (including derivatives) anywhere within the last step, optionally saving the state for later dense evaluation,
(2) A function that returns a solution object evaluable anywhere (to either side of the start point), in any order, and of any derivative up to the order of the method, and
(3) A function that simply returns the solution at a prescribed set of points.

New solver methods just need to be available to (1), which (2) and (3) use. Similarly, event detection can be implemented on top of (1) independently of solution methods. Ideally we could use generic root-finding code shared with scipy.optimize, but achieving compiled efficiency when calling back and forth is nontrivial.

I'm not sure how to allow users to control finer details of event detection - high-order methods like Bulirsch-Stoer may take quite large steps, possibly large enough to miss events if the event function is non-monotonic. Should event functions have access to derivatives of the solution? How easy is it to vectorize event detection - that is, detect when some function of the solution passes each of many values? (In my problem I had a "time" quantity that was a monotonically increasing function of the ODE time variable, and I wanted to know the values of another function when it passed each of thirty thousand values.)

Efficiency - handling compiled code - is also a major concern. To be fully useful, we need the solver, dense evaluation, and event detection to run in compiled code, and we also need users to be able to supply compiled right-hand sides or event functions and have them evaluated without passing through the interpreter.

Parallelization is another challenge. I don't think trying to squeeze parallel execution out of solving individual initial-value problems, but careful management of the GIL may allow parallelization of solving for many initial conditions. This won't work for FORTRAN solvers or right-hand sides, but there should be C or Cython solvers that make this possible.

Contributor

aarchiba commented Aug 14, 2016

I don't have any experience with stiff ODEs, but I have spent a lot of time on solving a very smooth problem to very high accuracy. I used Bulirsch-Stoer, and it is generally agreed to be a useful method for such a situation; I have a partially-written implementation intended for inclusion in scipy, along with three levels of interface to it:

(1) A self-contained class that takes one (adaptive) step at a time and allows evaluation (including derivatives) anywhere within the last step, optionally saving the state for later dense evaluation,
(2) A function that returns a solution object evaluable anywhere (to either side of the start point), in any order, and of any derivative up to the order of the method, and
(3) A function that simply returns the solution at a prescribed set of points.

New solver methods just need to be available to (1), which (2) and (3) use. Similarly, event detection can be implemented on top of (1) independently of solution methods. Ideally we could use generic root-finding code shared with scipy.optimize, but achieving compiled efficiency when calling back and forth is nontrivial.

I'm not sure how to allow users to control finer details of event detection - high-order methods like Bulirsch-Stoer may take quite large steps, possibly large enough to miss events if the event function is non-monotonic. Should event functions have access to derivatives of the solution? How easy is it to vectorize event detection - that is, detect when some function of the solution passes each of many values? (In my problem I had a "time" quantity that was a monotonically increasing function of the ODE time variable, and I wanted to know the values of another function when it passed each of thirty thousand values.)

Efficiency - handling compiled code - is also a major concern. To be fully useful, we need the solver, dense evaluation, and event detection to run in compiled code, and we also need users to be able to supply compiled right-hand sides or event functions and have them evaluated without passing through the interpreter.

Parallelization is another challenge. I don't think trying to squeeze parallel execution out of solving individual initial-value problems, but careful management of the GIL may allow parallelization of solving for many initial conditions. This won't work for FORTRAN solvers or right-hand sides, but there should be C or Cython solvers that make this possible.

@drhagen

This comment has been minimized.

Show comment
Hide comment
@drhagen

drhagen Aug 15, 2016

Contributor

I'm not sure how to allow users to control finer details of event detection - high-order methods like Bulirsch-Stoer may take quite large steps, possibly large enough to miss events if the event function is non-monotonic

These would be events that are changing faster than the states they are tracking, am I right? I've never seen methods for dealing with things like this. This would seem to be a good argument for having the event detection being built on top of the step-by-step solver rather than into it so you could swap out the standard detector for your own method.

Contributor

drhagen commented Aug 15, 2016

I'm not sure how to allow users to control finer details of event detection - high-order methods like Bulirsch-Stoer may take quite large steps, possibly large enough to miss events if the event function is non-monotonic

These would be events that are changing faster than the states they are tracking, am I right? I've never seen methods for dealing with things like this. This would seem to be a good argument for having the event detection being built on top of the step-by-step solver rather than into it so you could swap out the standard detector for your own method.

@aarchiba

This comment has been minimized.

Show comment
Hide comment
@aarchiba

aarchiba Aug 15, 2016

Contributor

I'm not sure how to allow users to control finer details of event detection - high-order methods like Bulirsch-Stoer may take quite large steps, possibly large enough to miss events if the event function is non-monotonic

These would be events that are changing faster than the states they are tracking, am I right? I've never seen methods for dealing with things like this. This would seem to be a good argument for having the event detection being built on top of the step-by-step solver rather than into it so you could swap out the standard detector for your own method.

Well, there are two situations I'm concerned about:

  • The hypothetical situation of using event detection on a non-monotonic function, for example, using an ODE solver to get sin(t) and event detection to find where it crosses a level y. Here the concern is that the solver might take big steps - B-S can manage good accuracy for this problem using step sizes of order unity - and you might miss an event because the function crosses the threshold twice. Here the user might need to specify a maximum step size, or better, a brute-force sub-sampling of the solver step for use in event detection. Derivative information might help here too, or the fact that for most solver methods the dense output actually comes from evaluating a polynomial.
  • The practical situation is that I have a very smooth integration problem where I need to know about the system state at each of a prescribed list of "times", where each "time" is computed from the system state. The "times" I'm looking for increase monotonically with the ODE time variable, so there's no concern about missing events, but I have thirty thousand different output values I want to receive information about. They clump, so there might be many per solver step, but more to the point, they all use the same event function but they're all looking for different values of it. Efficiently handling this case may be non-trivial, but I suspect it's very common - a simpler version would be simply needing to invert a (monotonic) function defined by an ODE, or detecting where an ODE solution crosses each of a grid of y values.

The existing scipy.optimize.brent does a perfectly good job of dealing with a single event, once you have a zero crossing to establish that there is an event. It can't really be used in fully-compiled mode, it doesn't use derivative information, and it doesn't support either lists of distinct event functions or looking for any one of many output values (in fact it only looks for zeros). But with some work, an event detector could be added to scipy.optimize for general use and then applied here. I do think it might make sense to integrate event detection into the interfaces of scipy.integrate since it is such a common application of ODE solving.

Contributor

aarchiba commented Aug 15, 2016

I'm not sure how to allow users to control finer details of event detection - high-order methods like Bulirsch-Stoer may take quite large steps, possibly large enough to miss events if the event function is non-monotonic

These would be events that are changing faster than the states they are tracking, am I right? I've never seen methods for dealing with things like this. This would seem to be a good argument for having the event detection being built on top of the step-by-step solver rather than into it so you could swap out the standard detector for your own method.

Well, there are two situations I'm concerned about:

  • The hypothetical situation of using event detection on a non-monotonic function, for example, using an ODE solver to get sin(t) and event detection to find where it crosses a level y. Here the concern is that the solver might take big steps - B-S can manage good accuracy for this problem using step sizes of order unity - and you might miss an event because the function crosses the threshold twice. Here the user might need to specify a maximum step size, or better, a brute-force sub-sampling of the solver step for use in event detection. Derivative information might help here too, or the fact that for most solver methods the dense output actually comes from evaluating a polynomial.
  • The practical situation is that I have a very smooth integration problem where I need to know about the system state at each of a prescribed list of "times", where each "time" is computed from the system state. The "times" I'm looking for increase monotonically with the ODE time variable, so there's no concern about missing events, but I have thirty thousand different output values I want to receive information about. They clump, so there might be many per solver step, but more to the point, they all use the same event function but they're all looking for different values of it. Efficiently handling this case may be non-trivial, but I suspect it's very common - a simpler version would be simply needing to invert a (monotonic) function defined by an ODE, or detecting where an ODE solution crosses each of a grid of y values.

The existing scipy.optimize.brent does a perfectly good job of dealing with a single event, once you have a zero crossing to establish that there is an event. It can't really be used in fully-compiled mode, it doesn't use derivative information, and it doesn't support either lists of distinct event functions or looking for any one of many output values (in fact it only looks for zeros). But with some work, an event detector could be added to scipy.optimize for general use and then applied here. I do think it might make sense to integrate event detection into the interfaces of scipy.integrate since it is such a common application of ODE solving.

@aarchiba

This comment has been minimized.

Show comment
Hide comment
@aarchiba

aarchiba Aug 15, 2016

Contributor

I think perhaps we need to distinguish event detection from root finding.

At the single-stepping level, there's no difference: event detection is root finding within the last step, so there's no real need to have a special interface for it.

At higher levels, where the user wants either an anywhere-evaluable solution object or a result at a list of time values, event detection signals the solver not to bother integrating any more steps past the occurrence of the first event. Root detection, as in my second case above, does not mean stopping the integrator, and so there's not really much need to integrate it into the higher-level interfaces; users can use the anywhere-evaluable solution objects as long as they provide a few samples to help bracketing. So I don't think it's necessary for the "when do I stop?" codes to handle large numbers of stopping points.

Contributor

aarchiba commented Aug 15, 2016

I think perhaps we need to distinguish event detection from root finding.

At the single-stepping level, there's no difference: event detection is root finding within the last step, so there's no real need to have a special interface for it.

At higher levels, where the user wants either an anywhere-evaluable solution object or a result at a list of time values, event detection signals the solver not to bother integrating any more steps past the occurrence of the first event. Root detection, as in my second case above, does not mean stopping the integrator, and so there's not really much need to integrate it into the higher-level interfaces; users can use the anywhere-evaluable solution objects as long as they provide a few samples to help bracketing. So I don't think it's necessary for the "when do I stop?" codes to handle large numbers of stopping points.

@drhagen

This comment has been minimized.

Show comment
Hide comment
@drhagen

drhagen Aug 19, 2016

Contributor

I refactored @nmayorov's solvers into step-by-step classes. With my branch, it should be easier to (1) add additional solvers and (2) add additional functions that work with all solvers (e.g. a function that only returns the solutions at select points). This refactoring also removed a lot of duplicated code between the solvers.

This is my first attempt to contribute to scipy. Should I open a pull request and reference this one?

Contributor

drhagen commented Aug 19, 2016

I refactored @nmayorov's solvers into step-by-step classes. With my branch, it should be easier to (1) add additional solvers and (2) add additional functions that work with all solvers (e.g. a function that only returns the solutions at select points). This refactoring also removed a lot of duplicated code between the solvers.

This is my first attempt to contribute to scipy. Should I open a pull request and reference this one?

@charris

This comment has been minimized.

Show comment
Hide comment
@charris

charris Aug 20, 2016

Member

@drhagen If your changes are based on top of this PR, then @nmayorov can merge them into this if that is appropriate after discussion.

Member

charris commented Aug 20, 2016

@drhagen If your changes are based on top of this PR, then @nmayorov can merge them into this if that is appropriate after discussion.

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Aug 20, 2016

Contributor

@drhagen

It looks very interesting. I'm pleased to see that you figured what was written in my code. Some things I noticed:

  1. Generall.
    1. I think that (t, y) is the best choice of argument names in ODEs, they can be referenced as "time" and "state". And t0, tF (or t0, t_final) is better than a, b. So your choices are good.
  2. For the function ivp_solution.
    1. I guess solve_ivp sounds better and more conventional, because it uses a verb.
    2. Passing class as an argument --- a bit unconventional, maybe better stick to string keywords.
    3. Docstring is not up to date.
    4. I think it might be better to have all options listed as individual arguments (not just **options). It just feels more transparent and consistent. I know it contradicts "plug-and-play" concept, but I think usability and readability is more important for scipy. That was partially the reason I wanted to implement new solvers --- to have all methods and options described thoroughly in one place.
  3. For classes.
    1. Sort of a wide question: do we really need classes? I guess the main reason to have a "stepping" ability. Is it essential? I'm not saying we don't need classes, I just want to have a clear idea what classes will achieve for us (because I don't think they are more convenient per se).
    2. I have a feeling that making equivalent (and faster) Cython code using classes will be more challenging (that's one of the reasons to avoid them). If we have functions that do batch processing (like integration from t0 to tf) the way to convert to Cython is conceptually simple: we have a Python wrapper which delegates the computations to Cython once, then it fetches the results and returns them in a user friendly fashion. If we have classes which are assumed to be called repeatedly, then it is not clear how to make it efficient with Cython.
    3. Python code with classes will be somewhat slower too. You also seem to overuse OOP design, like introducing dummy classes SolverDirection, SolverStatus, it all hurt performance just a little bit more. I'm not really expert in this, but it seems to me such things are not necessary in Python, I would just use strings for this.
    4. I remember that there is a special case for computing interpolants when we terminated by an event, see here for example. But I don't see such lines in the new code. Does your version pass unit tests?
    5. It should be "Runge-Kutta" (Google doesn't find anything for "Runga-Kutta").
Contributor

nmayorov commented Aug 20, 2016

@drhagen

It looks very interesting. I'm pleased to see that you figured what was written in my code. Some things I noticed:

  1. Generall.
    1. I think that (t, y) is the best choice of argument names in ODEs, they can be referenced as "time" and "state". And t0, tF (or t0, t_final) is better than a, b. So your choices are good.
  2. For the function ivp_solution.
    1. I guess solve_ivp sounds better and more conventional, because it uses a verb.
    2. Passing class as an argument --- a bit unconventional, maybe better stick to string keywords.
    3. Docstring is not up to date.
    4. I think it might be better to have all options listed as individual arguments (not just **options). It just feels more transparent and consistent. I know it contradicts "plug-and-play" concept, but I think usability and readability is more important for scipy. That was partially the reason I wanted to implement new solvers --- to have all methods and options described thoroughly in one place.
  3. For classes.
    1. Sort of a wide question: do we really need classes? I guess the main reason to have a "stepping" ability. Is it essential? I'm not saying we don't need classes, I just want to have a clear idea what classes will achieve for us (because I don't think they are more convenient per se).
    2. I have a feeling that making equivalent (and faster) Cython code using classes will be more challenging (that's one of the reasons to avoid them). If we have functions that do batch processing (like integration from t0 to tf) the way to convert to Cython is conceptually simple: we have a Python wrapper which delegates the computations to Cython once, then it fetches the results and returns them in a user friendly fashion. If we have classes which are assumed to be called repeatedly, then it is not clear how to make it efficient with Cython.
    3. Python code with classes will be somewhat slower too. You also seem to overuse OOP design, like introducing dummy classes SolverDirection, SolverStatus, it all hurt performance just a little bit more. I'm not really expert in this, but it seems to me such things are not necessary in Python, I would just use strings for this.
    4. I remember that there is a special case for computing interpolants when we terminated by an event, see here for example. But I don't see such lines in the new code. Does your version pass unit tests?
    5. It should be "Runge-Kutta" (Google doesn't find anything for "Runga-Kutta").
@drhagen

This comment has been minimized.

Show comment
Hide comment
@drhagen

drhagen Aug 20, 2016

Contributor
  1. General

    i. I think that (t, y) is the best choice of argument names in ODEs

There three major conventions out there:

  • x, y: Used by most of the examples in Wikipedia
  • t, x: Used by the first example in Wikipedia. This is the one I am most familiar with (control theory background), where x is the internal (ODE) state and y is the observable state of the system
  • t, y: Used by scipy right now, which is why I used that in my code

I'll leave it up to you which one we settle on. I'll make all the code consistent with this choice.

  1. For the function ivp_solution

    i. I guess solve_ivp sounds better and more conventional, because it uses a verb.

I agree. I was in the process of changing that back. I was using a different name so I could run your version and my version side by side to make sure I didn't mess something up.

Whatever is chosen for this name, it should meld nicely with a variant of the solver that returns the solution at discrete times. In other projects where I have this distinction, I use something like solve_ivp_complete and solve_ivp_discrete.

ii. Passing class as an argument --- a bit unconventional, maybe better stick to string keywords.

If we do this, then the user can never write his own ODE solver and use it with solve_ivp.

iii. Docstring is not up to date.

I'll take care of that if we decide that classes are a worthwhile addition.

iv. I think it might be better to have all options listed as individual arguments (not just **options).

How would you deal with solvers that had different options? You wouldn't be able to add forward and backward Euler because it does not have the concept of rtol and atol. You would be able to specify the max order on variable order solvers.

  1. For classes

    i. Sort of a wide question: do we really need classes?

This is the most important question. Whether or not we should continue down this path hangs on the answer. In my opinion, there are two main benefits to step-by-step solvers:

  1. Step-by-step solvers are useful if you want to do anything new with the existing solvers. For example, if @aarchiba wants to do her 10000 root findings, she can do it pretty easily writing a wrapper around the step solvers. Without them, she would have to write each solver that she wanted to use from scratch because the only interface exposed is solve_ivp. For another example, if we want to add a solver that only keeps the solution at specific points, we can do that pretty easily by wrapping the existing steps and only keeping the state requested. Without them, we have to make a copy of rk(...) and radau(...) because they only provide the entire solution.
  2. Step-by-step solvers are useful if you want to add new solvers to the library without having to rewrite the integration wrappers. You can see this in the duplicated code between rk.py and radau.py. The event detection logic is the same regardless of the underlying solver.

ii. I have a feeling that making equivalent (and faster) Cython code using classes will be more challenging (that's one of the reasons to avoid them).

I have no experience with Cython. If the goal is fast solvers and classes prevent fast solvers from being compiled, then my approach will have to be abandoned.

iii. You also seem to overuse OOP design, like introducing dummy classes SolverDirection, SolverStatus

I'll take that as a compliment. :-) I figured that Cython would handle them efficiently considering that they are static objects. I'll bet they are faster than strings (especially in Cython), but they are probably slower than integers. I'll remove them if that's what you decide.

iv. I remember that there is a special case for computing interpolants when we terminated by an event

I couldn't figure out how that code would change the answer in any meaningful way, so I deleted it. My OdeSolution class prevents getting values past the terminating event time, even if the underlying spline has values for the rest of the step. All the unit tests passed just fine.

v. It should be "Runge-Kutta"

You are correct.

Contributor

drhagen commented Aug 20, 2016

  1. General

    i. I think that (t, y) is the best choice of argument names in ODEs

There three major conventions out there:

  • x, y: Used by most of the examples in Wikipedia
  • t, x: Used by the first example in Wikipedia. This is the one I am most familiar with (control theory background), where x is the internal (ODE) state and y is the observable state of the system
  • t, y: Used by scipy right now, which is why I used that in my code

I'll leave it up to you which one we settle on. I'll make all the code consistent with this choice.

  1. For the function ivp_solution

    i. I guess solve_ivp sounds better and more conventional, because it uses a verb.

I agree. I was in the process of changing that back. I was using a different name so I could run your version and my version side by side to make sure I didn't mess something up.

Whatever is chosen for this name, it should meld nicely with a variant of the solver that returns the solution at discrete times. In other projects where I have this distinction, I use something like solve_ivp_complete and solve_ivp_discrete.

ii. Passing class as an argument --- a bit unconventional, maybe better stick to string keywords.

If we do this, then the user can never write his own ODE solver and use it with solve_ivp.

iii. Docstring is not up to date.

I'll take care of that if we decide that classes are a worthwhile addition.

iv. I think it might be better to have all options listed as individual arguments (not just **options).

How would you deal with solvers that had different options? You wouldn't be able to add forward and backward Euler because it does not have the concept of rtol and atol. You would be able to specify the max order on variable order solvers.

  1. For classes

    i. Sort of a wide question: do we really need classes?

This is the most important question. Whether or not we should continue down this path hangs on the answer. In my opinion, there are two main benefits to step-by-step solvers:

  1. Step-by-step solvers are useful if you want to do anything new with the existing solvers. For example, if @aarchiba wants to do her 10000 root findings, she can do it pretty easily writing a wrapper around the step solvers. Without them, she would have to write each solver that she wanted to use from scratch because the only interface exposed is solve_ivp. For another example, if we want to add a solver that only keeps the solution at specific points, we can do that pretty easily by wrapping the existing steps and only keeping the state requested. Without them, we have to make a copy of rk(...) and radau(...) because they only provide the entire solution.
  2. Step-by-step solvers are useful if you want to add new solvers to the library without having to rewrite the integration wrappers. You can see this in the duplicated code between rk.py and radau.py. The event detection logic is the same regardless of the underlying solver.

ii. I have a feeling that making equivalent (and faster) Cython code using classes will be more challenging (that's one of the reasons to avoid them).

I have no experience with Cython. If the goal is fast solvers and classes prevent fast solvers from being compiled, then my approach will have to be abandoned.

iii. You also seem to overuse OOP design, like introducing dummy classes SolverDirection, SolverStatus

I'll take that as a compliment. :-) I figured that Cython would handle them efficiently considering that they are static objects. I'll bet they are faster than strings (especially in Cython), but they are probably slower than integers. I'll remove them if that's what you decide.

iv. I remember that there is a special case for computing interpolants when we terminated by an event

I couldn't figure out how that code would change the answer in any meaningful way, so I deleted it. My OdeSolution class prevents getting values past the terminating event time, even if the underlying spline has values for the rest of the step. All the unit tests passed just fine.

v. It should be "Runge-Kutta"

You are correct.

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Apr 20, 2017

Contributor

@rgommers or @pv could you hit the merge button? I'm a bit hesitant do do that :)

Contributor

nmayorov commented Apr 20, 2017

@rgommers or @pv could you hit the merge button? I'm a bit hesitant do do that :)

@pv pv merged commit afc69db into scipy:master Apr 21, 2017

1 of 2 checks passed

continuous-integration/travis-ci/pr The Travis CI build failed
Details
ci/circleci Your tests passed on CircleCI!
Details
@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Apr 21, 2017

Member

Thanks a lot everyone! Let's get this in now, so that it is in early enough for 1.0

Member

pv commented Apr 21, 2017

Thanks a lot everyone! Let's get this in now, so that it is in early enough for 1.0

@pv pv added this to the 1.0 milestone Apr 21, 2017

@WarrenWeckesser

This comment has been minimized.

Show comment
Hide comment
@WarrenWeckesser

WarrenWeckesser Apr 30, 2017

Member

Argh... solve_ivp now shows up over in #7168.

Please include illustrative examples in the docstrings of any new functions added to the public API of scipy.

Member

WarrenWeckesser commented Apr 30, 2017

Argh... solve_ivp now shows up over in #7168.

Please include illustrative examples in the docstrings of any new functions added to the public API of scipy.

@WarrenWeckesser

This comment has been minimized.

Show comment
Hide comment
@WarrenWeckesser

WarrenWeckesser Apr 30, 2017

Member

Sorry--the first thing I meant to say was "Thanks for the all the hard work in getting this new feature into scipy!".

Member

WarrenWeckesser commented Apr 30, 2017

Sorry--the first thing I meant to say was "Thanks for the all the hard work in getting this new feature into scipy!".

@charris

This comment has been minimized.

Show comment
Hide comment
@charris

charris Apr 30, 2017

Member

Oops ;)

Member

charris commented Apr 30, 2017

Oops ;)

@drhagen drhagen referenced this pull request May 20, 2017

Merged

Add examples to solve_ivp #7422

@Juanlu001

This comment has been minimized.

Show comment
Hide comment
@Juanlu001

Juanlu001 Aug 3, 2017

Contributor

If I understood correctly, the best way to emulate the old ode.set_solout in the new API would be using events. Is that correct?

Contributor

Juanlu001 commented Aug 3, 2017

If I understood correctly, the best way to emulate the old ode.set_solout in the new API would be using events. Is that correct?

@drhagen

This comment has been minimized.

Show comment
Hide comment
@drhagen

drhagen Aug 3, 2017

Contributor

If I understood correctly, the best way to emulate the old ode.set_solout in the new API would be using events. Is that correct?

The event function is called at every step, so that is a fine solution. Alternatively, you could set dense_output to True and do whatever you want with the value at every step after completion.

Alternatively, if what you want is particular enough that it won't fit into either of those solutions, you can work with an object-oriented step solver (e.g. scipt.integrate.RK45) directly, taking steps manually and doing whatever you want with the state at each step. The new step solver classes differ from the old ode one in that step actually takes a single integration step. You can do what you want with the ODE state at each step without having to include a set_selout function to capture that for you. And if I remember correctly, set_selout only worked for a subset of the old solvers, while this will work for all the new solvers.

Contributor

drhagen commented Aug 3, 2017

If I understood correctly, the best way to emulate the old ode.set_solout in the new API would be using events. Is that correct?

The event function is called at every step, so that is a fine solution. Alternatively, you could set dense_output to True and do whatever you want with the value at every step after completion.

Alternatively, if what you want is particular enough that it won't fit into either of those solutions, you can work with an object-oriented step solver (e.g. scipt.integrate.RK45) directly, taking steps manually and doing whatever you want with the state at each step. The new step solver classes differ from the old ode one in that step actually takes a single integration step. You can do what you want with the ODE state at each step without having to include a set_selout function to capture that for you. And if I remember correctly, set_selout only worked for a subset of the old solvers, while this will work for all the new solvers.

@Juanlu001

This comment has been minimized.

Show comment
Hide comment
@Juanlu001

Juanlu001 Aug 16, 2017

Contributor

Another question: what would be your approach to integrate functions receiving extra parameters? Using functools.partial, other?

def func(t, y, extra_param):
    # ...
Contributor

Juanlu001 commented Aug 16, 2017

Another question: what would be your approach to integrate functions receiving extra parameters? Using functools.partial, other?

def func(t, y, extra_param):
    # ...
@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Aug 16, 2017

Contributor

Either that, or lambda t, y: func(t, y, extra_params).

I think extra_params could be added to class constructors, but its utility is borderline.

Does anyone think it's better to include args, kwargs for fun into API?

Contributor

nmayorov commented Aug 16, 2017

Either that, or lambda t, y: func(t, y, extra_params).

I think extra_params could be added to class constructors, but its utility is borderline.

Does anyone think it's better to include args, kwargs for fun into API?

@Juanlu001

This comment has been minimized.

Show comment
Hide comment
@Juanlu001

Juanlu001 Aug 16, 2017

Contributor

Does anyone think it's better to include args, kwargs for fun into API?

When I raised this very same question in #2818, there were both supporting and opposing opinions:

Contributor

Juanlu001 commented Aug 16, 2017

Does anyone think it's better to include args, kwargs for fun into API?

When I raised this very same question in #2818, there were both supporting and opposing opinions:

@nmayorov

This comment has been minimized.

Show comment
Hide comment
@nmayorov

nmayorov Aug 16, 2017

Contributor
Contributor

nmayorov commented Aug 16, 2017

@drhagen

This comment has been minimized.

Show comment
Hide comment
@drhagen

drhagen Aug 16, 2017

Contributor

Another question: what would be your approach to integrate functions receiving extra parameters?

I'm with @nmayorov in that I think that adding parameters to the constructor simply to be passed along to fun is of limited utility. Using lambda should be the way to go. If I understand this thread correctly, Matlab only has the ability to pass along extra arguments because it used to not have anonymous function handles. Now that it has them, they have deprecated the old extra-arguments way.

Contributor

drhagen commented Aug 16, 2017

Another question: what would be your approach to integrate functions receiving extra parameters?

I'm with @nmayorov in that I think that adding parameters to the constructor simply to be passed along to fun is of limited utility. Using lambda should be the way to go. If I understand this thread correctly, Matlab only has the ability to pass along extra arguments because it used to not have anonymous function handles. Now that it has them, they have deprecated the old extra-arguments way.

@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Aug 16, 2017

Member

The extra_args stuff in several scipy routines may be because PEP 227 was implemented only around the same time as scipy.integrate, and lambda was not that usable before it.

Member

pv commented Aug 16, 2017

The extra_args stuff in several scipy routines may be because PEP 227 was implemented only around the same time as scipy.integrate, and lambda was not that usable before it.

@Juanlu001

This comment has been minimized.

Show comment
Hide comment
@Juanlu001

Juanlu001 Aug 16, 2017

Contributor

Wow, that's an extremely interesting historical remark...!

Contributor

Juanlu001 commented Aug 16, 2017

Wow, that's an extremely interesting historical remark...!

@ilayn

This comment has been minimized.

Show comment
Hide comment
@ilayn

ilayn Oct 14, 2017

Member

Hi guys, I've just stumbled upon this blog post. I just wanted to put it on the ODE table 😄 just in case you would like to have a look at a nicely compiled summary (I'm out of my depth here so I don't know whether it is biased or not). Maybe it might help to set out some Roadmap items.

http://www.stochasticlifestyle.com/comparison-differential-equation-solver-suites-matlab-r-julia-python-c-fortran/

Member

ilayn commented Oct 14, 2017

Hi guys, I've just stumbled upon this blog post. I just wanted to put it on the ODE table 😄 just in case you would like to have a look at a nicely compiled summary (I'm out of my depth here so I don't know whether it is biased or not). Maybe it might help to set out some Roadmap items.

http://www.stochasticlifestyle.com/comparison-differential-equation-solver-suites-matlab-r-julia-python-c-fortran/

@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Oct 16, 2017

Member
Member

pv commented Oct 16, 2017

@pv

This comment has been minimized.

Show comment
Hide comment
@pv

pv Oct 16, 2017

Member
Member

pv commented Oct 16, 2017

@dimitribalasoiu

This comment has been minimized.

Show comment
Hide comment
@dimitribalasoiu

dimitribalasoiu Nov 10, 2017

Looking through the doc, I see all thoses fancy complicated solvers, but no symplectic one.
Is there a specific reason why there is no symplectic solver in scipy?
I'm writing a (very basic, fixed-step ...) symplectic Euler solver.
Should I make the effort to use OdeSolver? (i.e. would you want it in scipy?)

dimitribalasoiu commented Nov 10, 2017

Looking through the doc, I see all thoses fancy complicated solvers, but no symplectic one.
Is there a specific reason why there is no symplectic solver in scipy?
I'm writing a (very basic, fixed-step ...) symplectic Euler solver.
Should I make the effort to use OdeSolver? (i.e. would you want it in scipy?)

@Juanlu001

This comment has been minimized.

Show comment
Hide comment
@Juanlu001

Juanlu001 May 6, 2018

Contributor

@dimitribalasoiu I agree having a symplectic solver would be great. And while we are at it, I also miss the old dop853... (#6326 (comment))

Contributor

Juanlu001 commented May 6, 2018

@dimitribalasoiu I agree having a symplectic solver would be great. And while we are at it, I also miss the old dop853... (#6326 (comment))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment