Merge branch 'master' of github.com:headmyshoulder/odeint-v2

boost/numeric/odeint/stepper/bulirsch_stoer.hpp

@@ -197,13 +197,13 @@ public:
         value_vector h_opt( m_k_max+1 );
         value_vector work( m_k_max+1 );
 
-        //std::cout << "t=" << t <<", dt=" << dt << "(" << m_dt_last << ")" << ", k_opt=" << m_current_k_opt << std::endl;
+        std::cout << "t=" << t <<", dt=" << dt << "(" << m_dt_last << ")" << ", k_opt=" << m_current_k_opt << std::endl;
 
         time_type new_h = dt;
 
         for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ )
         {
-            //std::cout << "k=" << k <<": " << ", first: " << m_first << std::endl;
+            std::cout << "  k=" << k; //<<": " << ", first: " << m_first << std::endl;
             m_midpoint.set_steps( m_interval_sequence[k] );
             if( k == 0 )
             {
@@ -220,7 +220,7 @@ public:
                 h_opt[k] = calc_h_opt( dt , error , k );
                 work[k] = m_cost[k]/h_opt[k];
                 //std::cout << '\t' << "h_opt=" << h_opt[k] << ", work=" << work[k] << std::endl;
-                //std::cout << '\t' << "error: " << error << std::endl;
+                std::cout << '\t' << "error: " << error << std::endl;
 
                 if( (k == m_current_k_opt-1) || m_first )
                 { // convergence before k_opt ?
@@ -285,7 +285,6 @@ public:
                         new_h = h_opt[m_current_k_opt];
                     } else
                     {
-                        //std::cout << "REJECT!" << std::endl;
                         reject = true;
                         new_h = h_opt[m_current_k_opt];
                     }
@@ -295,7 +294,10 @@ public:
         }
 
         if( !reject )
+        {
             t += dt;
+        } else
+            std::cout << "REJECT!" << std::endl;
 
         if( !m_last_step_rejected || (new_h < dt) )
         {

libs/numeric/odeint/examples/2d_lattice/Jamfile

@@ -0,0 +1,14 @@
+# Copyright 2011 Karsten Ahnert and 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 
+      <include>../../../../..
+      <define>BOOST_ALL_NO_LIB=1
+    : build-dir .
+    : default-build release
+    ;
+    
+exe spreading : spreading.cpp ;

libs/numeric/odeint/examples/2d_lattice/lattice2d.hpp

@@ -0,0 +1,155 @@
+/* stronlgy nonlinear hamiltonian lattice in 2d */
+
+#ifndef LATTICE2D_HPP
+#define LATTICE2D_HPP
+
+#include <vector>
+
+#include <boost/math/special_functions/pow.hpp>
+
+using boost::math::pow;
+
+template< int Kappa , int Lambda >
+struct lattice2d {
+
+    const double m_beta;
+    std::vector< std::vector< double > > m_omega;
+
+    lattice2d( const double beta )
+	: m_beta( beta )
+    { }
+
+    template< class StateIn , class StateOut >
+    void operator()( const StateIn &q , StateOut &dpdt )
+    {
+        // q and dpdt are 2d
+        const int N = q.size();
+
+        int i;
+        for( i = 0 ; i < N ; ++i )
+        {
+            const int i_l = (i-1+N) % N;
+            const int i_r = (i+1) % N;
+            for( int j = 0 ; j < N ; ++j )
+            {
+            const int j_l = (j-1+N) % N;
+            const int j_r = (j+1) % N;
+            dpdt[i][j] = - m_omega[i][j] * pow<Kappa-1>( q[i][j] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i][j_l] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i][j_r] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i_l][j] )
+                - m_beta * pow<Lambda-1>( q[i][j] - q[i_r][j] );
+            }
+        }
+    }
+
+    template< class StateIn >
+    double energy( const StateIn &q , const StateIn &p )
+    {
+        // q and dpdt are 2d
+        const int N = q.size();
+        double energy = 0.0;
+        int i;
+        for( i = 0 ; i < N ; ++i )
+        {
+            const int i_l = (i-1+N) % N;
+            const int i_r = (i+1) % N;
+            for( int j = 0 ; j < N ; ++j )
+            {
+            const int j_l = (j-1+N) % N;
+            const int j_r = (j+1) % N;
+            energy += p[i][j]*p[i][j] / 2.0
+                        + m_omega[i][j] * pow<Kappa>( q[i][j] ) / Kappa
+                + m_beta * pow<Lambda>( q[i][j] - q[i][j_l] ) / Lambda / 2
+                + m_beta * pow<Lambda>( q[i][j] - q[i][j_r] ) / Lambda / 2
+                + m_beta * pow<Lambda>( q[i][j] - q[i_l][j] ) / Lambda / 2
+                + m_beta * pow<Lambda>( q[i][j] - q[i_r][j] ) / Lambda / 2;
+            }
+        }
+        return energy;
+    }
+
+
+    template< class StateIn , class StateOut >
+    double local_energy( const StateIn &q , const StateIn &p , StateOut &energy )
+    {
+        // q and dpdt are 2d
+        const int N = q.size();
+        double e = 0.0;
+        int i;
+        for( i = 0 ; i < N ; ++i )
+        {
+            const int i_l = (i-1+N) % N;
+            const int i_r = (i+1) % N;
+            for( int j = 0 ; j < N ; ++j )
+            {
+                const int j_l = (j-1+N) % N;
+                const int j_r = (j+1) % N;
+                energy[i][j] = p[i][j]*p[i][j] / 2.0
+                    + m_omega[i][j] * pow<Kappa>( q[i][j] ) / Kappa
+                    + m_beta * pow<Lambda>( q[i][j] - q[i][j_l] ) / Lambda / 2
+                    + m_beta * pow<Lambda>( q[i][j] - q[i][j_r] ) / Lambda / 2
+                    + m_beta * pow<Lambda>( q[i][j] - q[i_l][j] ) / Lambda / 2
+                    + m_beta * pow<Lambda>( q[i][j] - q[i_r][j] ) / Lambda / 2;
+                e += energy[i][j];
+            }
+        }
+        //rescale
+        e = 1.0/e;
+        for( i = 0 ; i < N ; ++i )
+            for( int j = 0 ; j < N ; ++j )
+                energy[i][j] *= e;
+        return 1.0/e;
+    }
+
+    void load_pot( const char* filename , const double W , const double gap , 
+                   const size_t dim )
+    {
+        std::ifstream in( filename , std::ios::in | std::ios::binary );
+        if( !in.is_open() ) {
+            std::cerr << "pot file not found: " << filename << std::endl;
+            exit(0);
+        } else {
+            std::cout << "using pot file: " << filename << std::endl;
+        }
+
+        m_omega.resize( dim );
+        for( int i = 0 ; i < dim ; ++i )
+        {
+            m_omega[i].resize( dim );
+            for( size_t j = 0 ; j < dim ; ++j )
+            {
+                if( !in.good() )
+                {
+                    std::cerr << "I/O Error: " << filename << std::endl;
+                    exit(0);
+                }
+                double d;
+                in.read( (char*) &d , sizeof(d) );
+                if( (d < 0) || (d > 1.0) )
+                {
+                    std::cerr << "ERROR: " << d << std::endl;
+                    exit(0);
+                }
+                m_omega[i][j] = W*d + gap;
+            }
+        }
+
+    }
+
+    void generate_pot( const double W , const double gap , const size_t dim )
+    {
+        m_omega.resize( dim );
+        for( size_t i = 0 ; i < dim ; ++i )
+        {
+            m_omega[i].resize( dim );
+            for( size_t j = 0 ; j < dim ; ++j )
+            {
+                m_omega[i][j] = W*static_cast<double>(rand())/RAND_MAX + gap;
+            }
+        }
+    }
+
+};
+
+#endif

libs/numeric/odeint/examples/2d_lattice/nested_range_algebra.hpp

@@ -0,0 +1,36 @@
+/* nested range algebra */
+
+#ifndef NESTED_RANGE_ALGEBRA
+#define NESTED_RANGE_ALGEBRA
+
+namespace detail {
+
+    template< class Iterator1 , class Iterator2 , class Iterator3 , class Operation , class Algebra >
+    void for_each3( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3, Operation op , Algebra &algebra )
+{
+    for( ; first1 != last1 ; )
+        algebra.for_each3( *first1++ , *first2++ , *first3++ , op );
+}
+}
+
+
+template< class InnerAlgebra >
+struct nested_range_algebra
+{
+    
+    nested_range_algebra() 
+	: m_inner_algebra()
+    { }
+
+    template< class S1 , class S2 , class S3 , class Op >
+    void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op )
+    {
+        detail::for_each3( boost::begin( s1 ) , boost::end( s1 ) , boost::begin( s2 ) ,	boost::begin( s3 ) , op	, m_inner_algebra );
+    }
+
+
+private:
+    InnerAlgebra m_inner_algebra;
+};
+
+#endif

libs/numeric/odeint/examples/2d_lattice/spreading.cpp

@@ -0,0 +1,112 @@
+/*
+ * Example of a 2D simulation of nonlinearly coupled oscillators.
+ * Program just prints final energy which should be close to the initial energy (1.0).
+ * No parallelization is employed here.
+ * Run time on a 2.3GHz Intel Core-i5: about 10 seconds for 100 steps.
+ * Compile simply via bjam or directly:
+ * g++ -O3 -I${BOOST_ROOT} -I../../../../.. spreading.cpp
+ */
+
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+#include <sys/time.h>
+
+#include <boost/ref.hpp>
+#include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp>
+
+// we use a vector< vector< double > > as state type,
+// for that some functionality has to be added for odeint to work
+#include "nested_range_algebra.hpp"
+#include "vector_vector_resize.hpp"
+
+// defines the rhs of our dynamical equation
+#include "lattice2d.hpp"
+/* dynamical equations (Hamiltonian structure):
+dqdt_{i,j} = p_{i,j}
+dpdt_{i,j} = - omega_{i,j}*q_{i,j} - \beta*[ (q_{i,j} - q_{i,j-1})^3
+                                            +(q_{i,j} - q_{i,j+1})^3
+                                            +(q_{i,j} - q_{i-1,j})^3
+                                            +(q_{i,j} - q_{i+1,j})^3 ]
+*/
+
+
+using namespace std;
+
+static const int MAX_N = 1024;//2048;
+
+static const size_t KAPPA = 2;
+static const size_t LAMBDA = 4;
+static const double W = 1.0;
+static const double gap = 0.0;
+static const size_t steps = 100;
+static const double dt = 0.1;
+
+double initial_e = 1.0;
+double beta = 1.0;
+int realization_index = 0;
+
+//the state type
+typedef vector< vector< double > > state_type;
+
+//the stepper, choose a symplectic one to account for hamiltonian structure
+//use nested_range_algebra for calculations on vector< vector< ... > >
+typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan< 
+    state_type , state_type , double , state_type , state_type , double , 
+    nested_range_algebra< boost::numeric::odeint::range_algebra > ,
+    boost::numeric::odeint::default_operations > stepper_type;
+
+double time_diff_in_ms( timeval &t1 , timeval &t2 )
+{ return (t2.tv_sec - t1.tv_sec)*1000.0 + (t2.tv_usec - t1.tv_usec)/1000.0 + 0.5; }
+
+
+int main( int argc, const char* argv[] ) {
+
+    srand( time(NULL) );
+
+    lattice2d< KAPPA , LAMBDA > lattice( beta );
+
+
+    lattice.generate_pot( W , gap , MAX_N );
+
+    state_type q( MAX_N , vector< double >( MAX_N , 0.0 ) );
+
+    state_type p( q );
+
+    state_type energy( q );
+
+    p[MAX_N/2][MAX_N/2] = sqrt( 0.5*initial_e );
+    p[MAX_N/2+1][MAX_N/2] = sqrt( 0.5*initial_e );
+    p[MAX_N/2][MAX_N/2+1] = sqrt( 0.5*initial_e );
+    p[MAX_N/2+1][MAX_N/2+1] = sqrt( 0.5*initial_e );
+
+    cout.precision(10);
+
+    lattice.local_energy( q , p , energy );
+    double e=0.0;
+    for( size_t i=0 ; i<energy.size() ; ++i )
+        for( size_t j=0 ; j<energy[i].size() ; ++j )
+        {
+            e += energy[i][j];
+        }
+
+    cout << "initial energy: " << lattice.energy( q , p ) << endl;
+
+    timeval elapsed_time_start , elapsed_time_end;
+    gettimeofday(&elapsed_time_start , NULL);
+
+    stepper_type stepper;
+
+    for( size_t step=0 ; step<=steps ; ++step )
+    {
+        stepper.do_step( boost::ref( lattice ) , 
+                         make_pair( boost::ref( q ) , boost::ref( p ) ) , 
+                         0.0 , 0.1 );
+    }
+
+    gettimeofday(&elapsed_time_end , NULL);
+    double elapsed_time = time_diff_in_ms( elapsed_time_start , elapsed_time_end );
+    cout << steps << " steps in " << elapsed_time/1000 << " s (energy: " << lattice.energy( q , p ) << ")" << endl;
+}

libs/numeric/odeint/examples/2d_lattice/vector_vector_resize.hpp

@@ -0,0 +1,95 @@
+/* reserved vector */
+
+#ifndef VECTOR_VECTOR_RESIZE_HPP
+#define VECTOR_VECTOR_RESIZE_HPP
+
+#include <vector>
+
+#include <boost/range.hpp>
+
+namespace boost { namespace numeric { namespace odeint {
+
+template<>
+struct is_resizeable< std::vector< std::vector< double > > >
+{
+    typedef boost::true_type type;
+    const static bool value = type::value;
+};
+
+template<>
+struct same_size_impl< std::vector< std::vector< double > > , std::vector< std::vector< double > > >
+{
+    typedef std::vector< std::vector< double > > state_type;
+
+    static bool same_size( const state_type &x1 ,
+                           const state_type &x2 )
+    {
+        bool same = ( boost::size( x1 ) == boost::size( x2 ) );
+        if( !same )
+            return false;
+        typename state_type::const_iterator begin1 = boost::begin( x1 );
+        typename state_type::const_iterator begin2 = boost::begin( x2 );
+        while( begin1 != boost::end( x1 ) )
+            same &= ( boost::size( *begin1++ ) == boost::size( *begin2++ ) );
+        return same;
+    }
+};
+
+template<>
+struct resize_impl< std::vector< std::vector< double > > , std::vector< std::vector< double > > >
+{
+    typedef std::vector< std::vector< double > > state_type;
+
+    static void resize( state_type &x1 , const state_type &x2 )
+    {
+        x1.resize( boost::size( x2 ) );
+        typename state_type::iterator begin1 = boost::begin( x1 );
+        typename state_type::const_iterator begin2 = boost::begin( x2 );
+        while( begin1 != boost::end( x1 ) )
+            (*begin1++).resize( boost::size( *begin2++ ) );
+    }
+};
+
+template<>
+struct state_wrapper< std::vector< std::vector< double > > , true >
+{
+    typedef std::vector< std::vector< double > > state_type;
+    typedef state_wrapper< state_type > state_wrapper_type;
+    typedef boost::true_type is_resizeable;
+
+    state_type m_v;
+
+    template< class State >
+    bool same_size( const State &x )
+    {
+        bool same = ( boost::size( m_v ) == boost::size( x ) );
+        if( !same )
+            return false;
+        typename state_type::iterator begin1 = boost::begin( m_v );
+        typename State::const_iterator begin2 = boost::begin( x );
+        while( begin1 != boost::end( m_v ) )
+            same &= ( boost::size( *begin1++ ) == boost::size( *begin2++ ) );
+        return same;
+    }
+
+    template< class State >
+    bool resize( const State &x )
+    {
+        if( !same_size( x ) )
+        {
+            m_v.resize( boost::size( x ) );
+            typename state_type::iterator begin1 = boost::begin( m_v );
+            typename State::const_iterator begin2 = boost::begin( x );
+            while( begin1 != boost::end( m_v ) )
+                (*begin1++).resize( boost::size( *begin2++ ) );
+
+            return true;
+        } else
+            return false;
+    }
+ 
+};
+
+} } }
+
+#endif

libs/numeric/odeint/examples/Jamfile

@@ -28,6 +28,7 @@ exe list_lattice : list_lattice.cpp ;
 exe stepper_details : stepper_details.cpp ;
 exe my_vector : my_vector.cpp ;
 exe lorenz_point : lorenz_point.cpp ;
+exe bs_bug : bs_bug.cpp ;
 
 # build-project mtl ;
 # build-project ublas ;

libs/numeric/odeint/examples/bs_bug.cpp

@@ -0,0 +1,39 @@
+#include <iostream>
+#include <boost/array.hpp>
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::array< double , 3 > state_type;
+
+typedef boost::numeric::odeint::bulirsch_stoer< state_type > bulirsch_stoer_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];
+}
+
+int main() {
+    
+    const double eps_abs = 10.0;//1E-4;
+    const double eps_rel = 1E-4;
+    const double t_end = 10.0;
+
+    bulirsch_stoer_type bulirsch_stoer( eps_abs , eps_rel );
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+
+    cout.precision(17);
+
+    size_t steps = integrate_adaptive( bulirsch_stoer , lorenz , x , 0.0 , t_end , 0.1 );
+
+    cout << steps << " steps" << endl;
+
+    return 0;
+}