Re: [Boost-users] Review of ODEInt

Subject: Re: [Boost-users] Review of ODEInt
From: Karsten Ahnert (karsten.ahnert_at_[hidden])
Date: 2012-09-28 11:44:24

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/dense_output_stepper.html
>  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/dense_output_stepper.html

> 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_in_detail/steppers.html
> 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-dynamic-cast

Thank you again for your review.