[Boost-commit] svn:boost r59730 - in sandbox/odeint: . boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/examples

Subject: [Boost-commit] svn:boost r59730 - in sandbox/odeint: . boost/numeric/odeint boost/numeric/odeint/detail libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2010-02-17 10:01:46

Author: mariomulansky
Date: 2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
New Revision: 59730
URL: http://svn.boost.org/trac/boost/changeset/59730

Log:
 cleaned up iterator algebra, stepsize controller now uses error checker
Text files modified: 
   sandbox/odeint/ToDo                                                        |    11 +++--                                   
   sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp              |     4 +-                                      
   sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp        |    57 +++++++++++-----------------            
   sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp            |    79 +++++++++++++++++++++------------------ 
   sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp             |    46 ++++++++++------------                  
   sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp             |     2                                         
   sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp                  |    18 ++++++--                                
   sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp              |    17 +++++++-                                
   sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp                 |     2                                         
   sandbox/odeint/libs/numeric/odeint/examples/test_container_and_stepper.cpp |    14 ++++---                                 
   10 files changed, 132 insertions(+), 118 deletions(-)

Modified: sandbox/odeint/ToDo
==============================================================================
--- sandbox/odeint/ToDo	(original)
+++ sandbox/odeint/ToDo	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -10,11 +10,12 @@
 
 Mario:
 * code cleanup:
-  * error_checker in controlled_stepper_standard
-  * all iterator functionality in iterator_algebra
-  * cleanup iterator algebra 
+  * error_checker in controlled_stepper_standard - DONE
+  * all iterator functionality in iterator_algebra - DONE (where I found it)
+  * cleanup iterator algebra - DONE
+  * add scale_sum_inplace for different iterators numbers
 
-* all std::vector should given explictly the allocator
+* all std::vector should given explictly the allocator - that requires 4 additional template parameters in stepper_rk_generic?
 
 * numeric_traits or something similar for pow, abs, norm, sqrt ...
 
@@ -34,7 +35,7 @@
 
 * boost::ublas  / funktioniert
 * mtl4, implementing real iterators  / funktioniert so gehackt, nicht submitten
-* blitz++ / funktioniert , nicht submitten
+* blitz++ / nicht submitten -- BUGGY, da copy constructor nur referenz erzeugt!
 * boost::multi_array / muss gemacht werden
 
 

Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_bs.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -57,7 +57,7 @@
     private:
 
         stepper_midpoint< container_type, time_type, traits_type > m_stepper_mp;
-        error_checker_standard< container_type, time_type > m_error_checker;
+        error_checker_standard< container_type, time_type , traits_type > m_error_checker;
         
         const unsigned short m_k_max;
 
@@ -172,7 +172,7 @@
                 //std::clog << "Error: " << k << '\t' << m_xerr[0] << '\t' << m_xerr[1] << std::endl;
                 if( k != 0 ) 
                 {
-                    value_type max_err = m_error_checker.get_max_error_ratio(m_xerr, m_x_scale);
+                    time_type max_err = m_error_checker.get_max_error_ratio(m_xerr, m_x_scale);
                     m_error[k-1] = std::pow( max_err/m_safety1, 1.0/(2*k+1) );
                     if( (k >= m_current_k_opt-1) || !continuous_call )
                     { //we're in the order window where convergence is expected

Modified: sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/controlled_stepper_standard.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -22,6 +22,7 @@
 
 #include <boost/numeric/odeint/detail/iterator_algebra.hpp>
 #include <boost/numeric/odeint/concepts/state_concept.hpp>
+#include <boost/numeric/odeint/error_checker_standard.hpp>
 #include <boost/numeric/odeint/container_traits.hpp>
 
 namespace boost {
@@ -78,6 +79,7 @@
     private:
 
         stepper_type &m_stepper;
+        error_checker_standard< container_type, time_type , traits_type > m_error_checker;
 
         time_type m_eps_abs;
         time_type m_eps_rel;
@@ -86,33 +88,7 @@
         container_type m_dxdt;
         container_type m_x_tmp;
         container_type m_x_err;
-
-
-        // private methods
-
-        time_type calc_max_rel_err(
-                iterator x_start ,
-                iterator x_end ,
-                iterator dxdt_start ,
-                iterator x_err_start ,
-                time_type dt
-            )
-        {
-	    time_type max_rel_err = 0.0;
-            time_type err;
-
-	    while( x_start != x_end )
-            {
-                // get the maximal value of x_err/D where 
-                // D = eps_abs + eps_rel * (a_x*|x| + a_dxdt*|dxdt|);
-                err = m_eps_abs + m_eps_rel * (
-                    m_a_x * std::abs(*x_start++) + 
-                    m_a_dxdt * dt * std::abs(*dxdt_start++) );
-                max_rel_err = std::max( std::abs(*x_err_start++)/err , max_rel_err );
-	    }
-            return max_rel_err;
-        }
-
+        container_type m_x_scale;
 
 
         // public functions
@@ -123,6 +99,7 @@
                 time_type abs_err, time_type rel_err, 
                 time_type factor_x, time_type factor_dxdt )
             : m_stepper(stepper), 
+              m_error_checker( abs_err, rel_err, factor_x, factor_dxdt ),
               m_eps_abs(abs_err),
               m_eps_rel(rel_err),
               m_a_x(factor_x),
@@ -142,20 +119,19 @@
         controlled_step_result try_step( 
                 DynamicalSystem &system, 
                 container_type &x, 
+                container_type &dxdt,
                 time_type &t, 
                 time_type &dt )
         {
             traits_type::adjust_size( x , m_x_err );
-            traits_type::adjust_size( x , m_dxdt );
+            traits_type::adjust_size( x , m_x_scale );
 
-	    m_x_tmp = x;
-	    system( x , m_dxdt , t ); 
-	    m_stepper.do_step( system , x , m_dxdt, t , dt , m_x_err );
+            m_error_checker.fill_scale( x , dxdt , dt , m_x_scale );
 
-            time_type max_rel_err = calc_max_rel_err(
-                m_x_tmp.begin() , m_x_tmp.end() ,
-                m_dxdt.begin() , m_x_err.begin() , dt );
+	    m_x_tmp = x;
+	    m_stepper.do_step( system , x , dxdt, t , dt , m_x_err );
 
+            time_type max_rel_err = m_error_checker.get_max_error_ratio(m_x_err, m_x_scale);
 
             if( max_rel_err > 1.1 )
             { 
@@ -185,6 +161,19 @@
                 }
             }
         }
+
+	template< class DynamicalSystem >
+	controlled_step_result try_step( 
+                DynamicalSystem &system, 
+                container_type &x, 
+                time_type &t, 
+                time_type &dt )
+	{
+            traits_type::adjust_size( x , m_dxdt );
+            system( x , m_dxdt , t );
+            return try_step( system , x , m_dxdt , t , dt );
+        }
+
     };
 
 } // namespace odeint

Modified: sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/detail/iterator_algebra.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -14,6 +14,8 @@
 #ifndef BOOST_NUMERIC_ODEINT_DETAIL_ACCUMULATORS_HPP
 #define BOOST_NUMERIC_ODEINT_DETAIL_ACCUMULATORS_HPP
 
+#include <cmath>
+
 #include <iostream>
 
 namespace boost {
@@ -40,42 +42,6 @@
             (*first1++) += alpha * (*first2++);
     }
 
-    // computes y = x1 - x2
-    template <
-        class OutputIterator ,
-        class InputIterator1 ,
-        class InputIterator2
-        >
-    void assign_diff(
-                     OutputIterator first1 ,
-                     OutputIterator last1 ,
-                     InputIterator1 first2 ,
-                     InputIterator2 first3 )
-    {
-        while( first1 != last1 )
-            (*first1++) = (*first2++) - (*first3++);
-    }
-
-
-    // computes y = x1 + alpha * x2
-    template <
-        class OutputIterator ,
-        class InputIterator1 ,
-        class InputIterator2 ,
-        class T
-        >
-    void assign_sum(
-                    OutputIterator first1 ,
-                    OutputIterator last1 ,
-                    InputIterator1 first2 ,
-                    InputIterator2 first3 ,
-                    T alpha )
-    {
-        while( first1 != last1 )
-            (*first1++) = (*first2++) + alpha * (*first3++);
-    }
-
-
 
     // computes y = alpha1 * ( x1 + x2 + alpha2*x3 )
     template <
@@ -561,6 +527,47 @@
           }
     }
 
+    template<
+        class OutputIterator,
+        class InputIterator1,
+        class InputIterator2,
+        class T >
+    void weighted_error( OutputIterator y_begin,
+                         OutputIterator y_end,
+                         InputIterator1 x1_begin,
+                         InputIterator2 x2_begin,
+                         T eps_abs,
+                         T eps_rel,
+                         T a_x,
+                         T a_dxdt )
+    {
+        while( y_begin != y_end ) 
+        {
+            *y_begin++ = eps_abs + 
+                eps_rel * (a_x * std::abs(*x1_begin++) + 
+                           a_dxdt * std::abs(*x2_begin++));
+        }
+    }
+
+    
+    template<
+        class InputIterator1,
+        class InputIterator2,
+        class T >
+    T max_ratio( InputIterator1 x1_begin,
+                 InputIterator1 x1_end,
+                 InputIterator2 x2_begin,
+                 T initial_max )
+    {
+        while( x1_begin != x1_end ) 
+        {
+            initial_max = std::max( static_cast<T>(std::abs(*x1_begin++)/(*x2_begin++)),
+                                    initial_max);
+        }
+        return initial_max;
+    }
+    
+
     
     
 

Modified: sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/error_checker_standard.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -23,15 +23,20 @@
 namespace numeric {
 namespace odeint {
 
-    template< class Container, class Time >
+    template< class Container, 
+              class Time, 
+              class Traits = container_traits< Container > >
     class error_checker_standard {
 
 
     public:
-        typedef Container container_type;
         typedef Time time_type;
-        typedef typename container_type::value_type value_type;
-        typedef typename container_type::iterator iterator;
+        typedef Traits traits_type;
+        typedef typename traits_type::container_type container_type;
+        typedef typename traits_type::value_type value_type;
+        typedef typename traits_type::iterator iterator;
+        typedef typename traits_type::const_iterator const_iterator;
+
 
     private:
                time_type m_eps_abs;
@@ -54,31 +59,22 @@
                 time_type dt, 
                 container_type &scale )
         {
-            iterator x_iter = x.begin();
-            const iterator x_end = x.end();
-            iterator dxdt_iter = dxdt.begin();
-            iterator scale_iter = scale.begin();
-            while( x_iter != x_end ) 
-            {
-                *scale_iter++ = m_eps_abs + 
-                    m_eps_rel * (m_a_x * std::abs(*x_iter++) + 
-                                 m_a_dxdt * dt * std::abs(*dxdt_iter++));
-            }
+            detail::it_algebra::weighted_error( traits_type::begin(scale),
+                                                traits_type::end(scale),
+                                                traits_type::begin(x),
+                                                traits_type::end(dxdt),
+                                                m_eps_abs,
+                                                m_eps_rel,
+                                                m_a_x,
+                                                m_a_x*dt );
         }
 
         time_type get_max_error_ratio( container_type &x_err, container_type &scale)
         {
-            iterator x_err_iter = x_err.begin();
-            const iterator x_err_end = x_err.end();
-            iterator scale_iter = scale.begin();
-            time_type max_rel_error = static_cast<time_type>(0.0);
-            while( x_err_iter != x_err_end ) 
-            {
-                max_rel_error = std::max(
-                        static_cast<time_type>(std::abs(*x_err_iter++)/(*scale_iter++)) , 
-                        max_rel_error);
-            }            
-            return max_rel_error;
+            return detail::it_algebra::max_ratio( traits_type::begin(x_err),
+                                                  traits_type::end(x_err),
+                                                  traits_type::begin(scale),
+                                                  static_cast<time_type>(0.0) );
         }
 
         const time_type get_epsilon() { return std::max(m_eps_rel, m_eps_abs); }

Modified: sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/hamiltonian_stepper_rk.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -91,7 +91,7 @@
 
             if( !traits_type::same_size( q , p ) )
             {
-		std::string msg( "hamiltonian_stepper_euler::do_step(): " );
+		std::string msg( "hamiltonian_stepper_rk::do_step(): " );
                 msg += "q and p have different sizes";
                 throw std::invalid_argument( msg );
             }

Modified: sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -18,6 +18,7 @@
 
 #include <boost/numeric/odeint/detail/iterator_algebra.hpp>
 
+#include <iostream>
 
 namespace boost {
 namespace numeric {
@@ -91,17 +92,24 @@
         {
             traits_type::adjust_size( x , xerr );
 
-            m_xtemp = x;
+            /*** BUG FOR BLITZ ARRAYS ***/
+            /* for blitzz arrays, the copy constructor creates a REFERENCE and not a copy!!! */
+            traits_type::adjust_size( x , m_xtemp );
+            m_xtemp = container_type(x);
+            /*** END BUG ***/
+
             time_type dt2 = static_cast<time_type>(0.5) * dt;
 
             do_step( system , m_xtemp , dxdt , t , dt );
             do_step( system , x , dxdt , t , dt2 );
             do_step( system , x , t+dt2 , dt2 );
 
-            detail::it_algebra::assign_diff( traits_type::begin(xerr) ,
-                                             traits_type::end(xerr) ,
-                                             traits_type::begin(m_xtemp) ,
-                                             traits_type::begin(x) );
+            detail::it_algebra::scale_sum( traits_type::begin(xerr) ,
+                                           traits_type::end(xerr) ,
+                                           static_cast< value_type >(1.0),
+                                           traits_type::begin(m_xtemp) ,
+                                           static_cast< value_type >(-1.0),
+                                           traits_type::begin(x) );
         }
 
 

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -84,13 +84,24 @@
             time_type th = t + dh;
 
             //m_xt = x + dh*dxdt
-            assign_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
-                        traits_type::begin(x) , traits_type::begin(dxdt) , dh );
+            /*assign_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
+              traits_type::begin(x) , traits_type::begin(dxdt) , dh );*/
+            scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
+                       static_cast< time_type >(1.0) ,
+                       traits_type::begin(x) , 
+                       dh, 
+                       traits_type::begin(dxdt) );
 
             system( m_xt , m_dxt , th );
             //m_xt = x + dh*m_dxdt
-            assign_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
+            /*assign_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
                         traits_type::begin(x) , traits_type::begin(m_dxt) , dh );
+            */
+            scale_sum( traits_type::begin(m_xt) , traits_type::end(m_xt) ,
+                        static_cast< time_type >(1.0),
+                        traits_type::begin(x) , 
+                        dh ,
+                        traits_type::begin(m_dxt) );
 
             system( m_xt , m_dxm , th );
             //m_xt = x + dt*m_dxm ; m_dxm += m_dxt

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -59,7 +59,7 @@
         // provide basic typedefs
     public:
 
-        typedef const unsigned short order_type;
+        typedef unsigned short order_type;
         typedef Container container_type;
         typedef Time time_type;
         typedef Traits traits_type;

Modified: sandbox/odeint/libs/numeric/odeint/examples/test_container_and_stepper.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/test_container_and_stepper.cpp	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/test_container_and_stepper.cpp	2010-02-17 10:01:45 EST (Wed, 17 Feb 2010)
@@ -179,7 +179,7 @@
     stepper_euler< state_type5 > s5;
 
     test_stepper( s1 , deriv1 , x1 , dxdt1 );
-
+    
 
     {
         stepper_euler< state_type1 > s1;
@@ -189,7 +189,7 @@
         stepper_euler< state_type5 > s5;
         test_steppers( s1 , s2 , s3 , s4 , s5 );
     }
-
+    cout << "Euler Stepper OK!" << endl;
 
 
     {
@@ -201,7 +201,7 @@
         test_steppers( s1 , s2 , s3 , s4 , s5 );
         test_error_steppers( s1 , s2 , s3 , s4 , s5 );
     }
-
+    cout << "Half Step Stepper OK!" << endl;
 
 
     {
@@ -212,7 +212,7 @@
         stepper_rk4< state_type5 > s5;
         test_steppers( s1 , s2 , s3 , s4 , s5 );
     }
-
+    cout << "RK4 Stepper OK!" << endl;
 
 
 
@@ -225,7 +225,7 @@
         stepper_rk4_classical< state_type5 > s5;
         test_steppers( s1 , s2 , s3 , s4 , s5 );
     }
-
+    cout << "RK4 Classical Stepper OK!" << endl;
 
 
 
@@ -237,6 +237,7 @@
         stepper_rk5_ck< state_type5 > s5;
         test_error_steppers( s1 , s2 , s3 , s4 , s5 );
     }
+    cout << "RK5 CK Stepper OK!" << endl;
 
 
     {
@@ -247,7 +248,7 @@
         stepper_midpoint< state_type5 > s5;
         test_steppers( s1 , s2 , s3 , s4 , s5 );
     }
-
+    cout << "Midpoint Stepper OK!" << endl;
 
 
     {
@@ -258,6 +259,7 @@
         stepper_rk78_fehlberg< state_type5 > s5;
         test_steppers( s1 , s2 , s3 , s4 , s5 );
     }
+    cout << "RK78 Stepper OK!" << endl;
 
 
     cout << "Everything is allright!" << endl;