libs/numeric/odeint/examples/thrust/phase_oscillator_ensemble.cu

   1 /*
   2  The example how the phase_oscillator ensemble can be implemented using CUDA and thrust
   3
   4  Copyright 2011-2013 Mario Mulansky
   5  Copyright 2011 Karsten Ahnert
   6
   7  Distributed under the Boost Software License, Version 1.0.
   8  (See accompanying file LICENSE_1_0.txt or
   9  copy at http://www.boost.org/LICENSE_1_0.txt)
  10  */
  11
  12 #include <iostream>
  13 #include <fstream>
  14 #include <cmath>
  15 #include <utility>
  16
  17 #include <thrust/device_vector.h>
  18 #include <thrust/reduce.h>
  19 #include <thrust/functional.h>
  20
  21 #include <boost/numeric/odeint.hpp>
  22 #include <boost/numeric/odeint/external/thrust/thrust.hpp>
  23
  24 #include <boost/timer.hpp>
  25 #include <boost/random/cauchy_distribution.hpp>
  26
  27 using namespace std;
  28
  29 using namespace boost::numeric::odeint;
  30
  31 /*
  32  * Sorry for that dirty hack, but nvcc has large problems with boost::random.
  33  *
  34  * Nevertheless we need the cauchy distribution from boost::random, and therefore
  35  * we need a generator. Here it is:
  36  */
  37 struct drand48_generator
  38 {
  39     typedef double result_type;
  40     result_type operator()( void ) const { return drand48(); }
  41     result_type min( void ) const { return 0.0; }
  42     result_type max( void ) const { return 1.0; }
  43 };
  44
  45 //[ thrust_phase_ensemble_state_type
  46 //change this to float if your device does not support double computation
  47 typedef double value_type;
  48
  49 //change this to host_vector< ... > of you want to run on CPU
  50 typedef thrust::device_vector< value_type > state_type;
  51 // typedef thrust::host_vector< value_type > state_type;
  52 //]
  53
  54
  55 //[ thrust_phase_ensemble_mean_field_calculator
  56 struct mean_field_calculator
  57 {
  58     struct sin_functor : public thrust::unary_function< value_type , value_type >
  59     {
  60         __host__ __device__
  61         value_type operator()( value_type x) const
  62         {
  63             return sin( x );
  64         }
  65     };
  66
  67     struct cos_functor : public thrust::unary_function< value_type , value_type >
  68     {
  69         __host__ __device__
  70         value_type operator()( value_type x) const
  71         {
  72             return cos( x );
  73         }
  74     };
  75
  76     static std::pair< value_type , value_type > get_mean( const state_type &x )
  77     {
  78         //[ thrust_phase_ensemble_sin_sum
  79         value_type sin_sum = thrust::reduce(
  80                 thrust::make_transform_iterator( x.begin() , sin_functor() ) ,
  81                 thrust::make_transform_iterator( x.end() , sin_functor() ) );
  82         //]
  83         value_type cos_sum = thrust::reduce(
  84                 thrust::make_transform_iterator( x.begin() , cos_functor() ) ,
  85                 thrust::make_transform_iterator( x.end() , cos_functor() ) );
  86
  87         cos_sum /= value_type( x.size() );
  88         sin_sum /= value_type( x.size() );
  89
  90         value_type K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum );
  91         value_type Theta = atan2( sin_sum , cos_sum );
  92
  93         return std::make_pair( K , Theta );
  94     }
  95 };
  96 //]
  97
  98
  99
 100 //[ thrust_phase_ensemble_sys_function
 101 class phase_oscillator_ensemble
 102 {
 103
 104 public:
 105
 106     struct sys_functor
 107     {
 108         value_type m_K , m_Theta , m_epsilon;
 109
 110         sys_functor( value_type K , value_type Theta , value_type epsilon )
 111         : m_K( K ) , m_Theta( Theta ) , m_epsilon( epsilon ) { }
 112
 113         template< class Tuple >
 114         __host__ __device__
 115         void operator()( Tuple t )
 116         {
 117             thrust::get<2>(t) = thrust::get<1>(t) + m_epsilon * m_K * sin( m_Theta - thrust::get<0>(t) );
 118         }
 119     };
 120
 121     // ...
 122     //<-
 123     phase_oscillator_ensemble( size_t N , value_type g = 1.0 , value_type epsilon = 1.0 )
 124         : m_omega() , m_N( N ) , m_epsilon( epsilon )
 125     {
 126         create_frequencies( g );
 127     }
 128
 129     void create_frequencies( value_type g )
 130     {
 131         boost::cauchy_distribution< value_type > cauchy( 0.0 , g );
 132 //        boost::variate_generator< boost::mt19937&, boost::cauchy_distribution< value_type > > gen( rng , cauchy );
 133         drand48_generator d48;
 134         vector< value_type > omega( m_N );
 135         for( size_t i=0 ; i<m_N ; ++i )
 136             omega[i] = cauchy( d48 );
 137 //        generate( omega.begin() , omega.end() , gen );
 138         m_omega = omega;
 139     }
 140
 141     void set_epsilon( value_type epsilon ) { m_epsilon = epsilon; }
 142
 143     value_type get_epsilon( void ) const { return m_epsilon; }
 144     //->
 145
 146     void operator() ( const state_type &x , state_type &dxdt , const value_type dt ) const
 147     {
 148         std::pair< value_type , value_type > mean_field = mean_field_calculator::get_mean( x );
 149
 150         thrust::for_each(
 151                 thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_omega.begin() , dxdt.begin() ) ),
 152                 thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_omega.end() , dxdt.end()) ) ,
 153                 sys_functor( mean_field.first , mean_field.second , m_epsilon )
 154                 );
 155     }
 156
 157     // ...
 158     //<-
 159 private:
 160
 161     state_type m_omega;
 162     const size_t m_N;
 163     value_type m_epsilon;
 164     //->
 165 };
 166 //]
 167
 168
 169 //[ thrust_phase_ensemble_observer
 170 struct statistics_observer
 171 {
 172     value_type m_K_mean;
 173     size_t m_count;
 174
 175     statistics_observer( void )
 176     : m_K_mean( 0.0 ) , m_count( 0 ) { }
 177
 178     template< class State >
 179     void operator()( const State &x , value_type t )
 180     {
 181         std::pair< value_type , value_type > mean = mean_field_calculator::get_mean( x );
 182         m_K_mean += mean.first;
 183         ++m_count;
 184     }
 185
 186     value_type get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / value_type( m_count ) : 0.0 ; }
 187
 188     void reset( void ) { m_K_mean = 0.0; m_count = 0; }
 189 };
 190 //]
 191
 192
 193
 194 // const size_t N = 16384 * 128;
 195 const size_t N = 16384;
 196 const value_type pi = 3.1415926535897932384626433832795029;
 197 const value_type dt = 0.1;
 198 const value_type d_epsilon = 0.1;
 199 const value_type epsilon_min = 0.0;
 200 const value_type epsilon_max = 5.0;
 201 const value_type t_transients = 10.0;
 202 const value_type t_max = 100.0;
 203
 204 int main( int arc , char* argv[] )
 205 {
 206     // initial conditions on host
 207     vector< value_type > x_host( N );
 208     for( size_t i=0 ; i<N ; ++i ) x_host[i] = 2.0 * pi * drand48();
 209
 210     //[ thrust_phase_ensemble_system_instance
 211     phase_oscillator_ensemble ensemble( N , 1.0 );
 212     //]
 213
 214
 215
 216     boost::timer timer;
 217     boost::timer timer_local;
 218     double dopri5_time = 0.0 , rk4_time = 0.0;
 219     {
 220         //[thrust_phase_ensemble_define_dopri5
 221         typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type > stepper_type;
 222         //]
 223
 224         ofstream fout( "phase_ensemble_dopri5.dat" );
 225         timer.restart();
 226         for( value_type epsilon = epsilon_min ; epsilon < epsilon_max ; epsilon += d_epsilon )
 227         {
 228             ensemble.set_epsilon( epsilon );
 229             statistics_observer obs;
 230             state_type x = x_host;
 231
 232             timer_local.restart();
 233
 234             // calculate some transients steps
 235             //[ thrust_phase_ensemble_integration
 236             size_t steps1 = integrate_const( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , boost::ref( ensemble ) , x , 0.0 , t_transients , dt );
 237             //]
 238
 239             // integrate and compute the statistics
 240             size_t steps2 = integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , boost::ref( ensemble ) , x , 0.0 , t_max , dt , boost::ref( obs ) );
 241
 242             fout << epsilon << "\t" << obs.get_K_mean() << endl;
 243             cout << "Dopri5 : " << epsilon << "\t" << obs.get_K_mean() << "\t" << timer_local.elapsed() << "\t" << steps1 << "\t" << steps2 << endl;
 244         }
 245         dopri5_time = timer.elapsed();
 246     }
 247
 248
 249
 250     {
 251         //[ thrust_phase_ensemble_define_rk4
 252         typedef runge_kutta4< state_type , value_type , state_type , value_type > stepper_type;
 253         //]
 254
 255         ofstream fout( "phase_ensemble_rk4.dat" );
 256         timer.restart();
 257         for( value_type epsilon = epsilon_min ; epsilon < epsilon_max ; epsilon += d_epsilon )
 258         {
 259             ensemble.set_epsilon( epsilon );
 260             statistics_observer obs;
 261             state_type x = x_host;
 262
 263             timer_local.restart();
 264
 265             // calculate some transients steps
 266             size_t steps1 = integrate_const( stepper_type() , boost::ref( ensemble ) , x , 0.0 , t_transients , dt );
 267
 268             // integrate and compute the statistics
 269             size_t steps2 = integrate_const( stepper_type() , boost::ref( ensemble ) , x , 0.0 , t_max , dt , boost::ref( obs ) );
 270             fout << epsilon << "\t" << obs.get_K_mean() << endl;
 271             cout << "RK4     : " << epsilon << "\t" << obs.get_K_mean() << "\t" << timer_local.elapsed() << "\t" << steps1 << "\t" << steps2 << endl;
 272         }
 273         rk4_time = timer.elapsed();
 274     }
 275
 276     cout << "Dopri 5 : " << dopri5_time << " s\n";
 277     cout << "RK4     : " << rk4_time << "\n";
 278
 279     return 0;
 280 }