1 /*
2 * adaptive_iterator.cpp
3 *
4 * Copyright 2012-2013 Karsten Ahnert
5 * Copyright 2012 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
12
13
14 #include <iostream>
15 #include <iterator>
16 #include <utility>
17 #include <algorithm>
18 #include <cassert>
19
20 #include <boost/array.hpp>
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/runge_kutta_cash_karp54.hpp>
29 #include <boost/numeric/odeint/stepper/generation.hpp>
30
31 #include <boost/numeric/odeint/iterator/adaptive_iterator.hpp>
32 #include <boost/numeric/odeint/iterator/adaptive_time_iterator.hpp>
33
34 #define tab "\t"
35
36 using namespace std;
37 using namespace boost::numeric::odeint;
38
39 const double sigma = 10.0;
40 const double R = 28.0;
41 const double b = 8.0 / 3.0;
42
43 struct lorenz
44 {
45 template< class State , class Deriv >
operator ()lorenz46 void operator()( const State &x , Deriv &dxdt , double t ) const
47 {
48 dxdt[0] = sigma * ( x[1] - x[0] );
49 dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
50 dxdt[2] = -b * x[2] + x[0] * x[1];
51 }
52 };
53
54 #include <typeinfo>
55
main(int argc,char ** argv)56 int main( int argc , char **argv )
57 {
58 typedef boost::array< double , 3 > state_type;
59
60 /*
61 * Controlled steppers with time iterator
62 */
63
64 // std::for_each
65 {
66 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
67 state_type x = {{ 10.0 , 10.0 , 10.0 }};
68 std::for_each( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
69 make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
70 []( const std::pair< const state_type&, double > &x ) {
71 std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
72 }
73
74 // std::copy_if
75 {
76 std::vector< pair< state_type , double > > res;
77 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
78 state_type x = {{ 10.0 , 10.0 , 10.0 }};
79 std::copy_if( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
80 make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
81 std::back_inserter( res ) ,
82 []( const pair< const state_type& , double > &x ) {
83 return ( x.first[0] > 0.0 ) ? true : false; } );
84 for( size_t i=0 ; i<res.size() ; ++i )
85 cout << res[i].first[0] << tab << res[i].first[1] << tab << res[i].first[2] << "\n";
86 }
87
88 // std::accumulate
89 {
90 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
91 state_type x = {{ 10.0 , 10.0 , 10.0 }};
92 double res = std::accumulate( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
93 make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
94 0.0 ,
95 []( double sum , const pair< const state_type& , double > &x ) {
96 return sum + x.first[0]; } );
97 cout << res << endl;
98 }
99
100
101 // std::transform
102 {
103 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
104 state_type x = {{ 10.0 , 10.0 , 10.0 }};
105 vector< double > weights;
106 std::transform( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
107 make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
108 back_inserter( weights ) ,
109 []( const pair< const state_type& , double > &x ) {
110 return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
111 for( size_t i=0 ; i<weights.size() ; ++i )
112 cout << weights[i] << "\n";
113 }
114
115
116
117
118
119
120
121
122
123
124
125
126
127 /*
128 * Boost.Range versions of controlled stepper with time iterator
129 */
130
131
132 // boost::range::for_each
133 {
134 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
135 state_type x = {{ 10.0 , 10.0 , 10.0 }};
136 boost::range::for_each( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
137 []( const std::pair< const state_type& , double > &x ) {
138 std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
139 }
140
141
142 // boost::range::copy with filtered adaptor (simulating std::copy_if)
143 {
144 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
145 std::vector< std::pair< state_type , double > > res;
146 state_type x = {{ 10.0 , 10.0 , 10.0 }};
147 boost::range::copy( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
148 boost::adaptors::filtered( [] ( const pair< const state_type& , double > &x ) { return ( x.first[0] > 0.0 ); } ) ,
149 std::back_inserter( res ) );
150 for( size_t i=0 ; i<res.size() ; ++i )
151 cout << res[i].first[0] << tab << res[i].first[1] << tab << res[i].first[2] << "\n";
152 }
153
154 // boost::range::accumulate
155 {
156 //[adaptive_time_iterator_accumulate_range
157 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
158 state_type x = {{ 10.0 , 10.0 , 10.0 }};
159 double res = boost::accumulate( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
160 []( double sum , const pair< const state_type& , double > &x ) {
161 return sum + x.first[0]; } );
162 cout << res << endl;
163 //]
164 }
165
166
167 // boost::range::transform
168 {
169 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
170 state_type x = {{ 10.0 , 10.0 , 10.0 }};
171 vector< double > weights;
172 boost::transform( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
173 []( const pair< const state_type& , double > &x ) {
174 return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
175 for( size_t i=0 ; i<weights.size() ; ++i )
176 cout << weights[i] << "\n";
177 }
178
179
180 // boost::range::find with time iterator
181 {
182 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
183 state_type x = {{ 10.0 , 10.0 , 10.0 }};
184 auto iter = boost::find_if( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
185 []( const std::pair< const state_type & , double > &x ) {
186 return ( x.first[0] < 0.0 ); } );
187 cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";
188 }
189
190
191
192
193
194
195
196
197
198
199
200
201 // /*
202 // * Boost.Range versions for dense output steppers
203 // */
204
205 // // boost::range::for_each
206 // {
207 // runge_kutta_dopri5< state_type > stepper;
208 // state_type x = {{ 10.0 , 10.0 , 10.0 }};
209 // boost::range::for_each( make_adaptive_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
210 // []( const state_type &x ) {
211 // std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
212 // }
213
214
215 // // boost::range::for_each with time iterator
216 // {
217 // runge_kutta_dopri5< state_type > stepper;
218 // state_type x = {{ 10.0 , 10.0 , 10.0 }};
219 // boost::range::for_each( make_adaptive_time_range( make_dense_output( 1.0e-6 , 1.0e-6 , stepper ) , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
220 // []( const std::pair< state_type& , double > &x ) {
221 // std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
222
223 // }
224
225
226
227
228
229 /*
230 * Pure iterators for controlled stepper without time iterator
231 */
232
233 // std::for_each
234 {
235 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
236 state_type x = {{ 10.0 , 10.0 , 10.0 }};
237 std::for_each( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
238 make_adaptive_iterator_end( stepper , lorenz() , x ) ,
239 []( const state_type& x ) {
240 std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
241 }
242
243 // std::copy_if
244 {
245 std::vector< state_type > res;
246 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
247 state_type x = {{ 10.0 , 10.0 , 10.0 }};
248 std::copy_if( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
249 make_adaptive_iterator_end( stepper , lorenz() , x ) ,
250 std::back_inserter( res ) ,
251 []( const state_type& x ) {
252 return ( x[0] > 0.0 ) ? true : false; } );
253 for( size_t i=0 ; i<res.size() ; ++i )
254 cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
255 }
256
257 // std::accumulate
258 {
259 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
260 state_type x = {{ 10.0 , 10.0 , 10.0 }};
261 double res = std::accumulate( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
262 make_adaptive_iterator_end( stepper , lorenz() , x ) ,
263 0.0 ,
264 []( double sum , const state_type& x ) {
265 return sum + x[0]; } );
266 cout << res << endl;
267 }
268
269
270 // std::transform
271 {
272 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
273 state_type x = {{ 10.0 , 10.0 , 10.0 }};
274 vector< double > weights;
275 std::transform( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
276 make_adaptive_iterator_end( stepper , lorenz() , x ) ,
277 back_inserter( weights ) ,
278 []( const state_type& x ) {
279 return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
280 for( size_t i=0 ; i<weights.size() ; ++i )
281 cout << weights[i] << "\n";
282 }
283
284
285
286
287
288
289
290
291
292
293 /*
294 * Boost.Range versions of controlled stepper WITHOUT time iterator
295 */
296
297
298 // boost::range::for_each
299 {
300 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
301 state_type x = {{ 10.0 , 10.0 , 10.0 }};
302 boost::range::for_each( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
303 []( const state_type &x ) {
304 std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
305 }
306
307
308 // boost::range::copy with filtered adaptor (simulating std::copy_if)
309 {
310 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
311 std::vector< state_type > res;
312 state_type x = {{ 10.0 , 10.0 , 10.0 }};
313 boost::range::copy( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
314 boost::adaptors::filtered( [] ( const state_type& x ) { return ( x[0] > 0.0 ); } ) ,
315 std::back_inserter( res ) );
316 for( size_t i=0 ; i<res.size() ; ++i )
317 cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
318 }
319
320 // boost::range::accumulate
321 {
322 //[adaptive_iterator_accumulate_range
323 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
324 state_type x = {{ 10.0 , 10.0 , 10.0 }};
325 double res = boost::accumulate( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
326 []( double sum , const state_type& x ) {
327 return sum + x[0]; } );
328 cout << res << endl;
329 //]
330 }
331
332
333 // boost::range::transform
334 {
335 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
336 state_type x = {{ 10.0 , 10.0 , 10.0 }};
337 vector< double > weights;
338 boost::transform( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
339 []( const state_type& x ) {
340 return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
341 for( size_t i=0 ; i<weights.size() ; ++i )
342 cout << weights[i] << "\n";
343 }
344
345
346 // boost::range::find
347 {
348 auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
349 state_type x = {{ 10.0 , 10.0 , 10.0 }};
350 auto iter = boost::find_if( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
351 []( const state_type &x ) {
352 return ( x[0] < 0.0 ); } );
353 cout << (*iter)[0] << "\t" << (*iter)[1] << "\t" << (*iter)[2] << "\n";
354 }
355
356
357
358
359
360 return 0;
361 }
362