Skip to content

simplify internal_itp #1127

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Mar 27, 2025
Merged

simplify internal_itp #1127

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
629aa01
simplify internal_itp
Mar 19, 2025
09ea438
improvements
Mar 26, 2025
f121a90
remove stats
Mar 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 35 additions & 58 deletions src/internal_itp.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,11 @@ end
simpler dependencies.
"""
struct InternalITP
k1::Float64
k2::Float64
scaled_k1::Float64
n0::Int
end

InternalITP() = InternalITP(0.007, 1.5, 10)
InternalITP() = InternalITP(0.2, 10)

function SciMLBase.solve(prob::IntervalNonlinearProblem{IP, Tuple{T, T}}, alg::InternalITP,
args...;
Expand All @@ -36,78 +35,56 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem{IP, Tuple{T, T}}, alg::I
ϵ = eps(T)
if iszero(fl)
return SciMLBase.build_solution(prob, alg, left, fl;
retcode = ReturnCode.ExactSolutionLeft, left = left,
right = right)
retcode = ReturnCode.ExactSolutionLeft, left, right)
elseif iszero(fr)
return SciMLBase.build_solution(prob, alg, right, fr;
retcode = ReturnCode.ExactSolutionRight, left = left,
right = right)
retcode = ReturnCode.ExactSolutionRight, left, right)
end
#defining variables/cache
k1 = T(alg.k1)
k2 = T(alg.k2)
span = abs(right - left)
k1 = T(alg.scaled_k1)/span
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do you know if the lack of scaling based on span was intentional here? The NonlinearSolveBase verison already had it.

n0 = T(alg.n0)
n_h = ceil(log2(abs(right - left) / (2 * ϵ)))
mid = (left + right) / 2
x_f = (fr * left - fl * right) / (fr - fl)
xt = left
xp = left
r = zero(left) #minmax radius
δ = zero(left) # truncation error
σ = 1.0
ϵ_s = ϵ * 2^(n_h + n0)
i = 0 #iteration
while i <= maxiters
#mid = (left + right) / 2
n_h = exponent(span / (2 * ϵ))
ϵ_s = ϵ * exp2(n_h + n0)
T0 = zero(fl)

i = 1
while i ≤ maxiters
span = abs(right - left)
mid = (left + right) / 2
r = ϵ_s - (span / 2)
δ = k1 * (span^k2)

## Interpolation step ##
x_f = left + (right - left) * (fl / (fl - fr))
x_f = left + span * (fl / (fl - fr)) # Interpolation Step

## Truncation step ##
σ = sign(mid - x_f)
if δ <= abs(mid - x_f)
xt = x_f + (σ * δ)
else
xt = mid
end
δ = max(k1 * span^2, eps(x_f))
diff = mid - x_f

## Projection step ##
if abs(xt - mid) <= r
xp = xt
else
xp = mid - (σ * r)
end
xt = ifelse(δ ≤ abs(diff), x_f + copysign(δ, diff), mid) # Truncation Step

## Update ##
tmin, tmax = minmax(left, right)
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ϵ
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good catch. This was based on SciML/NonlinearSolve.jl#548 where we have an abstol. I think this branch can just be deleted.

return SciMLBase.build_solution(
prob, alg, xt, f(xt); retcode = ReturnCode.Success, left, right
)
end
yp = f(xp)
yps = yp * sign(fr)
if yps > 0
right = xp
fr = yp
elseif yps < 0
left = xp
fl = yp
if yps > T0
right, fr = xp, yp
elseif yps < T0
left, fl = xp, yp
else
left = prevfloat_tdir(xp, prob.tspan...)
right = xp
return SciMLBase.build_solution(prob, alg, left, f(left);
retcode = ReturnCode.Success, left = left,
right = right)
return SciMLBase.build_solution(
prob, alg, xp, yps; retcode = ReturnCode.Success, left, right
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Copy link
Contributor

@DaniGlez DaniGlez Apr 1, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

)
end

i += 1
mid = (left + right) / 2
ϵ_s /= 2

if nextfloat_tdir(left, prob.tspan...) == right
return SciMLBase.build_solution(prob, alg, left, fl;
retcode = ReturnCode.FloatingPointLimit, left = left,
right = right)
if nextfloat_tdir(left, left, right) == right
return SciMLBase.build_solution(
prob, alg, right, fr; retcode = ReturnCode.FloatingPointLimit, left, right
)
end
end
return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
Expand Down
Loading