Subject: [Boost-commit] svn:boost r71323 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor: . taylor_v2
From: karsten.ahnert_at_[hidden]
Date: 2011-04-16 11:41:39
Author: karsten
Date: 2011-04-16 11:41:37 EDT (Sat, 16 Apr 2011)
New Revision: 71323
URL: http://svn.boost.org/trac/boost/changeset/71323
Log:
taylor stepper v2
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_controller.hpp (props changed)
- copied unchanged from r71248, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor_controller.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_lorenz_evol_adaptive_statistics.cpp (props changed)
- copied unchanged from r71248, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v1/taylor_lorenz_evol_adaptive_statistics.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_lorenz_evol_direct.cpp (contents, props changed)
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes | 11 +
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/Jamfile | 2
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor.hpp | 340 ++++++++++++++++++++-------------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_main.cpp | 2
4 files changed, 187 insertions(+), 168 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes 2011-04-16 11:41:37 EDT (Sat, 16 Apr 2011)
@@ -1,5 +1,16 @@
+Version 1:
+
+* Done
+* Improve Binomial coefficient calculation
+
+
Version 2:
+* Requiremnts:
+ * Solve Lorenz with help of transforms and compare performance with Version 1
+
+Version 3:
+
* Requirements:
* Optimize a * x
* Optimize a + x
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/Jamfile 2011-04-16 11:41:37 EDT (Sat, 16 Apr 2011)
@@ -11,3 +11,5 @@
;
exe taylor_main : taylor_main.cpp ;
+exe taylor_lorenz_evol_direct : taylor_lorenz_evol_direct.cpp ;
+exe taylor_lorenz_evol_adaptive_statistics : taylor_lorenz_evol_adaptive_statistics.cpp ;
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor.hpp 2011-04-16 11:41:37 EDT (Sat, 16 Apr 2011)
@@ -10,6 +10,9 @@
#include <tr1/array>
+#include <iostream>
+using namespace std;
+
// general boost includes
#include <boost/static_assert.hpp>
#include <boost/math/special_functions/binomial.hpp>
@@ -19,6 +22,7 @@
#include <boost/fusion/include/is_sequence.hpp>
#include <boost/fusion/include/size.hpp>
#include <boost/fusion/include/at.hpp>
+#include <boost/fusion/include/transform.hpp>
// mpl includes
#include <boost/mpl/range_c.hpp>
@@ -39,99 +43,207 @@
{
namespace proto = boost::proto;
+ namespace fusion = boost::fusion;
+ namespace proto = boost::proto;
- template<int I> struct placeholder {};
-
- proto::terminal< placeholder< 0 > >::type const arg1 = {{}};
- proto::terminal< placeholder< 1 > >::type const arg2 = {{}};
- proto::terminal< placeholder< 2 > >::type const arg3 = {{}};
-
+ template< typename I > struct placeholder : I {};
+ proto::terminal< placeholder< mpl::size_t< 0 > > >::type const arg1 = {};
+ proto::terminal< placeholder< mpl::size_t< 1 > > >::type const arg2 = {};
+ proto::terminal< placeholder< mpl::size_t< 2 > > >::type const arg3 = {};
- template< class State , size_t Order >
- struct context_derivs : proto::callable_context< context_derivs< State , Order > const >
+ template< typename I >
+ std::ostream& operator<<( std::ostream &s , const placeholder< I > &p )
{
- typedef double result_type;
+ s << "placeholder<" << I::value << ">";
+ return s;
+ }
- const State &m_x;
- std::tr1::array< State , Order > &m_derivs;
- size_t m_which;
- context_derivs( const State &x , std::tr1::array< State , Order > &derivs , size_t which )
- : m_x( x ) , m_derivs( derivs ) , m_which( which ) { }
+ /*
+ * eventuel array< double , n > dt mit potenzen anlegen
+ */
+ template< class State , size_t MaxOrder >
+ struct taylor_context
+ {
+ typedef State state_type;
+ typedef std::tr1::array< state_type , MaxOrder > deriv_type;
+
+ size_t which;
+ double dt;
+ const state_type &x;
+ deriv_type &derivs;
+ taylor_context( size_t _which , double _dt , const state_type &_x , deriv_type &_derivs )
+ : which( _which ) , dt( _dt ) , x( _x ) , derivs( _derivs ) { }
+ };
- template< int I >
- double operator()( proto::tag::terminal , placeholder<I> ) const
+ template< typename Grammar >
+ struct plus_transform : proto::transform< plus_transform< Grammar > >
+ {
+ template< typename Expr , typename State , typename Data >
+ struct impl : proto::transform_impl< Expr , State , Data >
{
- return ( m_which == 0 ) ? m_x[I] : m_derivs[m_which-1][I];
- }
+ typedef double result_type;
- double operator()( proto::tag::terminal , double x ) const
- {
- return ( m_which == 0 ) ? x : 0.0;
- }
+ result_type operator ()(
+ typename impl::expr_param expr ,
+ typename impl::state_param state ,
+ typename impl::data_param data ) const
+ {
+ return Grammar()( proto::left( expr ) , state , data )
+ + Grammar()( proto::right( expr ) , state , data );
+ }
+ };
+ };
- template< typename L , typename R >
- double operator ()( proto::tag::plus , L const &l , R const &r ) const
+ template< typename Grammar >
+ struct minus_transform : proto::transform< minus_transform< Grammar > >
+ {
+ template< typename Expr , typename State , typename Data >
+ struct impl : proto::transform_impl< Expr , State , Data >
{
- return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) +
- proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
- }
+ typedef double result_type;
+
+ result_type operator ()(
+ typename impl::expr_param expr ,
+ typename impl::state_param state ,
+ typename impl::data_param data ) const
+ {
+ return Grammar()( proto::left( expr ) , state , data )
+ - Grammar()( proto::right( expr ) , state , data );
+ }
+ };
+ };
- template< typename L , typename R >
- double operator ()( proto::tag::minus , L const &l , R const &r ) const
+
+ template< typename Grammar >
+ struct multiplies_transform : proto::transform< multiplies_transform< Grammar > >
+ {
+ template< typename Expr , typename State , typename Data >
+ struct impl : proto::transform_impl< Expr , State , Data >
{
- return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) -
- proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
- }
+ typedef double result_type;
+ result_type operator()(
+ typename impl::expr_param expr ,
+ typename impl::state_param state ,
+ typename impl::data_param data ) const
+ {
+ typedef typename impl::data data_type;
+
+ double tmp = 0.0;
+ for( size_t k=0 ; k<=data.which ; ++k )
+ {
+ data_type data1( k , data.dt , data.x , data.derivs );
+ data_type data2( data.which - k , data.dt , data.x , data.derivs );
+
+ tmp += boost::math::binomial_coefficient< double >( data.which , k )
+ * Grammar()( proto::left( expr ) , state , data1 )
+ * Grammar()( proto::right( expr ) , state , data2 );
+ }
- template< typename L , typename R >
- double operator ()( proto::tag::multiplies , L const &l , R const &r ) const
+ return tmp;
+ }
+ };
+ };
+
+
+ struct terminal_double_transform : proto::transform< terminal_double_transform >
+ {
+ template< typename Expr , typename State , typename Data >
+ struct impl : proto::transform_impl< Expr , State , Data >
{
- double tmp = 0.0;
- for( size_t k=0 ; k<=m_which ; ++k )
+ typedef double result_type;
+
+ result_type operator()(
+ typename impl::expr_param expr ,
+ typename impl::state_param state ,
+ typename impl::data_param data ) const
{
- tmp += boost::math::binomial_coefficient< double >( m_which , k )
- * proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , k ) )
- * proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which - k ) );
+ return ( data.which == 0 ) ? proto::value( expr ) : 0.0;
}
- return tmp;
- }
+ };
+ };
- template< typename L , typename R >
- double operator ()( proto::tag::divides , L const &l , R const &r ) const
+ template< typename Index >
+ struct terminal_placeholder_transform : proto::transform< terminal_placeholder_transform< Index > >
+ {
+ template< typename Expr , typename State , typename Data >
+ struct impl : proto::transform_impl< Expr , State , Data >
{
- return 0.0;
- }
+ typedef double result_type;
+
+ result_type operator()(
+ typename impl::expr_param expr ,
+ typename impl::state_param state ,
+ typename impl::data_param data ) const
+ {
+ typedef typename impl::expr expr_type;
+ typedef typename expr_type::proto_args args_type;
+ typedef typename args_type::child0 index_type;
+ const size_t index = index_type::value;
+ return ( data.which == 0 ) ? data.x[ index ] : data.derivs[ data.which - 1 ][ index ];
+ }
+ };
};
+
+ struct taylor_double_terminal : proto::when< proto::terminal< double > , terminal_double_transform > { };
+
+ template< typename Which >
+ struct taylor_placeholder_terminal : proto::when< proto::terminal< placeholder< Which > > , terminal_placeholder_transform< Which > > { };
+
+ template< typename Grammar >
+ struct taylor_plus : proto::when< proto::plus< Grammar , Grammar > , plus_transform< Grammar > > { };
+
+ template< typename Grammar >
+ struct taylor_minus : proto::when< proto::minus< Grammar , Grammar > , minus_transform< Grammar > > { };
+
+ template< typename Grammar >
+ struct taylor_multiplies : proto::when< proto::multiplies< Grammar , Grammar > , multiplies_transform< Grammar > > { };
+
+
+
+ struct taylor_transform :
+ proto::or_
+ <
+ taylor_double_terminal ,
+ taylor_placeholder_terminal< proto::_ > ,
+ taylor_plus< taylor_transform > ,
+ taylor_minus< taylor_transform > ,
+ taylor_multiplies< taylor_transform >
+ >
+ { };
+
+
+
template< class System , class State , size_t Order >
struct eval_derivs
{
- typedef std::tr1::array< State , Order > derivs_type;
+ typedef State state_type;
+ typedef taylor_context< state_type , Order > taylor_context_type;
+ typedef typename taylor_context_type::deriv_type deriv_type;
System m_sys;
- const State &m_x;
- derivs_type &m_derivs;
- size_t m_which;
- double m_dt;
+ taylor_context_type m_data;
- eval_derivs( System sys , const State &x , derivs_type &derivs , size_t which , double dt )
- : m_sys( sys ) , m_x( x ) , m_derivs( derivs ) , m_which( which ) , m_dt( dt ) { }
+ eval_derivs( System sys , const State &x , deriv_type &derivs , size_t which , double dt )
+ : m_sys( sys ) , m_data( which , dt , x , derivs ) { }
template< class Index >
void operator()( Index )
{
- m_derivs[m_which][ Index::value ] = m_dt / double( m_which + 1 ) * boost::proto::eval(
- boost::fusion::at< Index >( m_sys ) ,
- taylor_adf::context_derivs< State , Order >( m_x , m_derivs , m_which ) );
+ typedef typename fusion::result_of::at< System , Index >::type expr_type;
+ const expr_type &expr = boost::fusion::at< Index >( m_sys );
+
+ double deriv = taylor_transform()( expr , 0.0 , m_data );
+ m_data.derivs[ m_data.which ][ Index::value ] = m_data.dt / double( m_data.which + 1 ) * deriv;
}
};
@@ -142,7 +254,7 @@
}
-
+ struct optimize_tree : proto::or_< proto::_ > { };
}
template< size_t N , size_t Order >
@@ -193,9 +305,13 @@
BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+ typedef typename boost::fusion::result_of::transform< System , taylor_adf::optimize_tree >::type optimized_system_type;
+
+ optimized_system_type optimized_system = boost::fusion::transform( sys , taylor_adf::optimize_tree() );
+
for( size_t i=0 ; i<Order ; ++i )
{
- boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( sys , in , m_derivs , i , dt ) );
+ boost::mpl::for_each< boost::mpl::range_c< size_t , 0 , dim > >( make_eval_derivs( optimized_system , in , m_derivs , i , dt ) );
}
for( size_t i=0 ; i<N ; ++i )
@@ -252,116 +368,6 @@
-namespace taylor_adf2
-{
-
- namespace proto = boost::proto;
-
- template<int I> struct placeholder {};
-
- proto::terminal< placeholder< 0 > >::type const arg1 = {{}};
- proto::terminal< placeholder< 1 > >::type const arg2 = {{}};
- proto::terminal< placeholder< 2 > >::type const arg3 = {{}};
-
-
-
- template< class State , size_t Order >
- struct context_derivs : proto::callable_context< context_derivs< State , Order > const >
- {
- typedef double result_type;
-
- const State &m_x;
- std::tr1::array< State , Order > &m_derivs;
- size_t m_which;
-
- context_derivs( const State &x , std::tr1::array< State , Order > &derivs , size_t which )
- : m_x( x ) , m_derivs( derivs ) , m_which( which ) { }
-
-
-
- template< int I >
- double operator()( proto::tag::terminal , placeholder<I> ) const
- {
- return ( m_which == 0 ) ? m_x[I] : m_derivs[m_which-1][I];
- }
-
- double operator()( proto::tag::terminal , double x ) const
- {
- return ( m_which == 0 ) ? x : 0.0;
- }
-
-
- template< typename L , typename R >
- double operator ()( proto::tag::plus , L const &l , R const &r ) const
- {
- return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) +
- proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
- }
-
-
- template< typename L , typename R >
- double operator ()( proto::tag::minus , L const &l , R const &r ) const
- {
- return proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , m_which ) ) -
- proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which ) );
- }
-
-
- template< typename L , typename R >
- double operator ()( proto::tag::multiplies , L const &l , R const &r ) const
- {
- double tmp = 0.0;
- for( size_t k=0 ; k<=m_which ; ++k )
- {
- tmp += boost::math::binomial_coefficient< double >( m_which , k )
- * proto::eval( l , context_derivs< State , Order >( m_x , m_derivs , k ) )
- * proto::eval( r , context_derivs< State , Order >( m_x , m_derivs , m_which - k ) );
- }
- return tmp;
- }
-
-
- template< typename L , typename R >
- double operator ()( proto::tag::divides , L const &l , R const &r ) const
- {
- return 0.0;
- }
-
- };
-
- template< class System , class State , size_t Order >
- struct eval_derivs
- {
- typedef std::tr1::array< State , Order > derivs_type;
-
- System m_sys;
- const State &m_x;
- derivs_type &m_derivs;
- size_t m_which;
-
- eval_derivs( System sys , const State &x , derivs_type &derivs , size_t which )
- : m_sys( sys ) , m_x( x ) , m_derivs( derivs ) , m_which( which ) { }
-
- template< class Index >
- void operator()( Index )
- {
- m_derivs[m_which][ Index::value ] = boost::proto::eval(
- boost::fusion::at< Index >( m_sys ) ,
- taylor_adf::context_derivs< State , Order >( m_x , m_derivs , m_which ) );
- }
- };
-
- template< class System , class State , size_t Order >
- eval_derivs< System , State , Order > make_eval_derivs( System sys , const State &x , std::tr1::array< State , Order > &derivs , size_t i )
- {
- return eval_derivs< System , State , Order >( sys , x , derivs , i );
- }
-
-
-
-}
-
-
} // namespace odeint
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_lorenz_evol_direct.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_lorenz_evol_direct.cpp 2011-04-16 11:41:37 EDT (Sat, 16 Apr 2011)
@@ -0,0 +1,80 @@
+/*
+ * main.cpp
+ *
+ * Created on: Apr 2, 2011
+ * Author: karsten
+ */
+
+#include <iostream>
+
+#include "taylor.hpp"
+
+#include <boost/numeric/odeint/stepper/explicit_rk4.hpp>
+
+#include <boost/fusion/include/make_vector.hpp>
+
+template< typename T , size_t N >
+std::ostream& operator<<( std::ostream& out , const std::tr1::array< T , N > &x )
+{
+ if( !x.empty() ) out << x[0];
+ for( size_t i=1 ; i<x.size() ; ++i ) out << "\t" << x[i];
+ return out;
+}
+
+typedef boost::numeric::odeint::taylor< 3 , 11 > taylor_type;
+typedef taylor_type::state_type state_type;
+typedef boost::numeric::odeint::explicit_rk4< state_type > rk4_type;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+}
+
+
+
+
+
+namespace fusion = boost::fusion;
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+using boost::numeric::odeint::taylor_adf::arg1;
+using boost::numeric::odeint::taylor_adf::arg2;
+using boost::numeric::odeint::taylor_adf::arg3;
+
+int main( int argc , char **argv )
+{
+ taylor_type tay;
+ rk4_type rk4;
+
+
+ state_type x1 = {{ 10.0 , 10.0 , 10.0 }} , x2 = x1 , xerr = x1;
+
+ const double dt = 0.01;
+ double t = 0.0;
+ for( size_t i=0 ; i<10000 ; ++i , t += dt )
+ {
+ tay.do_step(
+ fusion::make_vector
+ (
+ sigma * ( arg2 - arg1 ) ,
+ R * arg1 - arg2 - arg1 * arg3 ,
+ arg1 * arg2 - b * arg3
+ ) ,
+ x1 , t , dt , xerr );
+
+ rk4.do_step( lorenz , x2 , t , dt );
+
+ cout << t << "\t" << x1 << "\t" << x2 << "\n";
+ }
+
+
+ return 0;
+}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_main.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_main.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v2/taylor_main.cpp 2011-04-16 11:41:37 EDT (Sat, 16 Apr 2011)
@@ -71,8 +71,8 @@
+ cout << endl;
cout << in << endl;
- cout << dxdt << endl;
cout << xerr << endl;
cout << out << endl << endl;
const derivs_type &derivs = t.get_last_derivs();