|
| 1 | +@testset "Root-finding full example" begin |
| 2 | + # Problem |
| 3 | + """ |
| 4 | + RootFindingProblem |
| 5 | +
|
| 6 | + Struct for problem ``f(x) = 0``. |
| 7 | + """ |
| 8 | + mutable struct RootFindingProblem{T} |
| 9 | + f |
| 10 | + x₀::T |
| 11 | + evals::Int |
| 12 | + end |
| 13 | + RootFindingProblem(f, x₀::T) where {T} = RootFindingProblem{T}(f, x₀, 0) |
| 14 | + function (rfp::RootFindingProblem)(x) |
| 15 | + rfp.evals += 1 |
| 16 | + rfp.f(x) |
| 17 | + end |
| 18 | + function SolverCore.reset!(rfp::RootFindingProblem{T}) where {T} |
| 19 | + rfp.evals = 0 |
| 20 | + end |
| 21 | + |
| 22 | + # Output |
| 23 | + mutable struct RFPOutput{T} <: AbstractSolverOutput{T} |
| 24 | + status::Symbol |
| 25 | + solution::T |
| 26 | + fx::T |
| 27 | + evals::Int |
| 28 | + elapsed_time::Float64 |
| 29 | + end |
| 30 | + function RFPOutput( |
| 31 | + status::Symbol, |
| 32 | + x::T; |
| 33 | + fx::T = T(Inf), |
| 34 | + evals::Int = -1, |
| 35 | + elapsed_time::Float64 = Inf, |
| 36 | + ) where {T} |
| 37 | + RFPOutput{T}(status, x, fx, evals, elapsed_time) |
| 38 | + end |
| 39 | + |
| 40 | + # Solver |
| 41 | + mutable struct Bissection{T} <: AbstractSolver{T} |
| 42 | + initialized::Bool |
| 43 | + params::Dict |
| 44 | + workspace |
| 45 | + end |
| 46 | + |
| 47 | + function Bissection( |
| 48 | + rfp::RootFindingProblem{T}; |
| 49 | + θ = one(T) / 2, |
| 50 | + δ = one(T), |
| 51 | + explorer = true, |
| 52 | + explorer_tries = 3, |
| 53 | + ) where {T} |
| 54 | + Bissection{T}(true, Dict(:θ => θ, :δ => δ, :explorer => true, :explorer_tries => 3), []) |
| 55 | + end |
| 56 | + |
| 57 | + SolverCore.parameters(::Type{Bissection{T}}) where {T} = ( |
| 58 | + θ = (default = one(T) / 2, type = T, scale = :linear, min = T(0.1), max = T(0.9)), |
| 59 | + δ = (default = one(T), type = T, scale = :log, min = √eps(T), max = T(10)), |
| 60 | + explorer = (default = true, type = Bool), |
| 61 | + explorer_tries = (default = 3, type = Int, min = 1, max = 3), |
| 62 | + ) |
| 63 | + function SolverCore.are_valid_parameters(::Type{Bissection{T}}, θ, δ, _, _2) where {T} |
| 64 | + return (0 ≤ θ ≤ 1) && δ > 0 |
| 65 | + end |
| 66 | + |
| 67 | + function SolverCore.solve!( |
| 68 | + solver::Bissection, |
| 69 | + f::RootFindingProblem; |
| 70 | + atol = 1e-3, |
| 71 | + rtol = 1e-3, |
| 72 | + kwargs..., |
| 73 | + ) |
| 74 | + for (k, v) in kwargs |
| 75 | + solver.params[k] = v |
| 76 | + end |
| 77 | + |
| 78 | + θ = solver.params[:θ] |
| 79 | + δ = solver.params[:δ] |
| 80 | + explorer = solver.params[:explorer] |
| 81 | + explorer_tries = solver.params[:explorer_tries] |
| 82 | + a, b = f.x₀ - δ, f.x₀ + δ |
| 83 | + t₀ = time() |
| 84 | + |
| 85 | + fa, fb = f(a), f(b) |
| 86 | + ϵ = atol + rtol * (abs(fa) + abs(fb)) / 2 |
| 87 | + solved = false |
| 88 | + if abs(fa) < ϵ |
| 89 | + return RFPOutput(:acceptable, a, evals = f.evals, fx = fa, elapsed_time = time() - t₀) |
| 90 | + elseif abs(fb) < ϵ |
| 91 | + return RFPOutput(:acceptable, b, evals = f.evals, fx = fb, elapsed_time = time() - t₀) |
| 92 | + elseif fa * fb > 0 |
| 93 | + error("Increasing coverage") |
| 94 | + end |
| 95 | + x = a * θ + b * (1 - θ) |
| 96 | + fx = f(x) |
| 97 | + while abs(fx) ≥ ϵ && f.evals < 100 |
| 98 | + (a, b, fa, fb) = fa * fx < 0 ? (a, x, fa, fx) : (x, b, fx, fb) |
| 99 | + x = a * θ + b * (1 - θ) |
| 100 | + fx = f(x) |
| 101 | + if explorer |
| 102 | + for k = 1:explorer_tries |
| 103 | + r = rand() |
| 104 | + xt = a * r + b * (1 - r) |
| 105 | + ft = f(xt) |
| 106 | + if abs(ft) < abs(fx) |
| 107 | + x, fx = xt, ft |
| 108 | + end |
| 109 | + end |
| 110 | + end |
| 111 | + end |
| 112 | + status = :unknown |
| 113 | + if abs(fx) < ϵ |
| 114 | + status = :acceptable |
| 115 | + elseif f.evals ≥ 100 |
| 116 | + status = :max_eval |
| 117 | + end |
| 118 | + return RFPOutput(status, x, evals = f.evals, fx = fx, elapsed_time = time() - t₀) |
| 119 | + end |
| 120 | + |
| 121 | + # Testing |
| 122 | + rfp_problems = [ |
| 123 | + RootFindingProblem(x -> x - π, 3.0), |
| 124 | + RootFindingProblem(x -> x^3 - 2, 1.0), |
| 125 | + RootFindingProblem(x -> x * exp(x) - 1, 1.0), |
| 126 | + RootFindingProblem(x -> x / (1 + x^2) - 0.5, 0.0), |
| 127 | + RootFindingProblem(x -> 1 / (1 + exp(-x)) - 0.5, 1.0), |
| 128 | + ] |
| 129 | + |
| 130 | + for f in rfp_problems |
| 131 | + solver = Bissection(f) |
| 132 | + output = solve!(solver, f) |
| 133 | + @test abs(output.fx) < 1e-2 |
| 134 | + |
| 135 | + solver = Bissection(f, θ = 0.1) |
| 136 | + output = solve!(solver, f) |
| 137 | + @test abs(output.fx) < 1e-2 |
| 138 | + end |
| 139 | + |
| 140 | + # Tuning |
| 141 | + Random.seed!(0) |
| 142 | + grid = with_logger(NullLogger()) do |
| 143 | + grid_search_tune( |
| 144 | + Bissection, |
| 145 | + rfp_problems, |
| 146 | + success = o -> o.status == :acceptable, |
| 147 | + costs = [(o -> o.elapsed_time, 100.0), (o -> o.evals, 100)], |
| 148 | + grid_length = 9, |
| 149 | + δ = 10.0 .^ (-5:1), |
| 150 | + ) |
| 151 | + end |
| 152 | + best = sort(grid[:], by = x -> x[2][2])[1][1] |
| 153 | + @test best[1] ≈ 0.3 |
| 154 | + @test best[2] ≈ 1 |
| 155 | + @test !best[3] |
| 156 | + @test best[4] == 1 |
| 157 | +end |
0 commit comments