Hi Rhys, thank your for your review and your comments. On 09/28/2012 03:59 PM, Rhys Ulerich wrote:
Here is my review of the proposed ODEInt library...
I'm abstaining on the yes/no library acceptance vote.
- What is your evaluation of the design?
As a user, I can find my way around. Several particular nits appear in my miscellaneous comments section. I'm not convinced that the interfaces, which differ from many other ODE integration packages, are a marked improvement from those other packages (e.g. GSL, PETSc, Trilinos, SUNDIALS, etc). I think that adding more than classical use cases could be difficult, in particular IMEX schemes. In particular, my personal research code (a turbulence DNS platform requiring semi-implicit time advance for a hyperbolic PDE system) could not be retrofit into this framework.
The libraries you mention have a lot of boilerplate which is actually not needed, look for example at PETSc van der Pol example which is 200 lines of code. In odeint you can have the same with 30 lines. The same holds for gsl. Furthermore, they follow more or less a C style coding approach. odeint uses and supports modern C++ techniques, like RAII, iterators, lambdas, ... Why should IMEX schemes not fit into the design? I think they would only require a new system concept. For example the symplectic steppers are semi-explicit methods which have their own concept and fit very well into odeint.
- What is your evaluation of the implementation?
I've not read through the implementation.
- What is your evaluation of the documentation?
Fairly good on the whole. I wish the math more clearly appeared using display-style equations rather than inline text. One table appears repeatedly and many links are missing. The applications section is quite nice.
Your are right. Display-style equations are much better then inline text. It is on our ToDo list. The table is duplicated, but I think it does not hurt.
- What is your evaluation of the potential usefulness of the library?
I think the collection of features are useful, but I'd not adopt this library over the GSL, PETSc's time steppers, or Trilino's Rythmos when working in C or C++. GSL's APIs are quite easy and have no learning curve. PETSc's steppers are very powerful and support many new, exotic schemes. Rythmos build atop clean C++ vector space abstractions from the rest of Trilinos. odeint-v2 doesn't seem to beat any of these three alternatives at their game nor does it feel like enough of a winning compromise to make it compelling to me.
- Did you try to use the library? With what compiler? Did you have any problems?
No.
- How much effort did you put into your evaluation? A glance? A quick reading? In-depth study?
45 minutes of reading.
- Are you knowledgeable about the problem domain?
Somewhat.
In conclusion, I'm thrilled to see someone put forward an ODE integration library for Boost but I feel the design doesn't address current application needs and best practices in many portions of the numerical PDE communities. That said, I applaud the authors' effort and I'm sure odeint-v2 will be useful for other folks.
- Rhys
P.S.: Miscellaneous, jumbled comments:
Within the documentation, Table 1.1's stepper algorithm classes should link to the class. Ditto for lists of models within, e.g., http://headmyshoulder.github.com/odeint-v2/doc/boost_numeric_odeint/concepts... Table 1.2 appears to duplicate Table 1.1 and should simply be linked. Ditto for Table 1.7.
Yes, we will do this.
Stray < and > appear within the generated documentation, e.g. a http://headmyshoulder.github.com/odeint-v2/doc/boost_numeric_odeint/concepts...
sys.first(...) and sys.second(...) don't describe what's happening within the ImplicitSystem interface. sys.calculate_dxdt(...) and sys.calculate_dfdx(...) describe the operations involved. That they're called in some sequence shouldn't determine the names. Moreover, sometimes big gains can be made by simultaneously computing the Jacobian and the function evaluation in a single call. ImplicitSystem should be a System, but it's not. System should use a calculate_dxdt(...) method rather than operator(...) so that one has a hierarchy of system concepts. Similar naming concerns apply to the symplectic system concepts.
sys.first and sys.second come from the underlying std::pair. Using std::pair for the system function is highly generic, it allows you to write stepper.do_step( make_pair( func , jac ) , ... ) where func and jac can be simple C function or functors with the appropriate signature. The same holds for the System concept. It allows you to use simple C-functions as systems or functore with operator() which might be stateful. You are also right, that a function calculating the rhs and the Jacobian is more performant. This is missing and will be fixed in near future.
In- versus out-of-place do_step(...) method signatures should differ by more than the method signatures. In my own codes I've used terminology like 'apply' and 'accumulate' and modeled the signatures off Level 1 and Level 2 BLAS calls.
Distinguishing out-of-place do_step() methods form in-place seems like this seems like a good idea for me. The out out place methods are not in the concept and are not used by integrate() or the iterators such that a change of their name will hopefully not break existing code.
WRT ImplicitSystem, either providing or recommending some mechanism for finite differenced or complex-stepped Jacobians would be useful.
What do you mean by complex-stepper Jacobians?
Naming inconsistencies in the dense output schemes: bulirsch_stoer_dense_out vs rosenbrock4_dense_output
Right, this will be fixed.
Unclear if the design permits low storage Runge Kutta schemes, singly diagonally implicit.schemes, and other fun implicit-explicit things. E.g., there are IMEX schemes with different formal orders of accuracy for the implicit vs explicit portions of the scheme.
I think the orders are not a problem, as well as low storage Runge-Kutta schemes.
Unclear if the design permits CFL-related maximum stable step sizes within f(u). The System interface would need to be adapted to have each method return maximum stable eigenvalue magnitudes for first and second order differential operators.`
I think in this case you need a new kind of stepper and a new kind of system concept. But there is no general problem.
Rather than limiting implicit systems to folks using uBLAS matrices, it may be better to define the abstract operations a class would need to support to plug into the ImplicitSystem infrastructure. Especially because getting a uBLAS "wrapper" around an existing chunk of memory is quite a PITA and relies upon non-standard parts of the uBLAS (See, e.g., http://agentzlerich.blogspot.com/2009/10/boostublasshallowarrayadaptor.html).
Yes, this is also a long wanted feature from us and it will result in a complete redesign of the implicit methods. We will work on this this kind of "ImplicitSystem" infrastructure, but this is a long way and will take some time. Furthermore, it does not affect the existing explicit and symplectic methods. We have also stopped implementing more implicit methods due to the lack of the abstraction.
Why in http://headmyshoulder.github.com/odeint-v2/doc/boost_numeric_odeint/odeint_i... are static_casts preferred to using constructors a la Value(2) / Value(3)?
I think static_cast's are more strict then C-style casts, which is exactly what we want: http://stackoverflow.com/questions/28002/regular-cast-vs-static-cast-vs-dyna... Thank you again for your review.