Skip to content
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.

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

Modify the benchmark and tests #1

Merged
merged 5 commits into from
Sep 23, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
1 change: 0 additions & 1 deletion src/NonlinearSolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@ module NonlinearSolve

# DiffEq styled algorithms
export Bisection, Falsi, NewtonRaphson
export ScalarBisection, ScalarNewton

export reinit!
end # module
25 changes: 5 additions & 20 deletions src/scalar.jl
Original file line number Diff line number Diff line change
@@ -1,11 +1,4 @@
"""
ScalarNewton

Fast Newton Raphson for scalar problems.
"""
struct ScalarNewton <: AbstractNonlinearSolveAlgorithm end

function DiffEqBase.solve(prob::NonlinearProblem{uType, false}, ::ScalarNewton, args...; xatol = nothing, xrtol = nothing, maxiters = 1000, kwargs...) where {uType}
function DiffEqBase.solve(prob::NonlinearProblem{<:Number}, ::NewtonRaphson, args...; xatol = nothing, xrtol = nothing, maxiters = 1000, kwargs...)
f = Base.Fix2(prob.f, prob.p)
x = float(prob.u0)
T = typeof(x)
Expand All @@ -15,26 +8,18 @@ function DiffEqBase.solve(prob::NonlinearProblem{uType, false}, ::ScalarNewton,
xo = oftype(x, Inf)
for i in 1:maxiters
fx, dfx = value_derivative(f, x)
iszero(fx) && return x
iszero(fx) && return NewtonSolution(x, :Default)
Δx = dfx \ fx
x -= Δx
if isapprox(x, xo, atol=atol, rtol=rtol)
return x
return NewtonSolution(x, :Default)
end
xo = x
end
return oftype(x, NaN)
return NewtonSolution(x, :MaxitersExceeded)
end

"""
ScalarBisection

Fast Bisection for scalar problems. Note that it doesn't returns exact solution, but returns
the best left limit of the exact solution.
"""
struct ScalarBisection <: AbstractNonlinearSolveAlgorithm end

function DiffEqBase.solve(prob::NonlinearProblem{uType, false}, ::ScalarBisection, args...; maxiters = 1000, kwargs...) where {uType}
function DiffEqBase.solve(prob::NonlinearProblem, ::Bisection, args...; maxiters = 1000, kwargs...)
f = Base.Fix2(prob.f, prob.p)
left, right = prob.u0
fl, fr = f(left), f(right)
Expand Down
27 changes: 17 additions & 10 deletions src/solve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,11 @@ end
function DiffEqBase.init(prob::NonlinearProblem{uType, iip}, alg::AbstractBracketingAlgorithm, args...;
alias_u0 = false,
maxiters = 1000,
# bracketing algorithms only solve scalar problems
immutable = (prob.u0 isa StaticArray || prob.u0 isa Number),
kwargs...
) where {uType, iip}

if !(prob.u0 isa Tuple)
error("You need to pass a tuple of u0 in bracketing algorithms.")
end
Expand Down Expand Up @@ -60,7 +61,7 @@ function DiffEqBase.init(prob::NonlinearProblem{uType, iip}, alg::AbstractNewton
fu = f(u, p)
end


sol = build_newton_solution(u, Val(iip))
if immutable
return NewtonImmutableSolver(1, f, alg, u, fu, p, nothing, false, maxiters, internalnorm, :Default, tol, sol)
Expand All @@ -81,7 +82,7 @@ function DiffEqBase.solve!(solver::AbstractNonlinearSolver)
if solver.iter == solver.maxiters
solver.retcode = :MaxitersExceeded
end
set_solution!(solver)
solver = set_solution(solver)
return solver.sol
end

Expand Down Expand Up @@ -151,19 +152,25 @@ function check_for_exact_solution!(solver::BracketingSolver)
return false
end

function set_solution!(solver::BracketingSolver)
solver.sol.left = solver.left
solver.sol.right = solver.right
solver.sol.retcode = solver.retcode
function set_solution(solver::BracketingSolver)
sol = solver.sol
@set! sol.left = solver.left
@set! sol.right = solver.right
@set! sol.retcode = solver.retcode
@set! solver.sol = sol
return solver
end

function get_solution(solver::BracketingImmutableSolver)
return (left = solver.left, right = solver.right, retcode = solver.retcode)
end

function set_solution!(solver::NewtonSolver)
solver.sol.u = solver.u
solver.sol.retcode = solver.retcode
function set_solution(solver::NewtonSolver)
sol = solver.sol
@set! sol.u = solver.u
@set! sol.retcode = solver.retcode
@set! solver.sol = sol
return solver
end

function get_solution(solver::NewtonImmutableSolver)
Expand Down
2 changes: 1 addition & 1 deletion src/types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ function build_solution(u_prototype, ::Val{false})
return BracketingSolution(zero(u_prototype), zero(u_prototype), :Default)
end

mutable struct NewtonSolution{uType}
struct NewtonSolution{uType}
u::uType
retcode::Symbol
end
Expand Down
39 changes: 24 additions & 15 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,26 +3,35 @@ using StaticArrays
using BenchmarkTools
using Test

function benchmark_immutable()
probN = NonlinearProblem((u,p) -> u .* u .- 2, @SVector[1.0, 1.0])
function benchmark_immutable(f, u0)
probN = NonlinearProblem{false}(f, u0)
solver = init(probN, NewtonRaphson(), immutable = true, tol = 1e-9)
sol = @btime solve!($solver)
@test all(sol.u .* sol.u .- 2 .< 1e-9)
sol = solve!(solver)
end

function benchmark_mutable()
probN = NonlinearProblem((u,p) -> u .* u .- 2, @SVector[1.0, 1.0])
function benchmark_mutable(f, u0)
probN = NonlinearProblem{false}(f, u0)
solver = init(probN, NewtonRaphson(), immutable = false, tol = 1e-9)
sol = @btime (reinit!($solver, $probN); solve!($solver))
@test all(sol.u .* sol.u .- 2 .< 1e-9)
sol = (reinit!(solver, probN); solve!(solver))
end

function benchmark_scalar()
probN = NonlinearProblem((u,p) -> u .* u .- 2, 1.0)
sol = @btime (solve($probN, ScalarNewton()))
@test sol * sol - 2 < 1e-9
function benchmark_scalar(f, u0)
probN = NonlinearProblem{false}(f, u0)
sol = (solve(probN, NewtonRaphson()))
end

benchmark_immutable()
benchmark_mutable()
benchmark_scalar()
f, u0 = (u,p) -> u .* u .- 2, @SVector[1.0, 1.0]
sf, su0 = (u,p) -> u * u - 2, 1.0
sol = benchmark_immutable(f, u0)
@test sol.retcode === :Default
@test all(sol.u .* sol.u .- 2 .< 1e-9)
sol = benchmark_mutable(f, u0)
@test sol.retcode === :Default
@test all(sol.u .* sol.u .- 2 .< 1e-9)
sol = benchmark_scalar(sf, su0)
@test sol.retcode === :Default
@test sol.u * sol.u - 2 < 1e-9

@test (@ballocated benchmark_immutable($f, $u0)) == 0
@test (@ballocated benchmark_mutable($f, $u0)) < 200
@test (@ballocated benchmark_scalar($sf, $su0)) == 0