Re: [boost] [odeint] Iterator semantics

Subject: Re: [boost] [odeint] Iterator semantics
From: Jeffrey Lee Hellrung, Jr. (jeffrey.hellrung_at_[hidden])
Date: 2012-07-11 18:39:30

On Wed, Jul 11, 2012 at 2:57 PM, Karsten Ahnert
<karsten.ahnert_at_[hidden]>wrote:

> On 07/11/2012 07:44 PM, Jeffrey Lee Hellrung, Jr. wrote:
> > On Wed, Jul 11, 2012 at 9:04 AM, Karsten Ahnert
> > <karsten.ahnert_at_[hidden]>wrote:
> >
> >> 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.
> >>
> >
> > Can you elaborate? I don't understand why this is what equality means.
>
> One needs to implement the != operator (or equal if you use
> Boost.Iterator). A naive implementation of operator!= (or the equal() if
> you use boost.iterator) could check that the current time of the ODE is
> smaller then the time of other iterator:
>
> bool equal( iterator other )
> {
>     return this->t < other->t;
> }
>

:: confused :: I would think iterator1 op iterator2 is equivalent to time1
op time2.

Of course, this is not a check for equality but a check for less_then.
> It will also destroy commutativity of the !=operation and will not work
> in general. So you need somehow introduce which iterator is the first
> and which one is the last. I have done this with a flag which says if
> the iterator is the first and which the last such that the operation is now
>
> bool equal( iterator other )
> {
>      return ( this->first ) ?
>          ( this->t < other->t ) :
>          ( this->t > other->t );
> }
>
> >
> >
> >> 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.
> >>
> >
> > I might understand your reasoning for doing this better if I understood
> > your rationale for equality semantics...sorry, maybe I'm being dense :/
> >
> >
> >> 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.
> >
> >
> > Ugh, yes, this is an unfortunate consequence of the design of
> > iterator-based, STL-like algorithms. Perhaps it would be better to simply
> > define range-equivalents of your iterators, and encourage use of the
> > Boost.Range algorithms (which, AFAIK, under the hood, use the STL
> > iterator-based algorithms)? That would eliminate a lot of the repetition
> > you're seeing in defining begin and end iterators separately.
>
> Yeah, we already have a factory function returning a range. You could
> also write the upper example as follows:
>
> double res = boost::range::accumulate(
>      make_const_step_range(stepper, lorenz(), x, 0.0, 1.0 , 0.01 ) ,
>      0.0 , []( double sum, const state_type &x) {
>            return sum + x[0]; } );
>
> But I fear, this does not solve this problem. Boost.Range also expects
> that begin and end are of the same type.
>

Sounds like you want your iterator to be a wrapper around optional< state
>. The state also encapsulates the end condition, and an end iterator is an
uninitialized optional. When incrementing an iterator, if you've reached
the end condition, uninitialize (i.e., reset) the underlying optional. ==
and != on the iterators are equivalent to those on the optionals.
Equivalently, if your state is completely external (and you don't need to
manage it within the iterator), you can use a state* and use its nullness
to distguish end iterators from not-end iterators.

> 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.
> >>
> >
> > I don't see type erasure solving your problem, but regarding the overhead
> > of type erasure only, it could be minimal relative to the cost of the
> > actual time stepper. But that observation probably isn't relevant to your
> > problem.
>
> With type erasure it is in principle possible to remove the dependency
> of iterator from the stepper and the system and put this in some member
> class, maybe something like
>
> template< class State , class StepperTag >
> class const_step_iterator : public boost::iterator_facade
> <
>     const_step_iterator< State , StepperTag > ,
>     State const ,
>     boost::single_pass_traversal_tag
> >
> {
> public:
>
>     template< class Stepper , class System >
>     const_step_iterator( Stepper stepper , System system ,
>         State &x , Time t_start , Time dt )
>     {
>         this->data_ = new DataStepper< Stepper , System >(
>             stepper , system , State , t , dt );
>     }
>
>     template< class Stepper , class System >
>     const_step_iterator( Time t_end )
>     {
>         this->data_ = new DataTime( t );
>     }
>
>     struct IData { ... };
>
>     template< class Stepper , class System >
>     struct DataStepper : public IData { ... };
>
>     struct DataTime : public IData { ... };
>
>      ...
>
> private:
>
>     std::unique_ptr< IData > data_;
> };
>

I don't see what this buys you relative to the original problem.

- Jeff