1 /*
2  * const_step_iterator.cpp
3  *
4  * Copyright 2012-2013 Karsten Ahnert
5  * Copyright 2013 Mario Mulansky
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  * several examples for using iterators
12  */
13 
14 
15 #include <iostream>
16 #include <iterator>
17 #include <utility>
18 #include <algorithm>
19 #include <array>
20 #include <cassert>
21 
22 #include <boost/range/algorithm.hpp>
23 #include <boost/range/adaptor/filtered.hpp>
24 #include <boost/range/numeric.hpp>
25 
26 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
27 #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
28 #include <boost/numeric/odeint/stepper/generation.hpp>
29 #include <boost/numeric/odeint/iterator/const_step_iterator.hpp>
30 #include <boost/numeric/odeint/iterator/const_step_time_iterator.hpp>
31 
32 #define tab "\t"
33 
34 using namespace std;
35 using namespace boost::numeric::odeint;
36 
37 const double sigma = 10.0;
38 const double R = 28.0;
39 const double b = 8.0 / 3.0;
40 
41 struct lorenz
42 {
43     template< class State , class Deriv >
operator ()lorenz44     void operator()( const State &x , Deriv &dxdt , double t ) const
45     {
46         dxdt[0] = sigma * ( x[1] - x[0] );
47         dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
48         dxdt[2] = -b * x[2] + x[0] * x[1];
49     }
50 };
51 
52 
53 
main(int argc,char ** argv)54 int main( int argc , char **argv )
55 {
56     typedef std::array< double , 3 > state_type;
57 
58     // std::for_each
59     {
60         runge_kutta4< state_type > stepper;
61         state_type x = {{ 10.0 , 10.0 , 10.0 }};
62         std::for_each( make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
63                        make_const_step_time_iterator_end( stepper , lorenz() , x ) ,
64                        []( const std::pair< const state_type&, double > &x ) {
65                            std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
66     }
67 
68     // std::copy_if
69     {
70         std::vector< state_type > res;
71         runge_kutta4< state_type > stepper;
72         state_type x = {{ 10.0 , 10.0 , 10.0 }};
73         std::copy_if( make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
74                       make_const_step_iterator_end( stepper , lorenz() , x ) ,
75                       std::back_inserter( res ) ,
76                       []( const state_type& x ) {
77                           return ( x[0] > 0.0 ) ? true : false; } );
78         for( size_t i=0 ; i<res.size() ; ++i )
79             cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
80     }
81 
82     // std::accumulate
83     {
84         //[ const_step_iterator_accumulate
85         runge_kutta4< state_type > stepper;
86         state_type x = {{ 10.0 , 10.0 , 10.0 }};
87         double res = std::accumulate( make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
88                                       make_const_step_iterator_end( stepper , lorenz() , x ) ,
89                                       0.0 ,
90                                       []( double sum , const state_type &x ) {
91                                           return sum + x[0]; } );
92         cout << res << endl;
93         //]
94     }
95 
96 
97     // std::transform
98     {
99         runge_kutta4< state_type > stepper;
100         state_type x = {{ 10.0 , 10.0 , 10.0 }};
101         vector< double > weights;
102         std::transform( make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
103                         make_const_step_iterator_end( stepper , lorenz() , x ) ,
104                         back_inserter( weights ) ,
105                         []( const state_type &x ) {
106                             return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
107         for( size_t i=0 ; i<weights.size() ; ++i )
108             cout << weights[i] << "\n";
109     }
110 
111 
112 
113     // std::transform with time_iterator
114     {
115         runge_kutta4< state_type > stepper;
116         state_type x = {{ 10.0 , 10.0 , 10.0 }};
117         vector< double > weights;
118         std::transform( make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
119                         make_const_step_time_iterator_end( stepper , lorenz() , x ) ,
120                         back_inserter( weights ) ,
121                         []( const std::pair< const state_type &, double > &x ) {
122                             return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
123         for( size_t i=0 ; i<weights.size() ; ++i )
124             cout << weights[i] << "\n";
125     }
126 
127 
128 
129 
130 
131 
132 
133 
134 
135 
136     // /*
137     //  * Boost.Range versions
138     //  */
139 
140 
141     // boost::range::for_each
142     {
143         runge_kutta4< state_type > stepper;
144         state_type x = {{ 10.0 , 10.0 , 10.0 }};
145         boost::range::for_each( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
146                                 []( const state_type &x ) {
147                                     std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
148     }
149 
150     // boost::range::for_each with time iterator
151     {
152         runge_kutta4< state_type > stepper;
153         state_type x = {{ 10.0 , 10.0 , 10.0 }};
154         boost::range::for_each( make_const_step_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
155                                 []( const std::pair< const state_type& , double > &x ) {
156                                     std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
157     }
158 
159 
160     // boost::range::copy with filtered adaptor (simulating std::copy_if)
161     {
162         std::vector< state_type > res;
163         runge_kutta4< state_type > stepper;
164         state_type x = {{ 10.0 , 10.0 , 10.0 }};
165         boost::range::copy( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
166                             boost::adaptors::filtered( [] ( const state_type &x ) { return ( x[0] > 0.0 ); } ) ,
167                             std::back_inserter( res ) );
168         for( size_t i=0 ; i<res.size() ; ++i )
169             cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
170     }
171 
172     // boost::range::accumulate
173     {
174         //[const_step_iterator_accumulate_range
175         runge_kutta4< state_type > stepper;
176         state_type x = {{ 10.0 , 10.0 , 10.0 }};
177         double res = boost::accumulate( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
178                                         []( double sum , const state_type &x ) {
179                                             return sum + x[0]; } );
180         cout << res << endl;
181         //]
182     }
183 
184     // boost::range::accumulate with time iterator
185     {
186         //[const_step_time_iterator_accumulate_range
187         runge_kutta4< state_type > stepper;
188         state_type x = {{ 10.0 , 10.0 , 10.0 }};
189         double res = boost::accumulate( make_const_step_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
190                                         []( double sum , const std::pair< const state_type &, double > &x ) {
191                                             return sum + x.first[0]; } );
192         cout << res << endl;
193         //]
194     }
195 
196 
197     //  boost::range::transform
198     {
199         runge_kutta4< state_type > stepper;
200         state_type x = {{ 10.0 , 10.0 , 10.0 }};
201         vector< double > weights;
202         boost::transform( make_const_step_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
203                           []( const state_type &x ) {
204                               return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
205         for( size_t i=0 ; i<weights.size() ; ++i )
206             cout << weights[i] << "\n";
207     }
208 
209 
210     // boost::range::find with time iterator
211     {
212         runge_kutta4< state_type > stepper;
213         state_type x = {{ 10.0 , 10.0 , 10.0 }};
214         auto iter = boost::find_if( make_const_step_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
215                                     []( const std::pair< const state_type & , double > &x ) {
216                                         return ( x.first[0] < 0.0 ); } );
217         cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";
218 
219     }
220 
221 
222 
223 
224 
225 
226 
227 
228 
229 
230 
231 
232     /*
233      * Boost.Range versions for dense output steppers
234      */
235 
236     // boost::range::for_each
237     {
238         runge_kutta_dopri5< state_type > stepper;
239         state_type x = {{ 10.0 , 10.0 , 10.0 }};
240         boost::range::for_each( make_const_step_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
241                                 []( const state_type &x ) {
242                                     std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
243     }
244 
245 
246     // boost::range::for_each with time iterator
247     {
248         runge_kutta_dopri5< state_type > stepper;
249         state_type x = {{ 10.0 , 10.0 , 10.0 }};
250         boost::range::for_each( make_const_step_time_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
251                                 []( const std::pair< const state_type& , double > &x ) {
252                                     std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
253 
254     }
255 
256 
257 
258 
259 
260     /*
261      * Pure iterators
262      */
263     {
264         runge_kutta4< state_type > stepper;
265         state_type x = {{ 10.0 , 10.0 , 10.0 }};
266         auto first = make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 );
267         auto last  = make_const_step_iterator_end( stepper , lorenz() , x );
268         while( first != last )
269         {
270             assert( last != first );
271             cout << (*first)[0] << tab << (*first)[1] << tab << (*first)[2] << "\n";
272             ++first;
273         }
274     }
275 
276     {
277         runge_kutta4< state_type > stepper;
278         state_type x = {{ 10.0 , 10.0 , 10.0 }};
279         auto first = make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 );
280         auto last  = make_const_step_time_iterator_end( stepper , lorenz() , x );
281         while( first != last )
282         {
283             assert( last != first );
284             cout << first->second << tab << first->first[0] << tab << first->first[1] << tab << first->first[2] << "\n";
285             ++first;
286         }
287     }
288 
289 
290 
291 
292 
293 
294 
295     return 0;
296 }
297