refactor all integrate functions

all integrate function now use the adaption to provide checking functionality
the implementations in detail::integrate_* are not aware of the checkers at
all.
This commit is contained in:
Mario Mulansky 2015-10-28 13:49:27 -06:00
parent 5efc181f4c
commit c3856f83a5
8 changed files with 109 additions and 90 deletions

View File

@ -43,11 +43,11 @@ namespace odeint {
namespace detail {
// forward declaration
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , StepOverflowChecker checker , stepper_tag );
Observer observer , stepper_tag );
/*
* integrate_adaptive for simple stepper is basically an integrate_const + some last step
@ -60,7 +60,7 @@ size_t integrate_adaptive(
)
{
size_t steps = detail::integrate_const( stepper , system , start_state , start_time ,
end_time , dt , observer , null_checker() , stepper_tag() );
end_time , dt , observer , stepper_tag() );
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
@ -76,18 +76,17 @@ size_t integrate_adaptive(
/*
* small helper function that calls a given checker - this is for use in integrate_const for ControlledStepper
* integrate adaptive for controlled stepper
*/
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_adaptive_checked(
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , StepOverflowChecker checker
Observer observer , controlled_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
typename odeint::unwrap_reference< StepOverflowChecker >::type &chk = checker;
failed_step_checker fail_checker; // to throw a runtime_error if step size adjustment fails
size_t count = 0;
@ -103,7 +102,6 @@ size_t integrate_adaptive_checked(
do
{
res = st.try_step( system , start_state , start_time , dt );
chk(); // check number of steps
fail_checker(); // check number of failed steps
}
while( res == fail );
@ -116,23 +114,6 @@ size_t integrate_adaptive_checked(
}
/*
* classical integrate adaptive
*/
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , controlled_stepper_tag
)
{
// call the checked version with a null_checker
return integrate_adaptive_checked(stepper, system, start_state, start_time, end_time,
dt, observer, null_checker());
}
/*
* integrate adaptive for dense output steppers
*

View File

@ -31,19 +31,19 @@ namespace odeint {
namespace detail {
// forward declaration
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
size_t integrate_adaptive_checked(
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , StepOverflowChecker checker
Observer observer , controlled_stepper_tag
);
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , StepOverflowChecker /* checker */ , stepper_tag
Observer observer , stepper_tag
)
{
@ -72,15 +72,14 @@ size_t integrate_const(
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , StepOverflowChecker checker , controlled_stepper_tag
Observer observer , controlled_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< StepOverflowChecker >::type &chk = checker;
Time time = start_time;
const Time time_step = dt;
@ -90,11 +89,10 @@ size_t integrate_const(
while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
{
obs( start_state , time );
chk.reset(); // reset after each observation
// integrate_adaptive_checked uses the given checker to throw if an overflow occurs
real_steps += detail::integrate_adaptive_checked(stepper, system, start_state, time,
static_cast<Time>(time + time_step), dt,
null_observer(), checker);
real_steps += detail::integrate_adaptive(stepper, system, start_state, time,
static_cast<Time>(time + time_step), dt,
null_observer(), controlled_stepper_tag());
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
step++;
@ -106,16 +104,15 @@ size_t integrate_const(
}
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_const(
Stepper stepper , System system , State &start_state ,
Time start_time , Time end_time , Time dt ,
Observer observer , StepOverflowChecker checker , dense_output_stepper_tag
Observer observer , dense_output_stepper_tag
)
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
typename odeint::unwrap_reference< StepOverflowChecker >::type &chk = checker;
Time time = start_time;
@ -132,7 +129,6 @@ size_t integrate_const(
{
st.calc_state( time , start_state );
obs( start_state , time );
chk.reset(); // reset checker at observation points
++obs_step;
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
@ -146,7 +142,6 @@ size_t integrate_const(
while( less_eq_with_sign( st.current_time() , time , dt ) )
{
st.do_step( system );
chk();
++real_step;
}
}
@ -154,7 +149,6 @@ size_t integrate_const(
{ // do the last step ending exactly on the end point
st.initialize( st.current_state() , st.current_time() , end_time - st.current_time() );
st.do_step( system );
chk();
++real_step;
}

View File

@ -31,20 +31,20 @@ namespace odeint {
namespace detail {
// forward declaration
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer >
size_t integrate_adaptive_checked(
Stepper stepper , System system , State &start_state ,
Time &start_time , Time end_time , Time &dt ,
Observer observer , StepOverflowChecker checker
Observer observer, controlled_stepper_tag
);
/* basic version */
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer>
Time integrate_n_steps(
Stepper stepper , System system , State &start_state ,
Time start_time , Time dt , size_t num_of_steps ,
Observer observer , StepOverflowChecker /* checker */ , stepper_tag )
Observer observer , stepper_tag )
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
@ -66,14 +66,13 @@ Time integrate_n_steps(
/* controlled version */
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer >
Time integrate_n_steps(
Stepper stepper , System system , State &start_state ,
Time start_time , Time dt , size_t num_of_steps ,
Observer observer , StepOverflowChecker checker , controlled_stepper_tag )
Observer observer , controlled_stepper_tag )
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< StepOverflowChecker >::type &chk = checker;
Time time = start_time;
Time time_step = dt;
@ -81,10 +80,9 @@ Time integrate_n_steps(
for( size_t step = 0; step < num_of_steps ; ++step )
{
obs( start_state , time );
chk.reset(); // reset after each observation
// integrate_adaptive_checked uses the given checker to throw if an overflow occurs
detail::integrate_adaptive_checked(stepper, system, start_state, time, static_cast<Time>(time + time_step), dt,
null_observer(), checker);
detail::integrate_adaptive(stepper, system, start_state, time, static_cast<Time>(time + time_step), dt,
null_observer(), controlled_stepper_tag());
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
time = start_time + static_cast< typename unit_value_type<Time>::type >(step+1) * time_step;
@ -96,15 +94,14 @@ Time integrate_n_steps(
/* dense output version */
template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
template< class Stepper , class System , class State , class Time , class Observer >
Time integrate_n_steps(
Stepper stepper , System system , State &start_state ,
Time start_time , Time dt , size_t num_of_steps ,
Observer observer , StepOverflowChecker checker , dense_output_stepper_tag )
Observer observer , dense_output_stepper_tag )
{
typename odeint::unwrap_reference< Observer >::type &obs = observer;
typename odeint::unwrap_reference< Stepper >::type &st = stepper;
typename odeint::unwrap_reference< StepOverflowChecker >::type &chk = checker;
Time time = start_time;
const Time end_time = start_time + static_cast< typename unit_value_type<Time>::type >(num_of_steps) * dt;
@ -119,7 +116,6 @@ Time integrate_n_steps(
{
st.calc_state( time , start_state );
obs( start_state , time );
chk.reset(); // reset checker after every observer call
++step;
// direct computation of the time avoids error propagation happening when using time += dt
// we need clumsy type analysis to get boost units working here
@ -132,13 +128,11 @@ Time integrate_n_steps(
st.current_time_step() ) )
{
st.do_step( system );
chk(); // potential check for too much work exception
}
else if( less_with_sign( st.current_time() , end_time , st.current_time_step() ) )
{ // do the last step ending exactly on the end point
st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
st.do_step( system );
chk();
}
}
@ -150,7 +144,6 @@ Time integrate_n_steps(
st.current_time_step() ) )
st.initialize( st.current_state() , st.current_time() , static_cast<Time>(end_time - st.current_time()) );
st.do_step( system );
chk();
}
// observation at final point

View File

@ -27,7 +27,6 @@
#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
namespace boost {

View File

@ -24,6 +24,7 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/integrate/null_observer.hpp>
#include <boost/numeric/odeint/integrate/null_checker.hpp>
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
@ -51,9 +52,16 @@ size_t integrate_const(
observer, stepper_category());
}
else {
return detail::integrate_const(stepper, system, start_state,
// unwrap references
typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
typedef typename odeint::unwrap_reference< Observer >::type observer_type;
typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
return detail::integrate_const(checked_stepper<stepper_type, checker_type>(stepper, checker),
system, start_state,
start_time, end_time, dt,
observer, checker, stepper_category());
checked_observer<observer_type, checker_type>(observer, checker),
stepper_category());
}
}
@ -69,7 +77,7 @@ size_t integrate_const(
) {
typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
// we want to get as fast as possible to the end
// ToDo: forward checker to this call
if (boost::is_same<null_observer, Observer>::value) {
return detail::integrate_adaptive(
stepper, system, start_state,
@ -77,9 +85,15 @@ size_t integrate_const(
observer, stepper_category());
}
else {
return detail::integrate_const(stepper, system, start_state,
typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
typedef typename odeint::unwrap_reference< Observer >::type observer_type;
typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
return detail::integrate_const(checked_stepper<stepper_type, checker_type>(stepper, checker),
system, start_state,
start_time, end_time, dt,
observer, checker, stepper_category());
checked_observer<observer_type, checker_type>(observer, checker),
stepper_category());
}
}
@ -90,10 +104,11 @@ size_t integrate_const(
template<class Stepper, class System, class State, class Time, class Observer>
size_t integrate_const(
Stepper stepper, System system, State &start_state,
Time start_time, Time end_time, Time dt, Observer observer
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, observer,
null_checker());
Time start_time, Time end_time, Time dt, Observer observer)
{
typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
return detail::integrate_const(stepper, system, start_state,
start_time, end_time, dt, observer, stepper_category());
}
/**
@ -105,8 +120,9 @@ size_t integrate_const(
Stepper stepper, System system, const State &start_state,
Time start_time, Time end_time, Time dt, Observer observer
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, observer,
null_checker());
typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
return detail::integrate_const(stepper, system, start_state,
start_time, end_time, dt, observer, stepper_category());
}
@ -118,8 +134,7 @@ size_t integrate_const(
Stepper stepper, System system, State &start_state,
Time start_time, Time end_time, Time dt
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer(),
null_checker());
return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer());
}
/**
@ -131,8 +146,7 @@ size_t integrate_const(
Stepper stepper, System system, const State &start_state,
Time start_time, Time end_time, Time dt
) {
return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer(),
null_checker());
return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer());
}

View File

@ -41,11 +41,18 @@ Time integrate_n_steps(
Time start_time , Time dt , size_t num_of_steps ,
Observer observer , StepOverflowChecker checker )
{
typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
// unwrap references
typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
typedef typename odeint::unwrap_reference< Observer >::type observer_type;
typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
typedef typename stepper_type::stepper_category stepper_category;
return detail::integrate_n_steps(
stepper , system , start_state ,
checked_stepper<stepper_type, checker_type>(stepper, checker),
system , start_state ,
start_time , dt , num_of_steps ,
observer , checker , stepper_category() );
checked_observer<observer_type, checker_type>(observer, checker),
stepper_category() );
}
/**
@ -57,11 +64,17 @@ Time integrate_n_steps(
Time start_time , Time dt , size_t num_of_steps ,
Observer observer , StepOverflowChecker checker )
{
typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
typedef typename odeint::unwrap_reference< Observer >::type observer_type;
typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
typedef typename stepper_type::stepper_category stepper_category;
return detail::integrate_n_steps(
stepper , system , start_state ,
start_time , dt , num_of_steps ,
observer , checker , stepper_category() );
checked_stepper<stepper_type, checker_type>(stepper, checker),
system , start_state ,
start_time , dt , num_of_steps ,
checked_observer<observer_type, checker_type>(observer, checker),
stepper_category() );
}
@ -73,7 +86,12 @@ Time integrate_n_steps(
Stepper stepper , System system , State &start_state ,
Time start_time , Time dt , size_t num_of_steps , Observer observer )
{
return integrate_n_steps(stepper, system, start_state, start_time, dt, num_of_steps, observer, null_checker());
typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
return detail::integrate_n_steps(
stepper , system , start_state ,
start_time , dt , num_of_steps ,
observer , stepper_category() );
}
/**
@ -84,7 +102,12 @@ Time integrate_n_steps(
Stepper stepper , System system , const State &start_state ,
Time start_time , Time dt , size_t num_of_steps , Observer observer )
{
return integrate_n_steps(stepper, system, start_state, start_time, dt, num_of_steps, observer, null_checker());
typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
return detail::integrate_n_steps(
stepper , system , start_state ,
start_time , dt , num_of_steps ,
observer , stepper_category() );
}
/**
@ -95,8 +118,8 @@ Time integrate_n_steps(
Stepper stepper , System system , State &start_state ,
Time start_time , Time dt , size_t num_of_steps )
{
return integrate_n_steps(stepper, system, start_state, start_time, dt, num_of_steps, null_observer(),
null_checker());
return integrate_n_steps(stepper, system, start_state, start_time,
dt, num_of_steps, null_observer());
}
/**
@ -107,8 +130,8 @@ Time integrate_n_steps(
Stepper stepper , System system , const State &start_state ,
Time start_time , Time dt , size_t num_of_steps )
{
return integrate_n_steps(stepper, system, start_state, start_time, dt, num_of_steps, null_observer(),
null_checker());
return integrate_n_steps(stepper, system, start_state, start_time,
dt, num_of_steps, null_observer());
}

View File

@ -25,6 +25,7 @@
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
#include <boost/numeric/odeint/integrate/null_observer.hpp>
#include <boost/numeric/odeint/integrate/null_checker.hpp>
#include <boost/numeric/odeint/integrate/check_adapter.hpp>
#include <boost/numeric/odeint/integrate/detail/integrate_times.hpp>
namespace boost {

View File

@ -65,6 +65,7 @@ struct push_back_time
};
typedef runge_kutta_dopri5<state_type> stepper_type;
typedef runge_kutta_cash_karp54<state_type> cash_karp_type;
BOOST_AUTO_TEST_SUITE( integrate_overflow )
@ -173,6 +174,19 @@ BOOST_AUTO_TEST_CASE( test_integrate_times )
t0, t1, 1.0 , push_back_time(times), max_step_checker(10)),
std::runtime_error);
// check cash karp controlled stepper
// no exception if no checker is provided
integrate_times(make_controlled<cash_karp_type>(1E-5, 1E-5), lorenz, x, t0, t1, 1.0 , push_back_time(times));
// very small error terms -> standard overflow threshold of 500 should fire an exception
BOOST_CHECK_THROW(integrate_times(make_controlled<cash_karp_type>(1E-15, 1E-15), lorenz, x,
t0, t1, 1.0 , push_back_time(times), max_step_checker()),
std::runtime_error);
// small threshold of 10 -> larger error still gives an exception
BOOST_CHECK_THROW(integrate_times(make_controlled<cash_karp_type>(1E-5, 1E-5), lorenz, x,
t0, t1, 1.0 , push_back_time(times), max_step_checker(10)),
std::runtime_error);
// check exceptions for dense output stepper
integrate_times(make_dense_output<stepper_type>(1E-5, 1E-5), lorenz, x, t0, t1, 1.0 ,
push_back_time(times));