[boost] [odeint] Iterator semantics

Subject: [boost] [odeint] Iterator semantics
From: Karsten Ahnert (karsten.ahnert_at_[hidden])
Date: 2012-07-11 12:04:02

Hi,

we have posted a review request for odeint (a library for solving
ordinary differential equations) on this list some days ago. An ODE is
dx/dt = f(x(t)) and one is looking for the solution x(t) which is a
function of t. Solving ODEs is usually an iterative task, so one starts
with x(0) and iteratively applies a solver to obtain a discretized (and
approximated) representation of the solution: x(0) -> x(dt) -> x(2*dt)
-> x(3*dt) -> ...

I played around with an iterator which is doing exactly this iterative
procedure. For example, it can be used via

odeint::runge_kutta4< state_type > stepper;             // ode solver
std::array< double , 3 > x = {{ 10.0 , 10.0 , 10.0 }};  // initial state

double res = std::accumulate(
    // parameters : solver , ode , state , time , time_step
    make_const_step_iterator_begin(stepper, lorenz(), x, 0.0, 0.01) ,
    make_const_step_iterator_end( stepper, lorenz(), x, 1.0, 0.01) ,
    0.0 , []( double sum, const state_type &x) {
          return sum + x[0]; } );

The first iterator increments the time until the time of the second
iterator is rearched t=0.0 -> t=0.01 -> t=0.02 -> .. t=1.0 and the first
component of the solution is accumulated.

In principle, this approach works quite well, but there are some
problems at least semantically. Maybe someone here has some comments or
ideas on this:

The defintion of this iterator is:

template< class Stepper , class System , class StepperTag >
class const_step_iterator : public boost::iterator_facade
<
    const_step_iterator< Stepper , System , StepperTag > ,
    typename Stepper::state_type const ,
    boost::single_pass_traversal_tag
>
{
...
};

It is a single pass iterator, so it can only check for equality of two
iterators. This is problematic since equality here means that the time
(t) of the begin iterator is smaller then the time of the end iterator.
To ensure commutativity of the != or == operation I needed to flag the
begin and end iterator explicitly. This is the reason for the two
factory functions make_const_step_iterator_begin and
make_const_step_iterator_end, which are doing this flagging.

The second semantically problem is that the end iterator in principle
does not need to know the stepper as well as the system (lorenz() in the
above example). But all algorithms from the standard library and
Boost.Range assume that the begin and end iterator are of same type.
Therefore you have to put the stepper and the system into the end
iterator too. Of course, there are techniques like type erasure which
might solve this problem, but they introduce some run-time overhead
which we want to avoid here.

Any ideas or comments about this?

P.S. The code can be found in the iterator branch in odeint:
boost/numeric/odeint/iterator/const_step_iterator.hpp and some examples
libs/numeric/odeint/examples/const_step_iterator.cpp

git clone -b iterator git_at_[hidden]:headmyshoulder/odeint-v2.git