[odeint] gsoc 2011: finalization of odeint
Dear boost community, I would like to propose a project for the GSoC 2011, namely to improve the odeint library currently living in the boost sandbox [1] with the (optimistic) goal to put it on a level suitable for a boost review process. I have been working on this library for about one year now, although I couldn't spent much time on it during the last months. The other author of odeint is my colleague Karsten. We have had contact to this mailing list before [2], getting few, but rather positive feedback. The odeint library is a modern, template based, generic, C++ library for integrating (solving) initial value problems of ordinary differential equations. Such an ode is usually given as dx/dt = f(x,t) with x(t=0) = x0. x here is generally vector valued, possibly high dimensional. The easiest method to solve this problem is the Euler-method, where you introduce a small value Dt and then iterate x(t+Dt) = f(x,t)*Dt. Obviously, in real applications one uses more sophisticated methods, a well known class are the Runge-Kutta methods [3], suitable for many cases. Additionally, there exist special methods for system where the time evolution has conserved quantities, called symplectic methods. Furthermore there are methods for stiff systems and many others. For odeint, our first focus lies on Runge-Kutta methods. We have implemented several such methods together with error control routines and lately added dense output functionality. Our very recent code is a generic Runge-Kutta stepper that allows the user to build any Runge-Kutta method at compile time by providing the appropriate constants and is based on fusion::vectors. Karsten wrote a good article on that [4], unfortunately in German. There are also symplectic routines for Hamiltonian systems and we started to add implicit methods as well. We also have unit tests based on the boost test framework and some basic documentation. The library is currently at a level where it can be used in scientific applications and it is indeed in use by us and several other people in an academic environment. The major enhancement of our library over existing ones like gsl or methods from the numerical recipes is that our routines do not depend on the container type used to represent the vector x. So our library works with the usual std::vector<double> and array<double,N> but also with vector< complex<double> > (already a major improvement over gsl and nr, which we get for free, basically) and blitz++ arrays. The major strength of this approach, however, is seen when you want to solve odes on complex networks where you have to implement your own, non-trivial container type representing some kind of network. With odeint you can fully separate the network-implementation from the ode-solving algorithm. Furthermore, we are able to easily employ the power of CUDA devices to solve odes by the help of the thrust library and using thrust::device_vector as our container type. It is also possible to use odeint with arbitrary precision types like the gmp library. Some more, unfortunately slightly out-dated, information about odeint can be found in [5]. What are my plans for a GSoC? 1. There are still some parts of the library which are not satisfactorily designed and need to be rethought. This is where we could really need the help of an experienced mentor to create a clean, well structured library. 2. Add methods with varying order (burlisch stoer). 3. Add dense output and error control to the generic runge-kutta implementation. 4. Add convenience function that allow users an easier access to the libraries capabilities. 5. Create an observer interface that allows the user to analyze properties of the system while integrating the equations. 6. Make everything ready for a boos review. Points 1-5 can surely be realized during the GSoC period, I believe. And whatever time is left I will spend on documentation and review preparation. A few words to myself: I'm currently a PhD student in physics at the University of Potsdam working on problems in the field of nonlinear dynamics. Solving odes is basically every-day-work for me and I have a used many of the routines in odeint for my scientific research. Thanks for reading. I greatly appreciate any feedback! Regards, Mario Mulansky [1] https://svn.boost.org/svn/boost/sandbox/odeint/ [2] http://boost.2283326.n4.nabble.com/library-proposal-odeint-td2670019.html [3] http://en.wikipedia.org/wiki/Runge–Kutta_methods [4] http://www.c-plusplus.de/forum/282091 [5] http://www.codeproject.com/KB/recipes/odeint.aspx
Hi Mario,
What are my plans for a GSoC? 1. There are still some parts of the library which are not satisfactorily designed and need to be rethought. This is where we could really need the help of an experienced mentor to create a clean, well structured library. 2. Add methods with varying order (burlisch stoer). 3. Add dense output and error control to the generic runge-kutta implementation. 4. Add convenience function that allow users an easier access to the libraries capabilities. 5. Create an observer interface that allows the user to analyze properties of the system while integrating the equations. 6. Make everything ready for a boos review.
I think that this sounds like a great project. I'm not sure how complex the various tasks are, so I couldn't say whether or not the plan is feasible, but I tend towards optimism :) I hope that we can find a mentor with some domain experience. If not, I might be able to help. I don't know a great deal about the domain, but I can help with design issues. I'm looking forward to seeing the full proposal. Andrew
Just a little remark, the current (development) version with all the features Mario described is in https://svn.boost.org/svn/boost/sandbox/odeint/branches/karsten Some Runge-Kutta steppers are missing but it is easy to add them. On 03/21/2011 11:13 PM, Mario Mulansky wrote:
Dear boost community,
I would like to propose a project for the GSoC 2011, namely to improve the odeint library currently living in the boost sandbox [1] with the (optimistic) goal to put it on a level suitable for a boost review process.
Just a little remark,
the current (development) version with all the features Mario described is in
https://svn.boost.org/svn/boost/sandbox/odeint/branches/karsten
Some Runge-Kutta steppers are missing but it is easy to add them.
I would love to see this in Boost. Of course, odeint needs to play well with Boost.Units. Here is how I would expect the code to look (stripped down version of numeric/odeint/examples/lorenz_integrator.cpp) : #include <iostream> #include <iterator> #include <boost/tr1/array.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&); 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; 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; } 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
I would love to see this in Boost. Of course, odeint needs to play well with Boost.Units. Here is how I would expect the code to look (stripped down version of numeric/odeint/examples/lorenz_integrator.cpp) :
It works with Boost.Units, have a look at https://svn.boost.org/svn/boost/sandbox/odeint/branches/karsten/libs/numeric... It also works well with fusion containers, where every entry can have a different dimension.
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 struct lorenz { template< class State , class Value > void operator()( const State &x , State &dxdt, Value t ) { ... } }; Karsten
It works with Boost.Units, have a look at
https://svn.boost.org/svn/boost/sandbox/odeint/branches/karsten/libs/numeric...
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
This isn't a general solution, but works to deduce the function return and argument types :
Ok, I see your point. It is an interesting ansatz and worth to consider. I think Boost.function_types is the boost library for that. In your suggestion one has a stepper for one particular function (signature). This might be suitable for the user, but the original ansatz is more direct, because your are passing the types (state and deriv type) which stores intermediate values. It has also the advantage that you can use templatized functors like struct lorenz { template< class State , class Value > void operator()( const State &x , State &dxdt , Value dxdt ) { ... } }; which allows you to use one system function with different state types (for example with std::vector, boost::range, ublas::vector). This is a feature, which I don't want to miss. Especially the possibility to use boost.range it very nice. Nevertheless, one could build something like your suggestion on top of the existing one, maybe with some metafunctions ... Best regards, Karsten
participants (4)
-
Andrew Sutton -
Karsten Ahnert -
Mario Mulansky -
Matthias Schabel