boost/libs/numeric/odeint/examples/adaptive_iterator.cpp

/*
 * adaptive_iterator.cpp
 *
 * Copyright 2012-2013 Karsten Ahnert
 * Copyright 2012 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)
 */



#include <iostream>
#include <iterator>
#include <utility>
#include <algorithm>
#include <cassert>

#include <boost/array.hpp>

#include <boost/range/algorithm.hpp>
#include <boost/range/adaptor/filtered.hpp>
#include <boost/range/numeric.hpp>

#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
#include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
#include <boost/numeric/odeint/stepper/generation.hpp>

#include <boost/numeric/odeint/iterator/adaptive_iterator.hpp>
#include <boost/numeric/odeint/iterator/adaptive_time_iterator.hpp>

#define tab "\t"

using namespace std;
using namespace boost::numeric::odeint;

const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;

struct lorenz
{
    template< class State , class Deriv >
    void operator()( const State &x , Deriv &dxdt , double t ) const
    {
        dxdt[0] = sigma * ( x[1] - x[0] );
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = -b * x[2] + x[0] * x[1];
    }
};

#include <typeinfo>

int main( int argc , char **argv )
{
    typedef boost::array< double , 3 > state_type;

    /*
     * Controlled steppers with time iterator
     */

    // std::for_each
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        std::for_each( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                       make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
                       []( const std::pair< const state_type&, double > &x ) {
                           std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
    }

    // std::copy_if
    {
        std::vector< pair< state_type , double > > res;
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        std::copy_if( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                      make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
                      std::back_inserter( res ) ,
                      []( const pair< const state_type& , double > &x ) {
                          return ( x.first[0] > 0.0 ) ? true : false; } );
        for( size_t i=0 ; i<res.size() ; ++i )
            cout << res[i].first[0] << tab << res[i].first[1] << tab << res[i].first[2] << "\n";
    }

    // std::accumulate
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        double res = std::accumulate( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                      make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
                                      0.0 ,
                                      []( double sum , const pair< const state_type& , double > &x ) {
                                          return sum + x.first[0]; } );
        cout << res << endl;
    }


    // std::transform
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        vector< double > weights;
        std::transform( make_adaptive_time_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                        make_adaptive_time_iterator_end( stepper , lorenz() , x ) ,
                        back_inserter( weights ) ,
                        []( const pair< const state_type& , double > &x ) {
                            return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
        for( size_t i=0 ; i<weights.size() ; ++i )
            cout << weights[i] << "\n";
    }













    /*
     * Boost.Range versions of controlled stepper with time iterator
     */


    // boost::range::for_each
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        boost::range::for_each( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                []( const std::pair< const state_type& , double > &x ) {
                                    std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );
    }


    // boost::range::copy with filtered adaptor (simulating std::copy_if)
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        std::vector< std::pair< state_type , double > > res;
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        boost::range::copy( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
                            boost::adaptors::filtered( [] ( const pair< const state_type& , double > &x ) { return ( x.first[0] > 0.0 ); } ) ,
                            std::back_inserter( res ) );
        for( size_t i=0 ; i<res.size() ; ++i )
            cout << res[i].first[0] << tab << res[i].first[1] << tab << res[i].first[2] << "\n";
    }

    // boost::range::accumulate
    {
        //[adaptive_time_iterator_accumulate_range
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        double res = boost::accumulate( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
                                        []( double sum , const pair< const state_type& , double > &x ) {
                                            return sum + x.first[0]; } );
        cout << res << endl;
        //]
    }


    //  boost::range::transform
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        vector< double > weights;
        boost::transform( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
                          []( const pair< const state_type& , double > &x ) {
                              return sqrt( x.first[0] * x.first[0] + x.first[1] * x.first[1] + x.first[2] * x.first[2] ); } );
        for( size_t i=0 ; i<weights.size() ; ++i )
            cout << weights[i] << "\n";
    }


    // boost::range::find with time iterator
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        auto iter = boost::find_if( make_adaptive_time_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                    []( const std::pair< const state_type & , double > &x ) {
                                        return ( x.first[0] < 0.0 ); } );
        cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";
    }












    // /*
    //  * Boost.Range versions for dense output steppers
    //  */

    // // boost::range::for_each
    // {
    //     runge_kutta_dopri5< state_type > stepper;
    //     state_type x = {{ 10.0 , 10.0 , 10.0 }};
    //     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 ) ,
    //                             []( const state_type &x ) {
    //                                 std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
    // }

    
    // // boost::range::for_each with time iterator
    // {
    //     runge_kutta_dopri5< state_type > stepper;
    //     state_type x = {{ 10.0 , 10.0 , 10.0 }};
    //     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 ) ,
    //                             []( const std::pair< state_type& , double > &x ) {
    //                                 std::cout << x.second << tab << x.first[0] << tab << x.first[1] << tab << x.first[2] << "\n"; } );

    // }





    /*
     * Pure iterators for controlled stepper without time iterator
     */

    // std::for_each
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        std::for_each( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                       make_adaptive_iterator_end( stepper , lorenz() , x ) ,
                       []( const state_type& x ) {
                           std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
    }

    // std::copy_if
    {
        std::vector< state_type > res;
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        std::copy_if( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                      make_adaptive_iterator_end( stepper , lorenz() , x ) ,
                      std::back_inserter( res ) ,
                      []( const state_type& x ) {
                          return ( x[0] > 0.0 ) ? true : false; } );
        for( size_t i=0 ; i<res.size() ; ++i )
            cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
    }

    // std::accumulate
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        double res = std::accumulate( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                      make_adaptive_iterator_end( stepper , lorenz() , x ) ,
                                      0.0 ,
                                      []( double sum , const state_type& x ) {
                                          return sum + x[0]; } );
        cout << res << endl;
    }


    // std::transform
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        vector< double > weights;
        std::transform( make_adaptive_iterator_begin( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                        make_adaptive_iterator_end( stepper , lorenz() , x ) ,
                        back_inserter( weights ) ,
                        []( const state_type& x ) {
                            return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
        for( size_t i=0 ; i<weights.size() ; ++i )
            cout << weights[i] << "\n";
    }










    /*
     * Boost.Range versions of controlled stepper WITHOUT time iterator
     */


    // boost::range::for_each
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        boost::range::for_each( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                []( const state_type &x ) {
                                    std::cout << x[0] << tab << x[1] << tab << x[2] << "\n"; } );
    }


    // boost::range::copy with filtered adaptor (simulating std::copy_if)
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        std::vector< state_type > res;
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        boost::range::copy( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) |
                            boost::adaptors::filtered( [] ( const state_type& x ) { return ( x[0] > 0.0 ); } ) ,
                            std::back_inserter( res ) );
        for( size_t i=0 ; i<res.size() ; ++i )
            cout << res[i][0] << tab << res[i][1] << tab << res[i][2] << "\n";
    }

    // boost::range::accumulate
    {
        //[adaptive_iterator_accumulate_range
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        double res = boost::accumulate( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , 0.0 ,
                                        []( double sum , const state_type& x ) {
                                            return sum + x[0]; } );
        cout << res << endl;
        //]
    }


    //  boost::range::transform
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        vector< double > weights;
        boost::transform( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) , back_inserter( weights ) ,
                          []( const state_type& x ) {
                              return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); } );
        for( size_t i=0 ; i<weights.size() ; ++i )
            cout << weights[i] << "\n";
    }


    // boost::range::find 
    {
        auto stepper = make_controlled( 1.0e-6 , 1.0e-6 , runge_kutta_cash_karp54< state_type >() );
        state_type x = {{ 10.0 , 10.0 , 10.0 }};
        auto iter = boost::find_if( make_adaptive_range( stepper , lorenz() , x , 0.0 , 1.0 , 0.01 ) ,
                                    []( const state_type &x ) {
                                        return ( x[0] < 0.0 ); } );
        cout << (*iter)[0] << "\t" << (*iter)[1] << "\t" << (*iter)[2] << "\n";
    }





    return 0;
}