[Boost-commit] svn:boost r71331 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor: . taylor_v3

Subject: [Boost-commit] svn:boost r71331 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor: . taylor_v3
From: karsten.ahnert_at_[hidden]
Date: 2011-04-16 17:03:54

Author: karsten
Date: 2011-04-16 17:03:53 EDT (Sat, 16 Apr 2011)
New Revision: 71331
URL: http://svn.boost.org/trac/boost/changeset/71331

Log:
taylor stepper v3, scalar multiplication optimization introduced
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/Jamfile   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor.hpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_controller.hpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_lorenz_evol_adaptive_statistics.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_lorenz_evol_direct.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_main.cpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/Notes |     7 ++++++-                                 
   1 files changed, 6 insertions(+), 1 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 17:03:53 EDT (Sat, 16 Apr 2011)
@@ -1,11 +1,12 @@
 Version 1:
 
 * Done
-* Improve Binomial coefficient calculation
+* Use contexts
 
 
 Version 2:
 
+* Done
 * Requiremnts:
   * Solve Lorenz with help of transforms and compare performance with Version 1
 
@@ -14,6 +15,10 @@
 * Requirements:
   * Optimize a * x 
   * Optimize a + x
+
+
+Version 4:
+  
   * Optimize x * y
   * Optimize g(x)
 

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/Jamfile	2011-04-16 17:03:53 EDT (Sat, 16 Apr 2011)
@@ -0,0 +1,15 @@
+# (C) Copyright 2010 : Karsten Ahnert, Mario Mulansky
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+
+
+project
+    : requirements
+      <define>BOOST_ALL_NO_LIB=1
+      <include>../../../../../..
+    ;
+
+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 ; 

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor.hpp	2011-04-16 17:03:53 EDT (Sat, 16 Apr 2011)
@@ -0,0 +1,482 @@
+/*
+ * taylor.hpp
+ *
+ *  Created on: Apr 2, 2011
+ *      Author: karsten
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_
+#define BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_
+
+#include <tr1/array>
+
+#include <iostream>
+using namespace std;
+
+// general boost includes
+#include <boost/static_assert.hpp>
+#include <boost/math/special_functions/binomial.hpp>
+#include <boost/math/special_functions/factorials.hpp>
+
+// fusion includes
+#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>
+#include <boost/mpl/for_each.hpp>
+
+// proto includes
+#include <boost/proto/proto.hpp>
+
+
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+namespace taylor_adf
+{
+
+	namespace proto = boost::proto;
+	namespace fusion = boost::fusion;
+	namespace proto = boost::proto;
+
+	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< typename I >
+	std::ostream& operator<<( std::ostream &s , const placeholder< I > &p )
+	{
+		s << "placeholder<" << I::value << ">";
+		return s;
+	}
+
+
+	/*
+	 * 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< 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 >
+		{
+			typedef double result_type;
+
+			result_type operator ()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Plus transform" << endl;
+				return Grammar()( proto::left( expr ) , state , data )
+                     + Grammar()( proto::right( expr ) , state , data );
+			}
+		};
+	};
+
+
+	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 >
+		{
+			typedef double result_type;
+
+			result_type operator ()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Minus transform" << endl;
+				return Grammar()( proto::left( expr ) , state , data )
+                     - Grammar()( proto::right( expr ) , state , data );
+			}
+		};
+	};
+
+
+
+	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 >
+		{
+			typedef double result_type;
+
+			result_type operator()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Multiplies transform" << endl;
+				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 );
+				}
+
+				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 >
+		{
+			typedef double result_type;
+
+			result_type operator()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Terminal double transform" << endl;
+				return ( data.which == 0 ) ? proto::value( expr ) : 0.0;
+			}
+		};
+	};
+
+
+	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 >
+		{
+			typedef double result_type;
+
+			result_type operator()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Terminal placeholder transform" << endl;
+				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 ];
+			}
+		};
+	};
+
+	template< typename Grammar >
+	struct right_shift_transform : proto::transform< right_shift_transform< Grammar > >
+	{
+		template< typename Expr , typename State , typename Data >
+		struct impl : proto::transform_impl< Expr , State , Data >
+		{
+			typedef double result_type;
+
+			result_type operator ()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Right shift transform" << endl;
+				return Grammar()( proto::left( expr ) , state , data )
+						+ ( data.which == 0 ) ? proto::value( proto::right( expr ) ) : 0.0;
+			}
+		};
+	};
+
+	template< typename Grammar >
+	struct left_shift_transform : proto::transform< left_shift_transform< Grammar > >
+	{
+		template< typename Expr , typename State , typename Data >
+		struct impl : proto::transform_impl< Expr , State , Data >
+		{
+			typedef double result_type;
+
+			result_type operator ()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Right shift transform" << endl;
+				return Grammar()( proto::right( expr ) , state , data )
+						+ ( data.which == 0 ) ? proto::value( proto::left( expr ) ) : 0.0;
+			}
+		};
+	};
+
+
+	template< typename Grammar >
+	struct right_scalar_multiplies_transform : proto::transform< right_scalar_multiplies_transform< Grammar > >
+	{
+		template< typename Expr , typename State , typename Data >
+		struct impl : proto::transform_impl< Expr , State , Data >
+		{
+			typedef double result_type;
+
+			result_type operator ()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Right scalar multiplies transform" << endl;
+				return Grammar()( proto::left( expr ) , state , data ) * proto::value( proto::right( expr ) );
+			}
+		};
+	};
+
+	template< typename Grammar >
+	struct left_scalar_multiplies_transform : proto::transform< left_scalar_multiplies_transform< Grammar > >
+	{
+		template< typename Expr , typename State , typename Data >
+		struct impl : proto::transform_impl< Expr , State , Data >
+		{
+			typedef double result_type;
+
+			result_type operator ()(
+					typename impl::expr_param expr ,
+					typename impl::state_param state ,
+					typename impl::data_param data ) const
+			{
+//				cout << "Left scalar multiplies transform" << endl;
+				return Grammar()( proto::right( expr ) , state , data ) * proto::value( proto::left( expr ) );
+			}
+		};
+	};
+
+
+
+
+
+	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 > >  { };
+
+	template< typename Grammar >
+	struct taylor_shift :
+		proto::or_<
+			proto::when< proto::plus< Grammar , proto::terminal< double > > , right_shift_transform< Grammar > > ,
+			proto::when< proto::plus< proto::terminal< double > , Grammar > , left_shift_transform< Grammar > >
+		> { };
+
+	template< typename Grammar >
+	struct taylor_scalar_multiplies :
+		proto::or_<
+			proto::when< proto::multiplies< Grammar , proto::terminal< double > > , right_scalar_multiplies_transform< Grammar > > ,
+			proto::when< proto::multiplies< proto::terminal< double > , Grammar > , left_scalar_multiplies_transform< Grammar > >
+		> { };
+
+
+
+
+
+	struct taylor_transform :
+	proto::or_
+	<
+		taylor_shift< taylor_transform > ,
+		taylor_scalar_multiplies< taylor_transform > ,
+		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 State state_type;
+		typedef taylor_context< state_type , Order > taylor_context_type;
+		typedef typename taylor_context_type::deriv_type deriv_type;
+
+		System m_sys;
+		taylor_context_type m_data;
+
+		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 )
+		{
+			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;
+		}
+	};
+
+	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 , double dt )
+	{
+		return eval_derivs< System , State , Order >( sys , x , derivs , i , dt );
+	}
+
+
+	struct optimize_tree : proto::or_< proto::_ > { };
+}
+
+template< size_t N , size_t Order >
+class taylor
+{
+public:
+
+	static const size_t dim = N;
+	static const size_t order_value = Order;
+	typedef double value_type;
+	typedef value_type time_type;
+	typedef unsigned short order_type;
+
+	typedef std::tr1::array< value_type , dim > state_type;
+	typedef state_type deriv_type;
+	typedef std::tr1::array< state_type , order_value > derivs_type;
+
+    order_type order( void ) const
+    {
+    	return order_value;
+    }
+
+    order_type stepper_order( void ) const
+    {
+    	return order_value;
+    }
+
+    order_type error_order( void ) const
+    {
+    	return order_value - 1;
+    }
+
+
+
+
+	template< class System >
+	void do_step( System sys , state_type &x , time_type t , time_type dt )
+	{
+		BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
+		BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+
+		do_step( sys , x , t , x , dt );
+	}
+
+	template< class System >
+	void do_step( System sys , const state_type &in , time_type t , state_type &out , time_type dt )
+	{
+		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( optimized_system , in , m_derivs , i , dt ) );
+		}
+
+		for( size_t i=0 ; i<N ; ++i )
+		{
+			out[i] = in[i];
+			for( size_t k=0 ; k<order_value ; ++k )
+				out[i] += m_derivs[k][i];
+		}
+	}
+
+	template< class System >
+	void do_step( System sys , state_type &x , time_type t , time_type dt , state_type &xerr )
+	{
+		BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
+		BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+		BOOST_STATIC_ASSERT(( order_value > 1 ));
+
+		do_step( sys , x , t , x , dt , xerr );
+	}
+
+	template< class System >
+	void do_step( System sys , const state_type &in , time_type t , state_type &out , time_type dt , state_type &xerr )
+	{
+		BOOST_STATIC_ASSERT(( boost::fusion::traits::is_sequence< System >::value ));
+		BOOST_STATIC_ASSERT(( size_t( boost::fusion::result_of::size< System >::value ) == dim ));
+		BOOST_STATIC_ASSERT(( order_value > 1 ));
+
+		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 ) );
+		}
+
+		for( size_t i=0 ; i<N ; ++i )
+		{
+			out[i] = in[i];
+			for( size_t k=0 ; k<order_value ; ++k )
+				out[i] += m_derivs[k][i];
+			xerr[i] = m_derivs[order_value-1][i];
+		}
+	}
+
+	const derivs_type& get_last_derivs( void ) const
+	{
+		return m_derivs;
+	}
+
+
+
+private:
+
+	derivs_type m_derivs;
+
+};
+
+
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif /* BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_HPP_ */

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_controller.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_controller.hpp	2011-04-16 17:03:53 EDT (Sat, 16 Apr 2011)
@@ -0,0 +1,148 @@
+/*
+ boost header: NUMERIC_ODEINT_STEPPER/controlled_error_stepper.hpp
+
+ Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
+ Copyright 2009 Andre Bergner
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_CONTROLLER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_CONTROLLER_HPP_INCLUDED
+
+#include <cmath>
+
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+#include <boost/numeric/odeint/stepper/controlled_error_stepper.hpp>
+
+
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+/*
+ * explicit stepper version
+ */
+template
+<
+	class TaylorStepper
+>
+class taylor_controller
+{
+
+public:
+
+	typedef TaylorStepper stepper_type;
+	typedef typename stepper_type::state_type state_type;
+	typedef typename stepper_type::value_type value_type;
+	typedef typename stepper_type::deriv_type deriv_type;
+	typedef typename stepper_type::time_type time_type;
+	typedef typename stepper_type::order_type order_type;
+	typedef default_error_checker< value_type > error_checker_type;
+	typedef controlled_stepper_tag stepper_category;
+
+
+	taylor_controller(
+			value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
+			value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
+			value_type a_x = static_cast< value_type >( 1.0 ) ,
+			value_type a_dxdt = static_cast< value_type >( 1.0 ) ,
+			const stepper_type &stepper = stepper_type()
+			)
+	: m_stepper( stepper ) ,
+	  m_error_checker( eps_abs , eps_rel , a_x , a_dxdt ) ,
+	  m_xerr() , m_xnew()
+	{
+	}
+
+
+	template< class System >
+	controlled_step_result try_step( System system , state_type &x , time_type &t , time_type &dt )
+	{
+		controlled_step_result res = try_step( system , x , t , m_xnew , dt );
+		if( ( res == success_step_size_increased ) || ( res == success_step_size_unchanged ) )
+		{
+			x = m_xnew;
+		}
+		return res;
+	}
+
+
+	template< class System >
+	controlled_step_result try_step( System system , state_type &in , time_type &t , state_type &out , time_type &dt )
+	{
+		using std::max;
+		using std::min;
+		using std::pow;
+
+		// do one step with error calculation
+		m_stepper.do_step( system , in , t , out , dt , m_xerr );
+
+		value_type max_rel_error = m_error_checker.error( in , m_stepper.get_last_derivs()[0] , m_xerr , dt );
+
+		if( max_rel_error > 1.1 )
+		{
+			// error too large - decrease dt ,limit scaling factor to 0.2 and reset state
+			dt *= max( 0.9 * pow( max_rel_error , -1.0 / ( m_stepper.error_order() - 1.0 ) ) , 0.2 );
+			return step_size_decreased;
+		}
+		else
+		{
+			if( max_rel_error < 0.5 )
+			{
+				//error too small - increase dt and keep the evolution and limit scaling factor to 5.0
+				t += dt;
+				dt *= min( 0.9 * pow( max_rel_error , -1.0 / m_stepper.stepper_order() ) , 5.0 );
+				return success_step_size_increased;
+			}
+			else
+			{
+				t += dt;
+				return success_step_size_unchanged;
+			}
+		}
+	}
+
+
+	stepper_type& stepper( void )
+	{
+		return m_stepper;
+	}
+
+	const stepper_type& stepper( void ) const
+	{
+		return m_stepper;
+	}
+
+
+private:
+
+
+
+	stepper_type m_stepper;
+	error_checker_type m_error_checker;
+	state_type m_xerr;
+	state_type m_xnew;
+};
+
+
+
+
+
+
+
+
+} // odeint
+} // numeric
+} // boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_STEPPER_TAYLOR_CONTROLLER_HPP_INCLUDED

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_lorenz_evol_adaptive_statistics.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_lorenz_evol_adaptive_statistics.cpp	2011-04-16 17:03:53 EDT (Sat, 16 Apr 2011)
@@ -0,0 +1,166 @@
+/*
+ * main.cpp
+ *
+ *  Created on: Apr 2, 2011
+ *      Author: karsten
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include "taylor.hpp"
+#include "taylor_controller.hpp"
+
+#include <boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp>
+#include <boost/numeric/odeint/stepper/controlled_error_stepper.hpp>
+#include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
+
+#include <boost/fusion/include/make_vector.hpp>
+#include <boost/spirit/include/phoenix.hpp>
+#include <boost/assign.hpp>
+#include <boost/timer.hpp>
+
+#include <boost/mpl/range_c.hpp>
+#include <boost/mpl/for_each.hpp>
+
+#define tab "\t"
+
+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 std::tr1::array< double , 3 > state_type;
+
+typedef boost::numeric::odeint::explicit_error_rk54_ck< state_type > rk54_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;
+namespace phoenix = boost::phoenix;
+namespace mpl = boost::mpl;
+
+using namespace std;
+using namespace boost::numeric::odeint;
+using namespace boost::assign;
+
+using boost::numeric::odeint::taylor_adf::arg1;
+using boost::numeric::odeint::taylor_adf::arg2;
+using boost::numeric::odeint::taylor_adf::arg3;
+
+struct run
+{
+	vector< double > m_eps_abs_values;
+	vector< double > m_eps_rel_values;
+	double m_t_end;
+
+	run( const vector< double > &eps_abs_values , const vector< double > &eps_rel_values , double t_end )
+	: m_eps_abs_values( eps_abs_values ) , m_eps_rel_values( eps_rel_values ) , m_t_end( t_end ){ }
+
+	template< class Order >
+	void operator()( Order ) const
+	{
+		static const size_t order = Order::value;
+		typedef boost::numeric::odeint::taylor< 3 , order > taylor_type;
+
+		boost::timer timer;
+
+		char str[512] = "";
+		sprintf( str , "dat/lorenz_taylor_%02d.dat" , int( order ) );
+		ofstream fout( str );
+
+		for( size_t i=0 ; i<m_eps_abs_values.size() ; ++i )
+		{
+			for( size_t j=0 ; j<m_eps_rel_values.size() ; ++j )
+			{
+				double eps_abs = m_eps_abs_values[i];
+				double eps_rel = m_eps_rel_values[j];
+
+				taylor_controller< taylor_type > taylor_controller( eps_abs , eps_rel );
+
+				state_type x = {{ 10.0 , 10.0 , 10.0 }};
+
+				timer.restart();
+				size_t steps_taylor = integrate_adaptive( taylor_controller ,
+						fusion::make_vector
+						(
+								sigma * ( arg2 - arg1 ) ,
+								R * arg1 - arg2 - arg1 * arg3 ,
+								arg1 * arg2 - b * arg3
+						) , x , 0.0 , m_t_end , 0.1 );
+				double time_taylor = timer.elapsed();
+
+				fout << i << tab << j << tab << eps_abs << tab << eps_rel << tab;
+				fout << steps_taylor << tab << time_taylor;
+				fout << endl;
+			}
+			fout << endl;
+		}
+
+
+	}
+};
+
+
+
+
+int main( int argc , char **argv )
+{
+	if( argc != 2 )
+	{
+		cerr << "usage taylor_lorenz_eval_adaptive_statistics t_end" << endl;
+		return -1;
+	}
+	double t_end = atof( argv[1] );
+
+	vector< double > eps_abs_values;
+	vector< double > eps_rel_values;
+
+	eps_abs_values += 1.0e1 , 1.0 , 1.0e-1 , 1.0e-2 , 1.0e-3 , 1.0e-4 , 1.0e-5 , 1.0e-6 , 1.0e-7 , 1.0e-8 , 1.0e-9 , 1.0e-10 , 1.0e-11 , 1.0e-12 , 1.0e-13 , 1.0e-14;
+	eps_rel_values += 1.0e-1 , 1.0e-2 , 1.0e-3 , 1.0e-4 , 1.0e-5 , 1.0e-6 , 1.0e-7 , 1.0e-8 , 1.0e-9 , 1.0e-10 , 1.0e-11 , 1.0e-12 , 1.0e-13 , 1.0e-14;
+
+
+	typedef mpl::range_c< size_t , 5 , 30 > order_values;
+	mpl::for_each< order_values >( run( eps_abs_values , eps_rel_values , t_end ) );
+
+
+	boost::timer timer;
+	ofstream fout( "dat/lorenz_rk54.dat" );
+	for( size_t i=0 ; i<eps_abs_values.size() ; ++i )
+	{
+		for( size_t j=0 ; j<eps_rel_values.size() ; ++j )
+		{
+			double eps_abs = eps_abs_values[i];
+			double eps_rel = eps_rel_values[j];
+
+			rk54_type rk54_plain;
+			controlled_error_stepper< rk54_type > rk54( rk54_plain , default_error_checker< double >( eps_abs , eps_rel ) );
+
+			state_type x = {{ 10.0 , 10.0 , 10.0 }};
+
+			timer.restart();
+			size_t steps_rk54 = integrate_adaptive( rk54 , lorenz , x , 0.0 , t_end , 0.1 );
+			double time_rk54 = timer.elapsed();
+
+			fout << i << tab << j << tab << eps_abs << tab << eps_rel << tab;
+			fout << steps_rk54 << tab << time_rk54 << tab;
+			fout << endl;
+		}
+		fout << endl;
+	}
+
+	return 0;
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_lorenz_evol_direct.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_lorenz_evol_direct.cpp	2011-04-16 17:03:53 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;
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_main.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/taylor/taylor_v3/taylor_main.cpp	2011-04-16 17:03:53 EDT (Sat, 16 Apr 2011)
@@ -0,0 +1,83 @@
+/*
+ * main.cpp
+ *
+ *  Created on: Apr 2, 2011
+ *      Author: karsten
+ */
+
+#include <iostream>
+
+#include "taylor.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 , 10 > taylor_type;
+typedef taylor_type::state_type state_type;
+typedef taylor_type::derivs_type derivs_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 )
+{
+	cout.precision( 14 );
+
+	taylor_type t;
+
+	state_type in = {{ 10.0 , 10.0 , 10.0 }} , dxdt = {{ 0.0 , 0.0 , 0.0 }} , xerr = {{ 0.0 , 0.0 , 0.0 }} , out = {{ 0.0 ,0.0 , 0.0 }};
+
+	lorenz( in , dxdt , 0.0 );
+
+	cout << in << endl;
+	cout << dxdt << endl << endl;
+
+	t.do_step(
+			fusion::make_vector
+			(
+					sigma * ( arg2 - arg1 ) ,
+					R * arg1 - arg2 - arg1 * arg3 ,
+					arg1 * arg2 - b * arg3
+			) ,
+			in , 0.0 , out , 0.1 , xerr );
+
+
+
+	cout << endl;
+	cout << in << endl;
+	cout << xerr << endl;
+	cout << out << endl << endl;
+	const derivs_type &derivs = t.get_last_derivs();
+	for( size_t i=0 ; i<derivs.size() ; ++i )
+		cout << derivs[i] << endl;
+
+	return 0;
+}