Skip to content

Commit

Permalink
Merge 1a6bc0f into 28a6bd0
Browse files Browse the repository at this point in the history
  • Loading branch information
@jverzani
jverzani committed Aug 15, 2018
2 parents 28a6bd0 + 1a6bc0f commit 234bb3a
Show file tree
Hide file tree
Showing 6 changed files with 193 additions and 165 deletions.
2 changes: 1 addition & 1 deletion src/Roots.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ else
using Missings
end

using Compat: @nospecialize, lastindex, range
using Compat: @nospecialize, lastindex, range, Nothing


export fzero,
Expand Down
11 changes: 6 additions & 5 deletions src/bracketing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -623,7 +623,7 @@ function assess_convergence(method::AbstractAlefeldPotraShi, state::UnivariateZe
u,fu = choose_smallest(state.xn0, state.xn1, state.fxn0, state.fxn1)
u, fu = choose_smallest(u, state.m[1], fu, state.fm[1])

if abs(fu) <= max(options.abstol, abs(u) * options.reltol)
if abs(fu) <= maximum(promote(options.abstol, abs(u) * oneunit(fu) / oneunit(u) * options.reltol))
state.f_converged = true
state.xn1=u
state.fxn1=fu
Expand All @@ -634,7 +634,7 @@ function assess_convergence(method::AbstractAlefeldPotraShi, state::UnivariateZe
end

a,b = state.xn0, state.xn1
tol = max(options.xabstol, max(abs(a),abs(b)) * options.xreltol)
tol = maximum(promote(options.xabstol, max(abs(a),abs(b)) * options.xreltol))

if abs(b-a) <= 2tol
# use smallest of a,b,m
Expand Down Expand Up @@ -943,28 +943,29 @@ FalsePosition(x=:anderson_bjork) = FalsePosition{x}()
function update_state(method::FalsePosition, fs, o::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}

a::T, b::T = o.xn0, o.xn1

fa::S, fb::S = o.fxn0, o.fxn1

lambda = fb / (fb - fa)
tau = eps(T)^(2/3) # some engineering to avoid short moves
tau = 1e-10 # some engineering to avoid short moves; still fails on some
if !(tau < abs(lambda) < 1-tau)
lambda = 1/2
end

x::T = b - lambda * (b-a)
fx::S = fs(x)
incfn(o)

if iszero(fx)
o.xn1 = x
o.fxn1 = fx
o.f_converged = true
return
end

if sign(fx)*sign(fb) < 0
a, fa = b, fb
else
fa = galdino_reduction(method, fa, fb, fx) #galdino[method.reduction_factor](fa, fb, fx)
fa = galdino_reduction(method, fa, fb, fx)
end
b, fb = x, fx

Expand Down
149 changes: 7 additions & 142 deletions src/derivative_free.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,150 +64,15 @@ to 3 times.
"""
struct Order0 <: AbstractSecant end

##################################################


## order 0
# goal: more robust to initial guess than higher order methods
# follows roughly algorithm described http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf, the SOLVE button from the HP-34C
# though some modifications were made.
# * use secant step
# * if along the way a bracket is found, switch to bracketing method.
# * if secant step fails to decrease, we use quadratic step up to 3 times
#
# Goal *was* to return a value `x` with either. This can be done if `Bisection` is used
# but for now we opt for a bit speedier using AlefeldPotraShi
# * `f(x) == 0.0` or
# * `f(prevfloat(x)) * f(nextfloat(x)) < 0`.
# if a bracket is found that can be done, otherwise secant step is used
# init_state/init_options try to call Order1 (AbstractSecant) methods with modifications
# we slip in quad_ctr into a value and keep xn1 with smaller norm
function init_state(M::Order0, f, x::Union{Tuple, Vector})
x0, x1 = promote(float.(x)...)
fx0, fx1 = promote(f(x0), f(x1))

# we keep xn1, fxn1 smallest
a,b,fa, fb = sort_smallest(x0, x1, fx0, fx1)


T, S = eltype(x1), eltype(fx1)
quad_ctr = one(x1)


state = UnivariateZeroState(b, a, [one(T)], ## x1, x0, quad_ctr
fb, fa, S[], ## fx1, fx0, fc, mflag
0, 2,
false, false, false, false,
"")
state
end

function init_state!(state::UnivariateZeroState{T,S}, M::Order0, f, x::Union{Tuple,Vector}) where {T,S}
x0::T, x1::T = x
fx0::S, fx1::S = promote(f(x0), f(x1))
a,b,fa, fb = sort_smallest(x0, x1, fx0, fx1)
quad_ctr = one(x1)
init_state!(state, b, a, T[quad_ctr], fb, fa, S[])
state.fnevals = 2
nothing
end



function init_options(::Order0,
state::UnivariateZeroState{T,S};
maxevals = nothing,
kwargs...) where {T, S}

# same as order 1 save maxevals
options = init_options(Order1(), state; kwargs...)
options.maxevals = maxevals == nothing ? 5 * ceil(Int, -log(eps(T))) : maxevals
options
end


function assess_convergence(method::Order0, state::UnivariateZeroState{T,S}, options) where {T,S}

quad_ctr = state.m[1] #add in test on quad qtr
if quad_ctr > 3
state.stopped = true
return true
end

# how to call super?
assess_convergence(Order1(), state, options)

end


function run_bisection(f, xs, state, options)
#M = Bisection() # exact for floating point
M = AlefeldPotraShi() # *usually* exact
#M = Brent() # a bit faster, but not always convergent, as implemented (cf. RootTesting)
steps = state.steps
init_state!(state, M, f, xs); state.steps += steps
init_options!(options, M)
find_zero(M, f, options, state)
a, b = xs
u,v = a > b ? (b, a) : (a, b)
state.message = "Bisection used over ($u, $v), steps not shown"
return nothing
end

# main algorithm
function update_state(method::Order0, f, state::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}

a::T, b::T, quad_ctr = state.xn0, state.xn1, state.m[1]
fa::S, fb::S = state.fxn0, state.fxn1

## we keep fb < fa Could do in init_state to save a conditional XXX
if abs(fa) < abs(fb)
a, b= b, a
fa,fb = fb, fa
end

isbracket(fa, fb) && return run_bisection(f, (a, b), state, options)

gamma::T = secant_step(a,b,fa,fb)
# modify if gamma is too small or too big
if isinf(gamma) || isnan(gamma) || iszero(abs(gamma-b))
gamma = b + sign(gamma-b) * 1/1000 * abs(b-a) # too small
elseif abs(gamma-b) >= 100 * abs(b-a)
gamma = b + sign(gamma-b) * 100 * abs(b-a) ## too big
end
fgamma::T = f(gamma)
incfn(state)


isbracket(fgamma, fb) && return run_bisection(f, (gamma, b), state, options)


# decreasing is good
if abs(fgamma) < abs(fb)
state.xn0, state.xn1 = b, gamma
state.fxn0, state.fxn1 = fb, fgamma
return nothing
end


# try quad step
gamma = quad_vertex(a,fa,b, fb, gamma, fgamma)
fgamma = f(gamma)
incfn(state)
state.m[1] += 1 # increment quad_ctr

if abs(fgamma) < abs(fb)
state.xn0, state.xn1 = b, gamma
state.fxn0, state.fxn1 = fb, fgamma
else
state.xn0 = gamma
state.fxn0 = fgamma
end

return nothing
function find_zero(fs, x0, method::Order0;
tracks::AbstractTracks=NullTracks(),
verbose=false,
kwargs...)
M = Order1()
N = AlefeldPotraShi()
_find_zero(fs, x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
end


##################################################

## Secant
Expand Down

0 comments on commit 234bb3a

Please sign in to comment.