Re: [boost] [odeint] Iterator semantics

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

On 07/12/2012 12:39 AM, Jeffrey Lee Hellrung, Jr. wrote:
> 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.

No, this is not the case iterator != iterator2 is equivalaent to time1 <
time2 . The first operation is commutative the second one not.

> 
> 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.

Yes, maybe optional might be a good idea. I will have a look at this.

[..]