simplify internal_itp #1127
simplify internal_itp #1127
Conversation
oscardssmith
force-pushed
the
os/simplify-itp
branch
from
e1f8a13 to
09ea438
Compare
| k1 = T(alg.k1) | ||
| k2 = T(alg.k2) | ||
| span = abs(right - left) | ||
| k1 = T(alg.scaled_k1)/span |
There was a problem hiding this comment.
do you know if the lack of scaling based on span was intentional here? The NonlinearSolveBase verison already had it.
|
I think this is ready now (but I have started nuking CI from orbit with ForwardDiff PRs so we'll see if CI here ever completes) |
| xp >= tmax && (xp = prevfloat(tmax)) | ||
| xp <= tmin && (xp = nextfloat(tmin)) | ||
| xp = ifelse(abs(xt - mid) ≤ r, xt, mid - copysign(r, diff)) # Projection Step | ||
| if span < 2ϵ |
There was a problem hiding this comment.
Is this correct? ϵ here is the machine epsilon eps(T) (not eps(left)), i.e. a relative difference, compared to an absolute value span. If left and right are >> 1 in absolute value and different this condition cannot be true, while if they are << 1 in absolute value this condition can be met even if the interval is not tight.
There was a problem hiding this comment.
good catch. This was based on SciML/NonlinearSolve.jl#548 where we have an abstol. I think this branch can just be deleted.
Added for similarity to BracketingNonlinearSolve, but not useful here since we always want ULP accurate solves here. Found by #1127 (comment)
| retcode = ReturnCode.Success, left = left, | ||
| right = right) | ||
| return SciMLBase.build_solution( | ||
| prob, alg, xp, yps; retcode = ReturnCode.Success, left, right |
There was a problem hiding this comment.
Another question: here we're no longer contracting the left - right interval for the solution object, does this change behavior in any case? I'm investigating a small outcome variation in a test in some internal code of ours that opts into rootfind=SciMLBase.RightRootFind for a callback and it seemingly appears between DiffEqBase v6.167.0 and v6.167.2 (though this is just a tentative assessment as of now).
There was a problem hiding this comment.
It's definitely possible that this changes behavior slightly. It's very unclear why the previous code thought that this was trying to do. Specifically, this case only is hit if the function hits exactly 0 (i.e. it will not hit this even if left==right which seems likely unintended).
There was a problem hiding this comment.
Yes, that makes sense; my concern is that while the actual crossing (xp) is provided as the "main" value of the solution sol.u, the way this solution is used ( https://github.com/SciML/DiffEqBase.jl/blob/master/src/callbacks.jl#L365 ) this value is apparently discarded; instead, the left and right fields of the solution seem to be the relevant ones (depending on the user-defined rootfinding side). However, in this file, the solution is built without updating the left and right variables, which could be arbitrarily far from xp (which is I guess the value we want) if the iteration hitting this branch is particularly lucky. The previous code basically used xp (the value producing the exact zero) to update the right value and its prevfloat as the left value, thus producing a tighter interval around the crossing.
There was a problem hiding this comment.
I guess another connected question (that is probably not very relevant practically except for sorting simultaneous events) is what should happen if there are multiple zeros in a row (e.g. for a sequence of adjacent floats for t the condition produces something like -a, 0, 0, 0, a where a is small). It's not as important, because for this to happen you need to have small absolute values of t (so that the expected difference between consecutive values is comparable to numerical error) and/or a flat-ish f, but maybe it's worth thinking about what's the desirable behavior in that case (it should be cheap to seek the first or last zero with a few refinement iterations, probably with straightforward bisection for smaller overhead).
There was a problem hiding this comment.
Upon further testing it appears that setting left and right to xp here reverts my test results to their pre-6.167.1 values, I will make a PR once I double-check.
There was a problem hiding this comment.
Upon further testing it appears that setting left and right to xp here reverts my test results to their pre-6.167.1 values, I will make a PR once I double-check.
to match SciML/NonlinearSolve.jl#548 (although this also fixes
k2to2since paying for a Floating point pow is expensive.Do not merge until SciML/NonlinearSolve.jl#548 is merged.