4

I am heading to understand odeint from c++ boost library and I need to know which part does what. In boost/numeric/odeint/integrate/integrate_adaptive.hpp, there is a function called integrate_adaptive. This function has a few overloads. The simplified file with my few manipulations is as following:

integrate_adaptive_minimal.hpp

#define BOOST_NUMERIC_ODEINT_INTEGRATE_INTEGRATE_ADAPTIVE_HPP_INCLUDED

#include <boost/type_traits/is_same.hpp>

#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/integrate/null_observer.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
using namespace std;
namespace boost {
namespace numeric {
namespace odeint {


/*
 * the two overloads are needed in order to solve the forwarding problem
 */
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
        Stepper stepper , System system , State &start_state ,
        Time start_time , Time end_time , Time dt ,
        Observer observer )
{
    cout<<"type one"<<endl;  //added by me ***************************************
    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
    return detail::integrate_adaptive(
            stepper , system , start_state ,
            start_time , end_time , dt ,
            observer , stepper_category() );
    /*
     * Suggestion for a new extendable version:
     *
     * integrator_adaptive< Stepper , System, State , Time , Observer , typename Stepper::stepper_category > integrator;
     * return integrator.run( stepper , system , start_state , start_time , end_time , dt , observer );
     */
}

/**
 * \brief Second version to solve the forwarding problem,
 * can be called with Boost.Range as start_state.
 */
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
        Stepper stepper , System system , const State &start_state ,
        Time start_time , Time end_time , Time dt ,
        Observer observer )
{
    cout<<"type two"<<endl;  //added by me ***************************************
    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
    return detail::integrate_adaptive(
            stepper , system , start_state ,
            start_time , end_time , dt ,
            observer , stepper_category() );
}


} // namespace odeint
} // namespace numeric
} // namespace boost

The difference between these two functions is related to parameter start_state. In the second function it contains a const keyword. I want to know which overload is called. So, I added a cout in each function. Then I call them from this test file:

minimal.cpp

#include "integrate_adaptive_minimal.hpp"
#include <boost/numeric/odeint.hpp>
#include <armadillo>

using namespace arma;
using namespace boost::numeric::odeint;

typedef vec::fixed<2> state_type;

void sys( const state_type &x , state_type &dxdt , const double t)
{
    mat A;
    A<<-1<<0<<endr<<0<<-1;
    mat B;
    B<<1<<endr<<1;

    dxdt=A*x+B*(t>0?1:0);
}

void observer( const state_type &x , const double t )
{
}

typedef runge_kutta_dopri5<state_type> stepper_type;

int main()
{
    state_type x;
    x(0) = 0.0;
    x(1) = 0.0;
    integrate_adaptive(make_controlled(1E-10,1E-10,stepper_type()),
                        sys,x,0.0,11.0,0.1,observer);
    return 0;
}

When I compile and run:

g++ -std=c++11 minimal.cpp  -larmadillo -lboost_thread -lboost_filesystem -lboost_system -Wfatal-errors

./a.out

The result shows that the first function integrate_adaptive is called which does not have const. Now, I want to know how is the proper way of calling the second function? The object x cannot be const!

If somebody even know about the mechanism of boost library, telling me the advantage of the second function is appreciated.

Update:

original code: github, doc

barej
  • 1,330
  • 3
  • 25
  • 56
  • Add `state_type const cx = x;` and then use `cx` instead of `x` in the call to `integrate_adaptive`. You can also initialize `cx` directly by calling a `vec::fixex<2>` constructor ... something like `state_type const cx(0.0, 0.0);` – Praetorian Jan 13 '15 at 22:20
  • "The result shows". Frankly, I don't see any result and I don't think that's because I'm doing it wrong. What result are you referring to? – sehe Jan 13 '15 at 23:43
  • @sehe `cout<<"type one"< – barej Jan 14 '15 at 08:15

3 Answers3

2

To be honest, I'm a bit confused about the premise of the question.

integrate_adaptive ends up calling that odeint::copy (from odeint/util/copy.hpp) from controlled_step_result try_step(System, const StateInOut&x, time_type&, time_type&). The documentation here states (among others)

  • \param system The system function to solve, hence the r.h.s. of the ODE. It must fulfill the Simple System concept.
  • \param x The state of the ODE which should be solved. Overwritten if the step is successful. Can be a boost range.

I simply don't think you want a const state. Because the function couldn't possibly fulfill it's promise if the state was, actually, const.

I'd assume the overload exists to accomodate state objects that represent mutable state even when const. This seems weird but could be conceivable with a double * const& (as opposed to double const* const&), if you will.

UPDATE

The documentation says for the const overloads:

/*
* the two overloads are needed in order to solve the forwarding problem
*/

Also, it says

/**
* \brief Second version to solve the forwarding problem,
* can be called with Boost.Range as start_state.
*/

So, indeed it could be just to allow calling with boost::make_iterator_range(a,b) where the Boost Range object itself would be a temporary (bound to const&) and the iterators still mutable.


Off topic:

In case you did want to pass it as const, here's how you'd do it (of course it doesn't compile for the reasons mentioned above):

state_type const x { 0, 0 };
integrate_adaptive(make_controlled(1E-10, 1E-10, stepper_type()), sys, x, 0.0, 11.0, 0.1/*, observer*/);

Or even

integrate_adaptive(make_controlled(1E-10, 1E-10, stepper_type()), sys, 
                 state_type { 0.0, 0.0 }, 0.0, 11.0, 0.1/*, observer*/);
sehe
  • 374,641
  • 47
  • 450
  • 633
  • Found some more evidence that indeed the overloads are meant for logically mutable state wrappers (where the wrappers are `const&` temporaries) – sehe Jan 14 '15 at 00:17
2

The state itself will be mutated. So passing a real const object which can not be mutated at all can not be passed to integrate adaptive. The const overload is here for allowing the creation of the state within the call of integrate_adaptive. As sehe already mentioned, it is possible to call integrate adaptive like

integrate_adaptive( stepper , system , make_pair( a.begin() , a.begin() + 3 ) , t0 , t1 , dt , obs );

Here, the state is range of three elements packed into a range. Without the second overload one needs to write

auto r = make_pair( a.begin() , a.begin() + 3 );    //  Without auto the type will be really complicated.
integrate_adaptive( stepper , system , r , t0 , t1 , dt , obs );

In pre C++11 times, the type of the r is really complicated, something like pair< vector< double >::iterator , vector< double >::iterator > and introducing the second overload simplifies a lot. Most of the steppers do_step methods have also similar overloads.

headmyshoulder
  • 6,240
  • 2
  • 20
  • 27
0

You need only to cast x as const as you call integrate_adaptive():

integrate_adaptive(make_controlled(1E-10,1E-10,stepper_type()),
                        sys, const_cast<const state_type &>(x),0.0,11.0,0.1,observer);

The const version of the function promises that it will not modify the contents of x during its scope, but main() can continue modifying it after integrate_adaptive() returns.

(The advantage of using the C++-style const_cast is to catch yourself sooner if you ever decide to change the type of x to some other class.)

Quip11
  • 194
  • 1
  • 7
  • Seems `const_cast(x)` causing so many errors: `In instantiation of ‘_OI std::__copy_move_a(_II, _II, _OI) [with bool _IsMove = false; _II = const double*; _OI = const double*]` in addition to an avalanche of other errors – barej Jan 13 '15 at 22:51
  • Yes, this will not work. integrate_adaptive modifies the state. – headmyshoulder Jan 14 '15 at 07:26
  • No worries; I don't have the declaration of detail::integrate_adaptive() handy, but if it can't take a const reference either, at least the compiler will tell you! :-) – Quip11 Jan 14 '15 at 16:24