Re: [Boost-users] [odeint] optional data

Subject: Re: [Boost-users] [odeint] optional data
From: Erik Schultheis (erik.schultheis_at_[hidden])
Date: 2014-07-22 08:44:16

Thanks for the fast response.
>
> Hi
>
> On 07/22/2014 01:24 AM, Erik Schultheis wrote:
> > Dear all,
> >
> > i am simulating particle trajectories, and depending on the data i want
> > to generate i may or may not have to track additional data [boost.ublas
> > c_matrix] along those trajectories. Right now i enable/disable that via
> > a template parameter (so, using std::tuple(vec, vec) or tuple(vec, vec,
> > c_matrix) and fusion_algebra), which has a few drawbacks.
> > The algorithm is already templated to support varying dimensions via a
> > template parameter, so the additional data practically doubles the
> > number of template instantiations for system and all possible observers.
> > This actually got so bad that gcc (mingw gcc 4.8) refuses to compile
> > with debugging enabled. Therefore, i would like to make the switch
> > between tracking the matrix and ignoring it a runtime choice instead of
> > a template parameter.
>
> I don't understand fully what you exactly want to do. Can you show us
> some more code? Can you put the tracking into an observer? I mean, do
> you need to solve an ODE to track the additional data?
>
Unfortunatly, the tracking is an ode itself, so the change depends on 
the position and the current matrix value, so
it cannot be done inside of an observer.

As for the current implementation:
  I have a class BasicTracer<S, I>, where S is the template parameter 
that decides the state type (currently one of std::array<vec<2>, 2>, 
std::array<vec<3>, 2>, std::tuple<vec<2>, <vec<2>, matrix<4>> and 
std::tuple<vec<3>, vec<3>, matrix<6>>) and I the integration startegy 
(either Integrator, that only operates on the vectors, or 
MonodromyIntegrator, that also operates on the matrix). The reason for 
this is that I have two different use cases for the program: 1) a 
human-viewable visualisation of the trajectories that does not need 
matrices, and 2) the actual data i want to generate. Since 1) needs far 
more trajectories to be calculated, it is really helpful that the 
computation time in this case is lower.
With the mentioned limitation (gcc refuses to compile with -g, 
"tracer.o: too many sections (50892)" and "ccQF7l7o.s: Fatal error: 
can't close .obj\Debug\src\tracer.o: File too big") the code works as it 
is now.

For now i have only worked with small problem sizes, where the wait time 
for a result is in the order of a few minutes, however, this is going to 
change and therefore i want to parallelize the code. The differentiation 
between he use cases via templates make the code quite complex at times, 
and i hope to reduce that complexity before i add the parallelization. 
Therefore, I want to transform from a templated approach to a runtime 
approach. I am not sure if it is possible to make the code select the 
dimension at runtime without significant performance impact, but at 
least the second parameter should be removable without increasing 
runtime significantly.

I am not really sure which parts of the code to post to significantly 
inprove my description.
The tracer itself is a type like

template class BasicTracer<std::tuple<vec<3>, vec<3>, monodromy_matrix<3>>, MonodromyIntegrator<std::tuple<vec<3>, vec<3>, monodromy_matrix<3>>>>;

or

template class BasicTracer<std::array<vec<2>, 2>, Integrator<std::array<vec<2>, 2>>>;

the integration is performed like

     double starte = mEnergy(p);

     try

     {

         mMasterObserver.startTrajectory(p);

         perform_integration(mIntegrator, p, 2.0, 1.0 / mPotential.getSize(), mMasterObserver, mAbsErrorBound, mRelErrorBound);

     } catch(int& i) {};

     /// \todo this is the only reason for getLastState's existence

     mEnergyError += std::pow( starte - mEnergy( mMasterObserver.getLastState() ), 2 );

where perform integration just sets parameters to the odeint function:

template<typename System, typename State, typename Observer>
std::size_t perform_integration(const System& system, State& start, double t_end, double dt, Observer& obs, double abs_err, double rel_err)
{
        using namespace boost::numeric::odeint;

        // when range access is possible, try range algebra, otherwise, try fusion
        using algebra_type = typename std::conditional<boost::has_range_iterator<State>::value, range_algebra, fusion_algebra>::type;

        typedef runge_kutta_cash_karp54<State, double, State, double, algebra_type> error_stepper_type;
        typedef controlled_runge_kutta<error_stepper_type> controlled_stepper_type;

        typedef default_error_checker<typename controlled_stepper_type::value_type, typename controlled_stepper_type::algebra_type> error_checker_type;

        controlled_stepper_type stepper(error_checker_type(abs_err, rel_err));
        return integrate_const(stepper, system, start, 0.0, t_end, dt, std::ref(obs) );
}

and the operator() of the integrators look like

template<class S>

void Integrator<S>::operator()( const S& state, S& deriv, double/*t*/ ) const

{

        using namespace system_state;

        auto position = get_position(state);

        auto p = position * mPotential.getSize();

        double x = p[0];

        double y = p[1];

        // prevent out of bounds access

        if(y < 0 || y > mPotential.getSize()-1 || x < 0 || x > mPotential.getSize()-1)

        {

                throw(1);

        };

        // calculate acceleration

        auto& acceleration = get_velocity(deriv);

        for(int i = 0; i < DIMENSION; ++i)

        {

                /// the compiler should be able to ompimize this

                std::array<int, DIMENSION> didx;

                didx.fill(0);

                didx[i] = 1;

                acceleration[i] = -linearInterpolate( mPotential.getDerivative(didx), p );

        }

        // change in position = current velocity

        get_position( deriv ) = get_velocity( state );

}

template<class S>
void MonodromyIntegrator<S>::operator()( const S& state, S& deriv, double/*t*/ ) const
{
        using namespace system_state;
        auto position = get_position(state);

        auto p = position * mPotential.getSize();
        double x = p[0];
        double y = p[1];

        // prevent out of bounds access
        if(y < 0 || y > mPotential.getSize()-1 || x < 0 || x > mPotential.getSize()-1)
        {
                throw(1);
        };

        // calculate acceleration
        auto& acceleration = get_velocity(deriv);
        for(int i = 0; i < DIMENSION; ++i)
        {
                std::array<int, DIMENSION> didx;
                didx.fill(0);
                didx[i] = 1;
                acceleration[i] = -linearInterpolate( mPotential.getDerivative(didx), p );
        }

        // change in position = current velocity
        get_position( deriv ) = get_velocity( state );

        // change in monodromy
        boost::numeric::ublas::axpy_prod(getMonodromyCoeff( mPotential, p ), get_data(state), get_data(deriv), true);
}

If there is some other part of the code that would be helpful, please 
let me know.

> > However, the cost of tracking that matrix is quite high (that is why it
> > is optional), but simply disableing the d matrix / dt calculation in the
> > system does not speed up the calculation much (86 s -> 58 s vs 27 s
> > without matrix). A lot of time is spent in adding and multiplying
> > numbers that are not used anywhere, or actually even copying them (e.g
> > std::swap(c_matrix, c_matrix) shows up quite high in gprof (even though
> > i am not sure how reliable the -pg option is in conjunction with -O3)).
> >
> > So what i am asking here is for suggestions how to best let odeint treat
> > the c_matrix as uninitialized memory and never perform any operations on
> > that depending on a runtime parameter.
> > Options i thought of were:
> > * deriving a class from c_matrix and providing my own versions of
> > operator overloads that check a flag:
> > as far as i know, boost.ublas matrices use some template magic to
> > let operators only return expression templates to optimize evaluation if
> > operations are chained, and i would not be sure about how to actually
> > retain that behaviour in case the matrix is tracked.
> > * creating my own algebra class:
> > from what i could tell from a few quick glances into the algebra
> > files, this would probably be conceptually easier than the other option,
> > but would require a lot of repetetive coding.
> >
> > Does anyone have suggestion for a more efficient implementation strategy?
>
> Have you considered to use a fusion view which only points to tuple( vec
> , vec ). This could be used as the state type you use. Nevertheless, you
> need to instantiate the stepper with a real resizable state_type which
> you could be again a tuple( vec , vec ).
>
No, I didn't know about fusion view up until now :), so i am going to 
take a look at it. Why would I need a resizable state_type?
Whether I track the matrix or not is decided before the call into the 
integrate function.
>
> >From the options you showed here, an own algebra looks most promising.
> This algebra you also be implemented in terms of the existing
> fusion_algebra and redirect all calls to fusion_algebra.