Workaround to "interval does not bracket a root f(X) = -0, f(Y) = -0,...7"

Hello, we are using a complex combination of detectors during our propagations. We often encounter an hipparchus exception “interval does not bracket a root” with one of two values which is 0 or -0 and the other of the same sign. Looking at BracketingNthOrderBrentSolver.doSolveInterval we can see that a strict comparison to 0 is performed. We introduced an epsilon for the comparison and the error is not raised anymore.
There is 2 ways of introducing this epsilon :

  • just before the exception is raised, check if one value is close to zero and if ot’s ok, return it. It does not modifiy the current algorithm except when an exception occurs. Code here : BracketingNthOrderBrentSolver.java (16.4 KB)
  • or in the algorithm, when comparing y[…] to 0. This could change slightly the behaviour, maybe by making it more faster but I don’t know if the impact can be bigger. Code here : BracketingNthOrderBrentSolver.java (15.9 KB)

So… here are my questions :

  • any reason to not do that ?
  • is one of the solution better than the other. I would prefer the second one but it can have an impact
  • the epsilon is arbitrary set to 1e-12 , is there a value that is already used somewhere which would be more significant or configurable ?
  • can this behaviour be included in a future release of hipparchus ? should I contact directly hipparchus ?

Thanks

I am on the fence on this. Looking at your proposal, it seems you want to accept a root when in fact it really does not belong to the interval. The fact the exception message shows 0 (or -0) does not mean the value is exactly zero, but rather a small number not displayed (like +1.0e-15 or something like that). The tests on the boundary are written in such a way they match “perfect root” (this is what is written in the comment), and by “perfect” we mean exact 0, not close to zero within some tolerance.

any reason to not do that ?

Yes, there are reasons not to do that, and these reasons are in fact related to event detection. In many cases, events are set up to separate processing of the “before” and “after” phases. This sometimes implies that a computation is changed and restarted or reset at the detected event, which implies that one of the boundary for the next phase (the after phase if propagating forward, the before phase if propagating backward) becomes by construction one of the detected zero. So basically, just when we restart, we are already at zero. If we consider that when we are close to zero we ignore the sign and consider we have found a root, then we will stall at the zero and never cross it. This was exactly what caused years ago the unstability on events detection we corrected since then. Now, we take great care to make sure we cross the zero and once we have crossed it we don’t want to detect the same zero again. If for example we know we have already considered a transition from negative to positive, if we start the next phase with a value y(t1) = +1.0e-300 and y(t2) = +0.2, then we consider this interval does not bracket a root. The root was in the previous interval and we consider a root at the boundary between two intervals is in one of the intervals only.

the epsilon is arbitrary set to 1e-12 , is there a value that is already used somewhere which would be more significant or configurable ?

There is already a configurable value that has this meaning. You can get it by calling getFunctionValueAccuracy(), the default value is 1.0e-14 and it can be set to another value by using one of the constructors.

can this behaviour be included in a future release of hipparchus ? should I contact directly hipparchus ?

This has to be discussed, because as I explained before, there is a reason for this very strict handling of bracketing intervals. For discussing this, yes the Hipparchus project is the place to go. You should open an issue here https://github.com/Hipparchus-Math/hipparchus/issues. You will stay in well-known community since the two Hipparchus developers that are more familiar with events detection and root solving are @evan.ward and myself (and I designed the algorithm behind BracketingNthOrderBrentSolver).

I would like to come back to your original problem. You say you encounter this exception quite often so there is definitely a problem there. Normally, as long as the function is really continuous in the neighbourhood of roots, the logic between EventState and the solver is very stable by now. With functions continuous around zero crossings, the solver should not be called with an interval where the zero is just before or after the boundaries, the solver should either have been called for the interval before or it will be called for the interval after, where the event really lies. We even handle properly functions that just touch zero without crossing it (and which don’t trigger events as there is no sign change). In many cases, problems in events detection were in fact problems in the function definition and coupling between the function definition and the dynamics changes triggered by the event. In most of these cases, the function was not really continuous and a zero crossing triggered by the before event logic changed the behaviour in the after event logic (say you switch off a device, or a force, or change an attitude) and as a side effect of this change the event function definition changed it sign again, hence breaking continuity. Are you sure this is not what happens in your case?

Thank you again for this complete answer. Starting by the end, I think that the root of the problem is the way we use detectors for modelling system : event detections change the g function by side effect (several detectors and systems combined). I don’t know if there is a better approach to do what we want.
Back to the hipparchus problem, I understand the risk of changing the algorithm, and to make old bugs come back… Definietely not a good idea. So I will just use the workaround to avoid the exception rise if we are close enough to the root. I will create an issue on the hipparchus to discuss this specific point.

Hi Mikael,

Good question, the event detection code can be tricky to understand and I’m very interested in any bugs you find in that code. Using complex event detectors that interact is supported, as long as the “rules” of event detection at [1] are followed. Specifically the definition of a g function can only change when an event is triggered. The event handler must return RESET_STATE or RESET_DERIVATIVES (or RESET_EVENTS in Orekit 10.0) if it changes the definition of the g function of another event detector.[2] This forces the event detection logic to recheck the g function of all the other event detectors. Often the error you are seeing is caused by an event detector that returns CONTINUE even though it changes the definition of the g function of another event detector.

One takeaway I have from this discussion is that we should improve the error message, because the current message about +/-0 not bracketing zero is misleading.

If you still believe it is an error could you attach a test case the demonstrates the error? I’ll need the test case to dissect exactly what went wrong and I’ll add it to our test suite to make sure the bug stays fixed.

Regards,
Evan

[1] https://hipparchus.org/apidocs/org/hipparchus/ode/events/package-summary.html
[2] https://hipparchus.org/apidocs/org/hipparchus/ode/events/ODEEventHandler.html#g-org.hipparchus.ode.ODEStateAndDerivative-

Hi Evan,
Thank you for this answer. I think you found our error, we do not use RESET_XXX.I will make the tests next week and I will let you know about the result.
Regards,

I made a check on our detectors and some RESET_DERIVATIVES where missing. It seems to fix the problem and the overload of hipparchus is not used anymore. Thanks a lot.