1 //==============================================================================
2 // Copyright 2011-2014 Karsten Ahnert
3 // Copyright 2011-2014 Mario Mulansky
4 // Copyright 2014 LRI UMR 8623 CNRS/Univ Paris Sud XI
5 // Copyright 2014 NumScale SAS
6 //
7 // Distributed under the Boost Software License, Version 1.0.
8 // See accompanying file LICENSE.txt or copy at
9 // http://www.boost.org/LICENSE_1_0.txt
10 //==============================================================================
11
12 #include <iostream>
13 #include <utility>
14
15 #include <boost/numeric/odeint.hpp>
16
17 #ifndef M_PI //not there on windows
18 #define M_PI 3.141592653589793 //...
19 #endif
20
21 #include <boost/random.hpp>
22 #include <boost/dispatch/meta/as_integer.hpp>
23
24 #include <nt2/include/functions/cos.hpp>
25 #include <nt2/include/functions/sin.hpp>
26 #include <nt2/include/functions/atan2.hpp>
27 #include <nt2/table.hpp>
28 #include <nt2/include/functions/zeros.hpp>
29 #include <nt2/include/functions/sum.hpp>
30 #include <nt2/include/functions/mean.hpp>
31 #include <nt2/arithmetic/include/functions/hypot.hpp>
32 #include <nt2/include/functions/tie.hpp>
33
34 #include <boost/numeric/odeint/external/nt2/nt2_algebra_dispatcher.hpp>
35
36
37 using namespace std;
38 using namespace boost::numeric::odeint;
39
40 template <typename container_type, typename T>
calc_mean_field(const container_type & x)41 pair< T, T > calc_mean_field( const container_type &x )
42
43 {
44 T cos_sum = 0.0 , sin_sum = 0.0;
45
46 nt2::tie(cos_sum,sin_sum) = nt2::tie(nt2::mean( nt2::cos(x) ), nt2::mean( nt2::sin(x) ));
47
48 T K = nt2::hypot(sin_sum,cos_sum);
49 T Theta = nt2::atan2( sin_sum , cos_sum );
50
51 return make_pair( K , Theta );
52 }
53
54 template <typename container_type, typename T>
55 struct phase_ensemble
56 {
57 typedef typename boost::dispatch::meta::as_integer<T,unsigned>::type int_type;
58 container_type m_omega;
59 T m_epsilon;
60
phase_ensemblephase_ensemble61 phase_ensemble( const int_type n , T g = 1.0 , T epsilon = 1.0 )
62 : m_epsilon( epsilon )
63 {
64 m_omega = nt2::zeros(nt2::of_size(n), nt2::meta::as_<T>());
65 create_frequencies( g );
66 }
67
create_frequenciesphase_ensemble68 void create_frequencies( T g )
69 {
70 boost::mt19937 rng;
71 boost::cauchy_distribution<> cauchy( 0.0 , g );
72 boost::variate_generator< boost::mt19937&, boost::cauchy_distribution<> > gen( rng , cauchy );
73 generate( m_omega.begin() , m_omega.end() , gen );
74 }
75
set_epsilonphase_ensemble76 void set_epsilon( T epsilon ) { m_epsilon = epsilon; }
77
get_epsilonphase_ensemble78 T get_epsilon( void ) const { return m_epsilon; }
79
operator ()phase_ensemble80 void operator()( const container_type &x , container_type &dxdt , T ) const
81 {
82 pair< T, T > mean = calc_mean_field<container_type,T>( x );
83 dxdt = m_omega + m_epsilon * mean.first * nt2::sin( mean.second - x );
84 }
85 };
86
87 template<typename T>
88 struct statistics_observer
89 {
90 typedef typename boost::dispatch::meta::as_integer<T,unsigned>::type int_type;
91 T m_K_mean;
92 int_type m_count;
93
statistics_observerstatistics_observer94 statistics_observer( void )
95 : m_K_mean( 0.0 ) , m_count( 0 ) { }
96
97 template< class State >
operator ()statistics_observer98 void operator()( const State &x , T t )
99 {
100 pair< T, T > mean = calc_mean_field<State,T>( x );
101 m_K_mean += mean.first;
102 ++m_count;
103 }
104
get_K_meanstatistics_observer105 T get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / T( m_count ) : 0.0 ; }
106
resetstatistics_observer107 void reset( void ) { m_K_mean = 0.0; m_count = 0; }
108 };
109
110 template<typename T>
111 struct test_ode_table
112 {
113 typedef nt2::table<T> array_type;
114 typedef void experiment_is_immutable;
115
116 typedef typename boost::dispatch::meta::as_integer<T,unsigned>::type int_type;
117
test_ode_tabletest_ode_table118 test_ode_table ( )
119 : size_(16384), ensemble( size_ , 1.0 ), unif( 0.0 , 2.0 * M_PI ), gen( rng , unif ), obs()
120 {
121 x.resize(nt2::of_size(size_));
122 }
123
operator ()test_ode_table124 void operator()()
125 {
126 for( T epsilon = 0.0 ; epsilon < 5.0 ; epsilon += 0.1 )
127 {
128 ensemble.set_epsilon( epsilon );
129 obs.reset();
130
131 // start with random initial conditions
132 generate( x.begin() , x.end() , gen );
133 // calculate some transients steps
134 integrate_const( runge_kutta4< array_type, T >() , boost::ref( ensemble ) , x , T(0.0) , T(10.0) , dt );
135
136 // integrate and compute the statistics
137 integrate_const( runge_kutta4< array_type, T >() , boost::ref( ensemble ) , x , T(0.0) , T(100.0) , dt , boost::ref( obs ) );
138 cout << epsilon << "\t" << obs.get_K_mean() << endl;
139 }
140 }
141
operator <<(std::ostream & os,test_ode_table<T> const & p)142 friend std::ostream& operator<<(std::ostream& os, test_ode_table<T> const& p)
143 {
144 return os << "(" << p.size() << ")";
145 }
146
sizetest_ode_table147 std::size_t size() const { return size_; }
148
149 private:
150 std::size_t size_;
151 phase_ensemble<array_type,T> ensemble;
152 boost::uniform_real<> unif;
153 array_type x;
154 boost::mt19937 rng;
155 boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen;
156 statistics_observer<T> obs;
157
158 static const T dt = 0.1;
159 };
160
main()161 int main()
162 {
163 std::cout<< " With T = [double] \n";
164 test_ode_table<double> test_double;
165 test_double();
166
167 std::cout<< " With T = [float] \n";
168 test_ode_table<float> test_float;
169 test_float();
170 }
171