Re: [boost] [odeint] gsoc 2011: finalization of odeint

Subject: Re: [boost] [odeint] gsoc 2011: finalization of odeint
From: Matthias Schabel (boost_at_[hidden])
Date: 2011-03-22 18:35:39

> It works with Boost.Units, have a look at
> 
> https://svn.boost.org/svn/boost/sandbox/odeint/branches/karsten/libs/numeric/odeint/test/stepper_with_units.cpp
> 
> It also works well with fusion containers, where every entry can have a
> different dimension.

Great...

>> In particular, the steppers should be able to take a function pointer or other function object and derive the container types directly from it rather than requiring that these be provided as template parameters. Then the user need only provide the function to integrate (optionally providing a Jacobian for stiff solvers).
>> 
>> Matthias
>> 
> 
> This seems to be difficult. The container type is used to store
> intermediate values within the steppers, such that objects like
> 
>   container_type m_dxdt;
> 
> exist. I don't see a way to avoid this (without loss of performance).
> But, it might be possible to use a default type for the internal
> containers (like vector< double >) and then use a templatized functors, like

This isn't a general solution, but works to deduce the function return and argument types :

#include <iostream>
#include <iterator>

#include <boost/tr1/array.hpp>

#include <boost/units/detail/utility.hpp>
#include <boost/units/systems/si.hpp>

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/container_traits_tr1_array.hpp>

#define tab "\t"

using namespace std;
using namespace boost;
using namespace boost::units;
using namespace boost::numeric::odeint;

typedef quantity<si::time>					time_type;
typedef array<quantity<si::length>, 3>		position_type;
typedef array<quantity<si::velocity>, 3>	velocity_type;

const double eps_abs = 1E-6;
const double eps_rel = 1E-7;

const size_t time_points = 101;

typedef velocity_type (function_type)(const position_type&,const time_type&);

template<class F> class getFunctionArgumentTypes;

template<class R,class A1,class A2> 
struct getFunctionArgumentTypes<R (*)(const A1&,const A2&)>
{
        typedef R (function_type)(const A1&,const A2&);
        
        typedef R	return_type;
        typedef A1	arg1_type;
        typedef A2	arg2_type;
        
        getFunctionArgumentTypes(function_type f)
        {
                return_type	ret;
                arg1_type	arg1;
                arg2_type	arg2;

                std::cout << boost::units::detail::demangle(typeid(f).name()) << std::endl
                                  << boost::units::detail::demangle(typeid(ret).name()) << std::endl
                                  << boost::units::detail::demangle(typeid(arg1).name()) << std::endl
                                  << boost::units::detail::demangle(typeid(arg2).name()) << std::endl;			  
        }
};

template<class F> 
void identifyFunction(F f)
{
        getFunctionArgumentTypes<F>	gfa(f);
}

velocity_type lorenz(const position_type &x,const time_type& t )
{
        const quantity<si::frequency>	sigma(10.0*si::hertz);
        const quantity<si::length>		rho(28.0*si::meter),
                                                                        beta((8.0/3.0)*si::meter);
        
    const array<quantity<si::velocity>,3>	dxdt = { (x[1] - x[0])*sigma,
                                                                                                         (x[0]*(rho - x[2])/(1.0*si::meter) - x[1])/si::second,
                                                                                                         (x[0]*x[1] - beta * x[2])/si::meter/si::second };
                                                                                                         
        return dxdt;
}

int main( int argc , char **argv )
{
    position_type x1;

    x1[0] = 1.0*si::meters;
    x1[1] = 0.0*si::meters;
    x1[2] = 20.0*si::meters;

    vector<position_type>	x_vec;
    vector<time_type>		t_vec;
        
        identifyFunction(lorenz);
        
//    stepper_half_step< stepper_euler< function_type > > euler(lorenz);
//
//    size_t steps = integrate( euler, 
//							  x1, 
//                              0.0*si::seconds, 
//							  10.0*si::seconds, 
//							  1E-4*si::seconds, 
//                              back_inserter(t_vec),
//                              back_inserter(x_vec));
//
//    clog << "Euler Half Step: #steps " << steps << endl;
//
//    cout.precision(16);
//    cout.setf(ios::fixed,ios::floatfield);
//
//    cout << t_vec.size() << endl;
//    
//    cout << (x_vec.back())[0] << tab << (x_vec.back())[1] << tab << (x_vec.back())[2] << tab;

    return 0;
}

I'm fairly certain that a more generic version of code for deducing function argument/return types exists in Boost, just not sure where...

Matthias