feat: implement all bracketing algorithms

Author: avik-pal

lib/SimpleBracketingNonlinearSolve/Project.toml renamed to lib/BracketingNonlinearSolve/Project.toml

@@ -1,13 +1,19 @@
-name = "SimpleBracketingNonlinearSolve"
+name = "BracketingNonlinearSolve"
 uuid = "70df07ce-3d50-431d-a3e7-ca6ddb60ac1e"
 authors = ["Avik Pal <[email protected]> and contributors"]
 version = "1.0.0"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [compat]
 CommonSolve = "0.2.4"
+ConcreteStructs = "0.2.3"
+NonlinearSolveBase = "1"
+PrecompileTools = "1.2.1"
 SciMLBase = "2.50"
 julia = "1.10"

lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl

@@ -0,0 +1,38 @@
+module BracketingNonlinearSolve
+
+using ConcreteStructs: @concrete
+
+using CommonSolve: CommonSolve
+using NonlinearSolveBase: NonlinearSolveBase
+using SciMLBase: SciMLBase, AbstractNonlinearAlgorithm, IntervalNonlinearProblem, ReturnCode
+
+using PrecompileTools: @compile_workload, @setup_workload
+
+abstract type AbstractBracketingAlgorithm <: AbstractNonlinearAlgorithm end
+
+include("common.jl")
+
+include("alefeld.jl")
+include("bisection.jl")
+include("brent.jl")
+include("falsi.jl")
+include("itp.jl")
+include("ridder.jl")
+
+@setup_workload begin
+    for T in (Float32, Float64)
+        prob_brack = IntervalNonlinearProblem{false}(
+            (u, p) -> u^2 - p, T.((0.0, 2.0)), T(2))
+        algs = (Alefeld(), Bisection(), Brent(), Falsi(), ITP(), Ridder())
+
+        @compile_workload begin
+            for alg in algs
+                CommonSolve.solve(prob_brack, alg; abstol = 1e-6)
+            end
+        end
+    end
+end
+
+export Alefeld, Bisection, Brent, Falsi, ITP, Ridder
+
+end

lib/BracketingNonlinearSolve/src/alefeld.jl

@@ -0,0 +1,135 @@
+"""
+    Alefeld()
+
+An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145/210089.210111).
+
+The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than
+algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
+"""
+struct Alefeld <: AbstractBracketingAlgorithm end
+
+function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Alefeld, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
+    f = Base.Fix2(prob.f, prob.p)
+    a, b = prob.tspan
+    c = a - (b - a) / (f(b) - f(a)) * f(a)
+
+    fc = f(c)
+    if a == c || b == c
+        return SciMLBase.build_solution(
+            prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit, left = a, right = b)
+    end
+
+    if iszero(fc)
+        return SciMLBase.build_solution(
+            prob, alg, c, fc; retcode = ReturnCode.Success, left = a, right = b)
+    end
+
+    a, b, d = Impl.bracket(f, a, b, c)
+    e = zero(a)   # Set e as 0 before iteration to avoid a non-value f(e)
+
+    for i in 2:maxiters
+        # The first bracketing block
+        f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
+        if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
+            c = Impl.newton_quadratic(f, a, b, d, 2)
+        else
+            c = Impl.ipzero(f, a, b, d, e)
+            if (c - a) * (c - b) ≥ 0
+                c = Impl.newton_quadratic(f, a, b, d, 2)
+            end
+        end
+
+        ē, fc = d, f(c)
+        if a == c || b == c
+            return SciMLBase.build_solution(
+                prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+                left = a, right = b)
+        end
+
+        if iszero(fc)
+            return SciMLBase.build_solution(
+                prob, alg, c, fc; retcode = ReturnCode.Success, left = a, right = b)
+        end
+
+        ā, b̄, d̄ = Impl.bracket(f, a, b, c)
+
+        # The second bracketing block
+        f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
+        if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
+            c = Impl.newton_quadratic(f, ā, b̄, d̄, 3)
+        else
+            c = Impl.ipzero(f, ā, b̄, d̄, ē)
+            if (c - ā) * (c - b̄) ≥ 0
+                c = Impl.newton_quadratic(f, ā, b̄, d̄, 3)
+            end
+        end
+        fc = f(c)
+
+        if ā == c || b̄ == c
+            return SciMLBase.build_solution(
+                prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+                left = ā, right = b̄)
+        end
+
+        if iszero(fc)
+            return SciMLBase.build_solution(
+                prob, alg, c, fc; retcode = ReturnCode.Success, left = ā, right = b̄)
+        end
+
+        ā, b̄, d̄ = Impl.bracket(f, ā, b̄, c)
+
+        # The third bracketing block
+        u = ifelse(abs(f(ā)) < abs(f(b̄)), ā, b̄)
+        c = u - 2 * (b̄ - ā) / (f(b̄) - f(ā)) * f(u)
+        if (abs(c - u)) > 0.5 * (b̄ - ā)
+            c = 0.5 * (ā + b̄)
+        end
+        fc = f(c)
+
+        if ā == c || b̄ == c
+            return SciMLBase.build_solution(
+                prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+                left = ā, right = b̄)
+        end
+
+        if iszero(fc)
+            return SciMLBase.build_solution(
+                prob, alg, c, fc; retcode = ReturnCode.Success, left = ā, right = b̄)
+        end
+
+        ā, b̄, d = Impl.bracket(f, ā, b̄, c)
+
+        # The last bracketing block
+        if b̄ - ā < 0.5 * (b - a)
+            a, b, e = ā, b̄, d̄
+        else
+            e = d
+            c = 0.5 * (ā + b̄)
+            fc = f(c)
+
+            if ā == c || b̄ == c
+                return SciMLBase.build_solution(
+                    prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+                    left = ā, right = b̄)
+            end
+            if iszero(fc)
+                return SciMLBase.build_solution(
+                    prob, alg, c, fc; retcode = ReturnCode.Success, left = ā, right = b̄)
+            end
+            a, b, d = Impl.bracket(f, ā, b̄, c)
+        end
+    end
+
+    # Reassign the value a, b, and c
+    if b == c
+        b = d
+    elseif a == c
+        a = d
+    end
+    fc = f(c)
+
+    # Reuturn solution when run out of max interation
+    return SciMLBase.build_solution(
+        prob, alg, c, fc; retcode = ReturnCode.MaxIters, left = a, right = b)
+end

lib/BracketingNonlinearSolve/src/bisection.jl

@@ -0,0 +1,81 @@
+"""
+    Bisection(; exact_left = false, exact_right = false)
+
+A common bisection method.
+
+### Keyword Arguments
+
+  - `exact_left`: whether to enforce whether the left side of the interval must be exactly
+    zero for the returned result. Defaults to false.
+  - `exact_right`: whether to enforce whether the right side of the interval must be exactly
+    zero for the returned result. Defaults to false.
+
+!!! danger "Keyword Arguments"
+
+    Currently, the keyword arguments are not implemented.
+"""
+@kwdef struct Bisection <: AbstractBracketingAlgorithm
+    exact_left::Bool = false
+    exact_right::Bool = false
+end
+
+function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Bisection,
+        args...; maxiters = 1000, abstol = nothing, kwargs...)
+    @assert !SciMLBase.isinplace(prob) "`Bisection` only supports out-of-place problems."
+
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+
+    abstol = NonlinearSolveBase.get_tolerance(
+        abstol, promote_type(eltype(left), eltype(right)))
+
+    if iszero(fl)
+        return SciMLBase.build_solution(
+            prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right)
+    end
+
+    if iszero(fr)
+        return SciMLBase.build_solution(
+            prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right)
+    end
+
+    i = 1
+    while i ≤ maxiters
+        mid = (left + right) / 2
+
+        if mid == left || mid == right
+            return SciMLBase.build_solution(
+                prob, alg, left, fl; retcode = ReturnCode.FloatingPointLimit, left, right)
+        end
+
+        fm = f(mid)
+        if abs((right - left) / 2) < abstol
+            return SciMLBase.build_solution(
+                prob, alg, mid, fm; retcode = ReturnCode.Success, left, right)
+        end
+
+        if iszero(fm)
+            right = mid
+            break
+        end
+
+        if sign(fl) == sign(fm)
+            fl = fm
+            left = mid
+        else
+            fr = fm
+            right = mid
+        end
+
+        i += 1
+    end
+
+    sol, i, left, right, fl, fr = Impl.bisection(
+        left, right, fl, fr, f, abstol, maxiters - i, prob, alg)
+
+    sol !== nothing && return sol
+
+    return SciMLBase.build_solution(
+        prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right)
+end

lib/BracketingNonlinearSolve/src/brent.jl

@@ -0,0 +1,109 @@
+"""
+    Brent()
+
+Left non-allocating Brent method.
+"""
+struct Brent <: AbstractBracketingAlgorithm end
+
+function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Brent, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
+    @assert !SciMLBase.isinplace(prob) "`Brent` only supports out-of-place problems."
+
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+    ϵ = eps(convert(typeof(fl), 1))
+
+    abstol = NonlinearSolveBase.get_tolerance(
+        abstol, promote_type(eltype(left), eltype(right)))
+
+    if iszero(fl)
+        return SciMLBase.build_solution(
+            prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right)
+    end
+
+    if iszero(fr)
+        return SciMLBase.build_solution(
+            prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right)
+    end
+
+    if abs(fl) < abs(fr)
+        left, right = right, left
+        fl, fr = fr, fl
+    end
+
+    c = left
+    d = c
+    i = 1
+    cond = true
+
+    while i < maxiters
+        fc = f(c)
+
+        if fl != fc && fr != fc
+            # Inverse quadratic interpolation
+            s = left * fr * fc / ((fl - fr) * (fl - fc)) +
+                right * fl * fc / ((fr - fl) * (fr - fc)) +
+                c * fl * fr / ((fc - fl) * (fc - fr))
+        else
+            # Secant method
+            s = right - fr * (right - left) / (fr - fl)
+        end
+
+        if (s < min((3 * left + right) / 4, right) ||
+            s > max((3 * left + right) / 4, right)) ||
+           (cond && abs(s - right) ≥ abs(right - c) / 2) ||
+           (!cond && abs(s - right) ≥ abs(c - d) / 2) ||
+           (cond && abs(right - c) ≤ ϵ) ||
+           (!cond && abs(c - d) ≤ ϵ)
+            # Bisection method
+            s = (left + right) / 2
+            if s == left || s == right
+                return SciMLBase.build_solution(prob, alg, left, fl;
+                    retcode = ReturnCode.FloatingPointLimit, left, right)
+            end
+            cond = true
+        else
+            cond = false
+        end
+
+        fs = f(s)
+        if abs((right - left) / 2) < abstol
+            return SciMLBase.build_solution(prob, alg, s, fs;
+                retcode = ReturnCode.Success, left, right)
+        end
+
+        if iszero(fs)
+            if right < left
+                left = right
+                fl = fr
+            end
+            right = s
+            fr = fs
+            break
+        end
+
+        if fl * fs < 0
+            d, c, right = c, right, s
+            fr = fs
+        else
+            left = s
+            fl = fs
+        end
+
+        if abs(fl) < abs(fr)
+            d = c
+            c, right, left = right, left, c
+            fc, fr, fl = fr, fl, fc
+        end
+        i += 1
+    end
+
+    sol, i, left, right, fl, fr = Impl.bisection(
+        left, right, fl, fr, f, abstol, maxiters - i, prob, alg)
+
+    sol !== nothing && return sol
+
+    return SciMLBase.build_solution(
+        prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right)
+end