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

Subject: Re: [Boost-users] [odeint] optional data
From: Karsten Ahnert (karsten.ahnert_at_[hidden])
Date: 2014-07-22 17:29:12

On 07/22/2014 02:44 PM, Erik Schultheis wrote:
> 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.

Ok, very interesting problem. I had a project where we had a given fluid
flow and a system of ODEs describing interacting particles moving in the
fluid. This sounds similar to your project.

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

Is BasicTracer the "main" class and the function you showed here are
methods of BasicTracer? In this case you could remove the dependency of
the integrator into a boost::function or std::function object, like

std::function< void( S const& , S const& , double > mIntegrator;

This remove the integrator template parameter.

If BasicTracer is not the "main" class, can you explain a bit more how
BasicTracer works and what it should do?