more ITP simplification

Author: oscarddssmith

lib/BracketingNonlinearSolve/src/itp.jl

@@ -49,7 +49,7 @@ end
 function ITP(; scaled_k1::Real = 0.2, k2::Real = 2, n0::Int = 10)
     scaled_k1 < 0 && error("Hyper-parameter κ₁ should not be negative")
     n0 < 0 && error("Hyper-parameter n₀ should not be negative")
-    if k2 < 1 || k2 > (1.5 + sqrt(5) / 2)
+    if !(1 <= k2 <= 1.5 + sqrt(5) / 2)
         throw(ArgumentError("Hyper-parameter κ₂ should be between 1 and 1 + ϕ where \
                              ϕ ≈ 1.618... is the golden ratio"))
     end
@@ -63,22 +63,23 @@ function CommonSolve.solve(
     @assert !SciMLBase.isinplace(prob) "`ITP` only supports out-of-place problems."
 
     f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.tspan
+    left, right = minmax(promote(prob.tspan...)...)
     fl, fr = f(left), f(right)
 
     abstol = NonlinearSolveBase.get_tolerance(
         left, abstol, promote_type(eltype(left), eltype(right))
     )
 
+    stats = SciMLBase.NLStats(2,0,0,0,0)
     if iszero(fl)
         return SciMLBase.build_solution(
-            prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right
+            prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right, stats
         )
     end
 
     if iszero(fr)
         return SciMLBase.build_solution(
-            prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right
+            prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right, stats
         )
     end
 
@@ -87,73 +88,63 @@ function CommonSolve.solve(
             @warn "The interval is not an enclosing interval, opposite signs at the \
                    boundaries are required."
         return SciMLBase.build_solution(
-            prob, alg, left, fl; retcode = ReturnCode.InitialFailure, left, right
+            prob, alg, left, fl; retcode = ReturnCode.InitialFailure, left, right, stats
         )
     end
 
     ϵ = abstol
     k2 = alg.k2
-    k1 = alg.scaled_k1 * abs(right - left)^(1 - k2)
+    span = right - left
+    k1 = alg.scaled_k1 * span^(1 - k2) # k1 > 0
     n0 = alg.n0
-    mid = (left + right) / 2
-    x_f = left + (right - left) * (fl / (fl - fr))
-    xt = left
-    xp = left
-    r = zero(left) # minmax radius
-    δ = zero(left) # truncation error
-    σ = one(mid)
-    n_h = exponent(abs(right - left) / (2 * ϵ))
+    n_h = exponent(span / (2 * ϵ))
     ϵ_s = ϵ * exp2(n_h + n0)
+    T0 = zero(fl)
 
     i = 1
     while i ≤ maxiters
-        span = abs(right - left)
+        stats.nsteps += 1
+        span = right - left
+        mid = (left + right) / 2
         r = ϵ_s - (span / 2)
-        δ = k1 * span^k2
 
-        x_f = left + (right - left) * (fl / (fl - fr))  # Interpolation Step
+        x_f = left + span * (fl / (fl - fr))  # Interpolation Step
 
+        δ = max(k1 * span^k2, eps(x_f))
         diff = mid - x_f
-        σ = sign(diff)
-        xt = ifelse(δ ≤ diff, x_f + σ * δ, mid)  # Truncation Step
 
-        xp = ifelse(abs(xt - mid) ≤ r, xt, mid - σ * r)  # Projection Step
+        xt = ifelse(δ ≤ abs(diff), x_f + copysign(δ, diff), mid)  # Truncation Step
 
-        if abs((left - right) / 2) < ϵ
+        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
+                prob, alg, xt, f(xt); retcode = ReturnCode.Success, left, right, stats
             )
         end
-
-        # update
-        tmin, tmax = minmax(xt, xp)
-        xp ≥ tmax && (xp = prevfloat(tmax))
-        xp ≤ tmin && (xp = nextfloat(tmin))
         yp = f(xp)
+        stats.nf += 1
         yps = yp * sign(fr)
-        T0 = zero(yps)
         if yps > T0
             right, fr = xp, yp
         elseif yps < T0
             left, fl = xp, yp
         else
             return SciMLBase.build_solution(
-                prob, alg, xp, yps; retcode = ReturnCode.Success, left, right
+                prob, alg, xp, yps; retcode = ReturnCode.Success, left, right, stats
             )
         end
 
         i += 1
-        mid = (left + right) / 2
         ϵ_s /= 2
 
-        if Impl.nextfloat_tdir(left, prob.tspan...) == right
+        if nextfloat(left) == right
             return SciMLBase.build_solution(
-                prob, alg, right, fr; retcode = ReturnCode.FloatingPointLimit, left, right
+                prob, alg, right, fr; retcode = ReturnCode.FloatingPointLimit, left, right, stats
             )
         end
     end
 
     return SciMLBase.build_solution(
-        prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right
+        prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right, stats
     )
 end