ODE System exceptions

[ph03]

Hi! I just found your greate lib and I am just figuring out, if I can use it for my project. I'm doing flow visualization of dataset with a finite domain. The odeint lib could be used to integrate stream / path-lines in this domains. However, since the domains are only finite it will always be possible that a streamline integration leaves the domain. This is easily determined in a FlowField class that samples the field. However, I don't know how to communicate this event to the odeint integrators. I think the best way would be to integrate some kind of exception system to the odeint lib, for example for out of bounds sampling, critical point detection or simply for generic errors related to ODE system evaluation. Would it be hard to implement support for this in a generic way or is the library not designed to handle this kind of ODE systems?

Greets Janick

[headmyshoulder]

Hi Janick,

yesterday we discussed a similar problem. I think this kind of application needs a special integrate function which handles the boundaries. This function might be implemented be the use of a general policy class capable of checking the bounds, for example

integrate_const( stepper , system , x , t_start , t_end , dt , boundary_handler );

boundary_handler is then a functor being called after each step, similar to the observers but which might change the state of the ODE. I think this is doable, although it will need some time for us to implement this.

Just two questions: What kind of steppers do you use? And what happens if a particle or streamline leaves the domain?

[ph03]

Hi! Right, a bounday_handler concept sounds like a very clean solution for this problem. I haven't used odeint yet because of this current limitation. In our implementations most of the time a RungeKutta 45 scheme with adaptive stepsize control is more than sufficient (I think we have it from "Numerical Recipies", yours is called Dopri5 if I'm correct). I think the best thing to happen when a streamline reaches a point outside of the boundary is to discard this point outside of the boundary (not even tell it to the observer) and instead redo a step from the previous point until the boundary is hit. I think depending on the used stepper there can be more specialized ways to solve this problem (actually we sample the polynomial curve that is implicitly generated by the integrator to obtain a point on the boundary), but I think a general viable way is to do some kind of binary search that reduces the stepsize of new points such that the boundary is hit "exactly" / in an epsilon distance. It would be also nice if the information of this "boundary hit" event will be retrievable, since the t_end value will not be reached. Particle movement outside the domain is undefined as there is no vector field data given. However, I'm wondering if a boundary_handler concept is too restictive as there are perhaps more general exceptions possible (e.g. the observer notices some event and wants the integration to stop / redo the step with different stepsize because an error criteria is not met). I therefore think a general policy class implemented by the user that offers a high level of control of the integrator (maybe by a hierachy of different exceptions) and that communicates it to the user should be a good and general solution.

Thanks for your effort in this project! I would love to test any solution you come up with in our project, just tell me what to clone :)

Greets J

[headmyshoulder]

Hmm, after reading about your problem I think that other solutions might also fit well for you.

You might implement a step-size controller proxy which uses a controlled_runge_kutta stepper with a dopri5 stepper and implement your boundary checking logic within this proxy. So, the proxy calls your controlled stepper, performs the boundary checks and rejects or accepts the current step size and adjust the step size for the next step.

Another point you might take into account is dense output which interpolates your solution at arbitrary points between two steps. So, if a streamline leaves your domain in a step you can interpolate the value of streamline at the boundary. But you must be careful, since the current step might have been obtained from informations from the exterior of your range.

[gauteh]

Hi,

is there any news on allowing events to control the integration process? I am currently using the adaptive integrator for tracing rays and I would like to both stop the integration at the boundaries of my model and to catch other non-terminal events, like an intercept of an interface. The observer should allow me to trigger on or at least register non-terminal events, but there is no way as far as I can see to indicate that the integration should stop.

Thanks for the great software.

Update: I realized from a discussion on the mailing list[1] that I could probably do something with iterators.

Best regards, Gaute

[1] http://boost.2283326.n4.nabble.com/A-couple-of-odeint-questions-tt4636453.html#none

[headmyshoulder]

Hi,

sorry there are no news about conditional integrators. It is still on our list. Actually I first want to finalize the iterators before implementing the conditional integrate. The Iterators could also solve your problem, at least in a limited fashion. But they are currently redesgined, because they had some serious problems with memory management. I think the iterators will be finished within the next week. I will keep you informed here.

[hazelnusse]

First of all, thanks for the excellent library.

Second, I'm also looking for a simple way to end integration when some function of the states undergoes a zero crossing. The other two ODE integrators that I am familiar with that have this functionality are Matlab (ode45 et al.) and the LSODAR function that is part of ODEPACK. There are lots of uses cases for this functionality and it would be great if this were part of odeint. I know a number of people who use Matlab solely for this events functionality but don't have the time/inclination/skills to write this functionality.

Are there any efforts at the moment to add this functionality?

[headmyshoulder]

On 03.08.2013 22:18, Dale Lukas Peterson wrote:

First of all, thanks for the excellent library.

Second, I'm also looking for a simple way to end integration when some function of the states undergoes a zero crossing. The other two ODE integrators that I am familiar with that have this functionality are Matlab (ode45 et al.) and the LSODAR function that is part of ODEPACK. There are lots of uses cases for this functionality and it would be great if this were part of odeint. I know a number of people who use Matlab solely for this events functionality but don't have the time/inclination/skills to write this functionality.

Are there any efforts at the moment to add this functionality?

We are working on this, but we are slow at the moment. But you can already use iterators to mimick the same behavior [1,2]. They are in the github version of odeint, some docs are already present, but you have to generate them by hand since these documentation pages are neither in the boost version nor on the webpage. I will make them available very soon on our webpage.

[1] https://github.com/headmyshoulder/odeint-v2/blob/master/libs/numeric/odeint/examples/adaptive_iterator.cpp [2] https://github.com/headmyshoulder/odeint-v2/blob/master/libs/numeric/odeint/examples/const_step_iterator.cpp

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-22060906.

[hazelnusse]

Thanks for pointing me to those examples.

Is there an example which illustrates how to find (to some specified accuracy) the time at which the zero crossing occurs? Also, it would be nice to have an example which allows for setting up the function whose zero crossings you want to find and/or terminate integration on, though I guess the example which uses x[0] > 0.0 ? true : false could be expanded to have more sophisticated logic than just checking one of the states.

Being able to specify the direction of the zero crossing of the function of the states is handy. For example allowing g(x, t) (the function whose zeros you want to find during integration) to go from positive to negative and not trigger the event, but when it goes from negative to positive, have it trigger the event.

[hazelnusse]

Also, if there is anything I can do to help, let me know. I'm still getting comfortable with the library, my knowledge of C++ template is fairly low, but if you have things to test out or just want to bounce ideas around, feel free to do so.

[gauteh]

I achieved events/exceptions (in ray tracing) by re-implementing the adaptive integrator with the event/exception logic in between, it is then fairly easy to re-iterate on an event to reach a given accuracy for the event.

[headmyshoulder]

On 04.08.2013 00:12, Dale Lukas Peterson wrote:

Thanks for pointing me to those examples.

Is there an example which illustrates how to find (to some specified accuracy) the time at which the zero crossing occurs? Also, it would be nice to have an example which allows for setting up the function whose zero crossings you want to find and/or terminate integration on, though I guess the example which uses |x[0] > 0.0 ? true : false| could be expanded to have more sophisticated logic than just checking one of the states.

These are the iterator with a time in the name. Their value type is a pair of state (or reference to the state) and the current time point. I think there are also example how to use them.

Being able to specify the direction of the zero crossing of the function of the states is handy. For example allowing |g(x, t)| (the functionoud to contribute a high quality, modern and stable numerical library into boost. whose zeros you want to find during integration) to go from positive to negative and not trigger the event, but when it goes from negative to positive, have it trigger the event.

Do you also need the possibility to iterate to an give time point / event with a given precision? Or is it sufficient for you to stop the integration if an event has occurred? Without getting exactly to the specific point where the event has occurred.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-22062595.

[headmyshoulder]

On 04.08.2013 04:10, Dale Lukas Peterson wrote:

Also, if there is anything I can do to help, let me know. I'm still getting comfortable with the library, my knowledge of C++ template is fairly low, but if you have things to test out or just want to bounce ideas around, feel free to do so.

Ok, this is good to know. I think we will come back to you if we have any ideas how to implement the condition integration.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-22065467.

[headmyshoulder]

On 04.08.2013 12:44, Gaute Hope wrote:

I achieved events/exceptions (in ray tracing) by re-implementing the adaptive integrator with the event/exception logic in between, it is then fairly easy to re-iterate on an event to reach a given accuracy for the event.

Can you give us your code? We have also some ideas how to implement such steppers. Maybe your existing class can be integration into our ideas and into odeint.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-22069771.

[gauteh]

I will factor it out, but it is not very generic.

[gauteh]

Check out: https://github.com/gauteh/odeint_integrate_event

When an event is detected, e.g.: https://github.com/gauteh/odeint_integrate_event/blob/master/ray.cpp#L260 I passed the stepper and state to a function which would then continue to step it (depending on the problem) to get the desired accuracy for the event. This is because, in my case, I could often quickly determine whether an event had happened between two steps with the overall stepper accuracy, but needed more time to determine exactly when (and at which state) the event happened in the range between the two steps.

When reading, bear in mind that this code has been relatively quickly cut out from a larger code base and some comments/stuff doesn't necessarily seem logical here. The most interesting part is Ray::solve () in ray.cpp.

[dyspropos]

I also made an implementation which required the use of a stepper with dense output. I then used bisection to find the zero-crossing. Here is an example based on the lorenz-example from the docs using to stop the first time the trajectory crosses x = 0: https://github.com/sgammelmark/odeint-v2/blob/6ecbb29271d4819c5246a92a3e8fbe8a5e90fb1f/libs/numeric/odeint/examples/lorenz_event.cpp

Søren

[headmyshoulder]

Hi guys,

Hehe, I have already implemented a first version of integrate_conditional some month ago. It is already quite advanced but not completed yet. Have a look into the integrate_conditional branch. It provides an integrate_conditional function which takes an additional controller argument. The controller will take care about stepping, stopping condition, and what will happen if the condition is fulfilled.

Have a look at the https://github.com/headmyshoulder/odeint-v2/blob/integrate_conditional/libs/numeric/odeint/examples/integrate_conditional.cpp. There are two controllers: conditional_stop which simply stops if a give condition has been fulfilled. The second is approximate which halves the stepsize if a given condition is met until a specified precision is met.

I think all your use case can easily be mapped into specific controllers, especially the case with the bisection. Can you check if this model would work for you? I would then start to finish the implementation of conditional_integrate, create some often used controllers, and write documentation. Of course, any help is appreciated.

Best regards,

Karsten

[arjungm]

Hi Karsten,

Can this handle multiple events in the integrate_conditional, like if you're detecting zero crossings in two state variables?

Thanks for the great functionality on this, it makes switching from MATLAB way easier!

Best, Arjun

[headmyshoulder]

On 01/03/2014 11:44 PM, Arjun Menon wrote:

Hi Karsten,

Can this handle multiple events in the integrate_conditional, like if you're detecting zero crossings in two state variables?

This is definitely possible but you have to write the condition for yourself.

Btw. how do you use integrate_conditional right now? We decided to not go for the integrate_conditional solution to solve the problem of stopping the integration at custom points. It has some real problems which can not be solved generically. Instead we decided to go for the iterator solution an which already works in the master branch. Can you sent us an short example how you use integrate_conditional? Then we can see if and how the iterator approach solves your use case.

Thank you,

Karsten

Thanks for the great functionality on this, it makes switching from MATLAB way easier!

Best, Arjun

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-31560204.

[arjungm]

Can you send us an short example how you use integrate_conditional?

I actually didn't attempt that yet. I just started porting most of my simulations, so sorry I didn't succeed in using integrate_conditional!

Instead we decided to go for the iterator solution which already works in the master branch.

I guess I will have to implement it this way then. Do you have example code somewhere that demonstrates this? I found this example but I am confused how stop the integration. Can you elaborate on that?

EDIT: Hi I found this useful example on how to use pure iterators on the master branch but I am currently stuck trying on a compiler error:

/usr/include/boost/iterator/iterator_facade.hpp:327:49: error: no matching function for call to ‘implicit_cast(std::pair<const std::array<double, 2ul>&, const double&>*)’ /usr/include/boost/iterator/iterator_facade.hpp:327:49: note: candidate is: /usr/include/boost/implicit_cast.hpp:18:10: note: template T boost::implicit_cast(typename boost::mpl::identity::type)

whenever I try to dereference the iterator. Do you know what might be wrong here?

Apologies if this is the wrong place to ask that.

Thanks

[headmyshoulder]

On 01/06/2014 09:09 PM, Arjun Menon wrote:

Can you send us an short example how you use integrate_conditional?

I actually didn't attempt that yet. I just started porting most of my simulations, so sorry I didn't succeed in using integrate_conditional!

Instead we decided to go for the iterator solution which already
works in the master branch.

I guess I will have to implement it this way then. Do you have example code somewhere that demonstrates this? I found this example https://github.com/headmyshoulder/odeint-v2/blob/master/libs/numeric/odeint/examples/adaptive_iterator.cpp but I am confused how stop the integration. Can you elaborate on that?

The last block shows a use case where the integration stops, or where one tries to find one particular "event":

// 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";

}

This snippet integrates the ODE (the Lorenz system) until x[0]>0.0.

I hope this helps! If it is still unclear just give us a small hint. I think we will write more in the tutorial about this topic very soon.

[arjungm]

Thanks again for the snippet!

I have just one more question. Is there a way to add an observer to make_adaptive_range ? I want to get all the steps in the integration as well.

Maybe I need to use the example in line 160? https://github.com/headmyshoulder/odeint-v2/blob/master/libs/numeric/odeint/examples/const_step_iterator.cpp#L160

Thanks again!

[headmyshoulder]

On 08.01.2014 18:22, Arjun Menon wrote:

Thanks again for the snippet!

I have just one more question. Is there a way to add an observer to make_adaptive_range ? I want to get all the steps in the integration as well.

Hmm, I am not sure right now. Maybe you can implement such an algorithm for yourself

Template< typename Rng , typename Cond , typename F > void iterate_until( range , Cond condition , F func ) { for( auto iter = boost::begin( range ) ; ! condition( iter ) ; ) { func( iter ); } }

Otherwise you could also add the observer into the condition.

Maybe I need to use the example in line 160? https://github.com/headmyshoulder/odeint-v2/blob/master/libs/numeric/odeint/examples/const_step_iterator.cpp#L160

Thanks again!

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-31856293.

[gauteh]

My apologies if this is already supported by the current implementation; it would be nice if the conditional integrator accepted input from the return value of the condition checker so that a distance or direction could be returned. If, for instance, I want to trace a ray untill a current interface the condition could more easily be achieved if the integrator knows whether it should go backwards or forwards, or how far it is from the interface. In MATLAB [0] this is implemented by using a return value (positive or negative) where 0 is condition met (with a given error, of course). In cases where distance or direction does not apply, returning 1 or 0 (condition met) would be sufficient. The integrator should also return from what direction the condition was met.

[0] http://www.mathworks.se/help/matlab/math/ordinary-differential-equations.html#f1-669698

[headmyshoulder]

Your problem sounds like you need a custom controller which controls the ode not with respect for step size adaption and optimal performance but for finding a specific point where a condition is fulfilled.

Maybe an abstraction of a controlled stepper would be sufficient and more elegant for your case? This stepper could adjust the step size in a way such that the surface is not reached but approximated. I think about a wrapping any controlled (or even any regular) stepper by such an controlled stepper. The question is then, how could this be put into a regular integrate function or how could the integration stop? I will think about this.

On 01/09/2014 11:42 AM, Gaute Hope wrote:

My apologies if this is already supported by the current implementation; it would be nice if the conditional integrator accepted input from the return value of the condition checker so that a distance or direction could be returned. If, for instance, I want to trace a ray untill a current interface the condition could more easily be achieved if the integrator knows whether it should go backwards or forwards, or how far it is from the interface. In MATLAB this is implemented by using a return value (positive or negative) where 0 is condition met (with a given error, of course). In cases where distance or direction does not apply, returning 1 or 0 (condition met) would be sufficient. The integrator should also return from what direction the condition was met.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-31919307.

[gauteh]

Yes, in https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-22071174 I wrote an integrator based on the integrate_adaptive code. But the controlling code was all included in the integrating function. This integrator detects if the event happens between two points, then uses the integrator to further approach the event. A similar approach could probably be generalized to solve my use case above for more general problems.

[Rischaard]

For a project I am working on I need to use an iterator-based solution in order to stop the integration when certain conditions are met. It looks like the following code, when applied to our specific situation, will work.

 {
      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";
 }

We are also interested in using the Bulirsch-Stoer method, but when I try to implement the Bulirsch-Stoer stepper as follows:

{
    bulirsch_stoer<state_type> stepper(1e-9, 0.0, 0.0, 0.0);        
    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";
}

I get the following compilation error.

1>c:\program files (x86)\microsoft visual studio 10.0\vc\include\xutility(275): error C2582: 'operator =' function is unavailable in 'boost::numeric::odeint::adaptive_time_iterator<Stepper,System,State,StepperTag>' 1> with 1> [ 1> Stepper=boost::numeric::odeint::bulirsch_stoer, 1> System=lorenz, 1> State=state_type, 1> StepperTag=boost::numeric::odeint::controlled_stepper_tag 1> ] Etc.

Could anyone help me implement the Bulirsch-Stoer stepper in make_adaptive_time_range?

Thanks for the help!

[headmyshoulder]

Hmm, your code works with gcc. Seems like a Visual Studio problem. I don't have visual studio at hand right now, maybe anyone else here?

On 19.01.2015 20:09, Rischaard wrote:

For a project I am working on I need to use an iterator-based solution in order to stop the integration when certain conditions are met. It looks like the following code, when applied to our specific situation, will work.

{ 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 ) , { return ( x.first[0] < 0.0 ); } );

cout << iter->second << "\t" << iter->first[0] << "\t" << iter->first[1] << "\t" << iter->first[2] << "\n";

}

We are also interested in using the Bulirsch-Stoer method, but when I try to implement the Bulirsch-Stoer stepper as follows:

|{ bulirsch_stoer stepper(1e-9, 0.0, 0.0, 0.0);

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";

} |

I get the following compilation error.

1>c:\program files (x86)\microsoft visual studio 10.0\vc\include\xutility(275): error C2582: 'operator =' function is unavailable in 'boost::numeric::odeint::adaptive_time_iterator' 1> with 1> [ 1> Stepper=boost::numeric::odeint::bulirsch_stoer, 1> System=lorenz, 1> State=state_type, 1> StepperTag=boost::numeric::odeint::controlled_stepper_tag 1> ] Etc.

Could anyone help me implement the Bulirsch-Stoer stepper in make_adaptive_time_range?

Thanks for the help!

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/9#issuecomment-70544570.

[mariomulansky]

Hi everyone,

After some discussion I've revisited this problem and added a short example based on dopri5, find_if and a final bisection: 766d2cb0631b89f8f7f5c8187b792a216e691bd9 Quite straight-forward and shouldnt be too difficult to adapt. I hope this is a helpful example for some.