From 629aa01c60640205538e4f97d31e0b8e09263dab Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 19 Mar 2025 18:25:23 -0400 Subject: [PATCH 1/3] simplify internal_itp --- src/internal_itp.jl | 89 ++++++++++++++++++--------------------------- 1 file changed, 35 insertions(+), 54 deletions(-) diff --git a/src/internal_itp.jl b/src/internal_itp.jl index ebdeb9292..bc7665615 100644 --- a/src/internal_itp.jl +++ b/src/internal_itp.jl @@ -36,78 +36,59 @@ 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.alg_k1) 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 + stats.nsteps += 1 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^k2, 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ϵ + return SciMLBase.build_solution( + prob, alg, xt, f(xt); retcode = ReturnCode.Success, left, right + ) + end yp = f(xp) + stats.nf += 1 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 + ) 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, From 09ea438c5f2f902878aaa0d58df21d40fb0b39a5 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 26 Mar 2025 16:44:42 -0400 Subject: [PATCH 2/3] improvements --- src/internal_itp.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/internal_itp.jl b/src/internal_itp.jl index bc7665615..554b7102d 100644 --- a/src/internal_itp.jl +++ b/src/internal_itp.jl @@ -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...; @@ -41,9 +40,8 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem{IP, Tuple{T, T}}, alg::I return SciMLBase.build_solution(prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right) end - k2 = T(alg.k2) span = abs(right - left) - k1 = T(alg.alg_k1) + k1 = T(alg.scaled_k1)/span n0 = T(alg.n0) n_h = exponent(span / (2 * ϵ)) ϵ_s = ϵ * exp2(n_h + n0) @@ -58,7 +56,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem{IP, Tuple{T, T}}, alg::I x_f = left + span * (fl / (fl - fr)) # Interpolation Step - δ = max(k1 * span^k2, eps(x_f)) + δ = max(k1 * span^2, eps(x_f)) diff = mid - x_f xt = ifelse(δ ≤ abs(diff), x_f + copysign(δ, diff), mid) # Truncation Step From f121a90e59361d084b6f0bdb0ffe01cd925b4781 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 27 Mar 2025 09:57:52 -0400 Subject: [PATCH 3/3] remove stats --- src/internal_itp.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/internal_itp.jl b/src/internal_itp.jl index 554b7102d..8e3c36bd5 100644 --- a/src/internal_itp.jl +++ b/src/internal_itp.jl @@ -49,7 +49,6 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem{IP, Tuple{T, T}}, alg::I i = 1 while i ≤ maxiters - stats.nsteps += 1 span = abs(right - left) mid = (left + right) / 2 r = ϵ_s - (span / 2) @@ -68,7 +67,6 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem{IP, Tuple{T, T}}, alg::I ) end yp = f(xp) - stats.nf += 1 yps = yp * sign(fr) if yps > T0 right, fr = xp, yp