Skip to content

Commit 9f212fb

Browse files
lbenetdpsanders
lbenet
authored and
dpsanders
committed
Move general methods of IntervalBox to IntervalArithmetic (#83)
* Remove src/bisect.jl and its tests * Remove guarded_mid (not used) and update where usage * Remove functs moved to IntervalArithmetic, and update where's usage * Update REQUIRE * Make IntervalLike{T} and NewtonLike constants * Fix problems with tests with current master of IntervalArithmetics * Minor cleanup * Fixes
1 parent 25ebe4c commit 9f212fb

File tree

9 files changed

+56
-153
lines changed

9 files changed

+56
-153
lines changed

‎REQUIRE

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
julia 0.6
2-
IntervalArithmetic 0.12
3-
ForwardDiff 0.3
4-
StaticArrays 0.5
2+
IntervalArithmetic 0.13
3+
ForwardDiff 0.7.5
4+
StaticArrays 0.7.1
55
Polynomials

‎src/IntervalRootFinding.jl

Lines changed: 11 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -30,26 +30,10 @@ import IntervalArithmetic: interval, wideinterval
3030
const derivative = ForwardDiff.derivative
3131
const D = derivative
3232

33-
# Common functionality:
34-
35-
doc"""Returns the midpoint of the interval x, slightly shifted in case
36-
the midpoint is an exact root"""
37-
function guarded_mid{T}(f::Function, x::Interval{T})
38-
m = mid(x)
39-
40-
if f(m) == zero(T) # midpoint exactly a root
41-
α = convert(T, 0.46875) # close to 0.5, but exactly representable as a floating point
42-
m = α*x.lo + (one(T)-α)*x.hi # displace to another point in the interval
43-
end
44-
45-
@assert m ∈ x
46-
47-
return m
48-
end
33+
const where_bisect = IntervalArithmetic.where_bisect ## 127//256
4934

5035

5136
include("root_object.jl")
52-
include("bisect.jl")
5337

5438
include("newton.jl")
5539
include("krawczyk.jl")
@@ -61,22 +45,25 @@ include("newton1d.jl")
6145
include("quadratic.jl")
6246

6347

64-
gradient(f) = X -> ForwardDiff.gradient(f, SVector(X))
48+
gradient(f) = X -> ForwardDiff.gradient(f, X[:])
49+
50+
ForwardDiff.jacobian(f, X::IntervalBox) = ForwardDiff.jacobian(f, X.v)
51+
6552
const ∇ = gradient
6653

6754
export ∇
6855

6956

7057

7158

72-
function find_roots{T}(f::Function, a::Interval{T}, method::Function = newton;
73-
tolerance = eps(T), debug = false, maxlevel = 30)
59+
function find_roots(f::Function, a::Interval{T}, method::Function = newton;
60+
tolerance = eps(T), debug = false, maxlevel = 30) where {T}
7461

7562
method(f, a; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
7663
end
7764

78-
function find_roots{T}(f::Function, f_prime::Function, a::Interval{T}, method::Function=newton;
79-
tolerance=eps(T), debug=false, maxlevel=30)
65+
function find_roots(f::Function, f_prime::Function, a::Interval{T}, method::Function=newton;
66+
tolerance=eps(T), debug=false, maxlevel=30) where {T}
8067

8168
method(f, f_prime, a; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
8269
end
@@ -127,7 +114,7 @@ function find_roots_midpoint(f::Function, a::Real, b::Real, method::Function=new
127114

128115
end
129116

130-
function Base.lexcmp{T}(a::Interval{T}, b::Interval{T})
117+
function Base.lexcmp(a::Interval{T}, b::Interval{T}) where {T}
131118
#@show a, b
132119
if a.lo < b.lo
133120
return -1
@@ -142,7 +129,7 @@ function Base.lexcmp{T}(a::Interval{T}, b::Interval{T})
142129

143130
end
144131

145-
Base.lexcmp{T}(a::Root{T}, b::Root{T}) = lexcmp(a.interval, b.interval)
132+
Base.lexcmp(a::Root{T}, b::Root{T}) where {T} = lexcmp(a.interval, b.interval)
146133

147134

148135
function clean_roots(f, roots)

‎src/bisect.jl

Lines changed: 0 additions & 42 deletions
This file was deleted.

‎src/complex.jl

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -6,22 +6,22 @@ struct Compl{T} <: Number
66
im::T
77
end
88

9-
Base.promote_rule{T, S<:Real}(::Type{Compl{T}}, ::Type{S}) = Compl{T}
10-
Base.convert{T}(::Type{Compl{T}}, x::Real) = Compl{T}(x, 0)
9+
Base.promote_rule(::Type{Compl{T}}, ::Type{S}) where {T, S<:Real} = Compl{T}
10+
Base.convert(::Type{Compl{T}}, x::Real) where {T} = Compl{T}(x, 0)
1111

1212
Base.show(io::IO, c::Compl) = println(io, c.re, " + ", c.im, "im")
1313

14-
Base.:+{T}(a::Compl{T}, b::Compl{T}) = Compl{T}(a.re+b.re, a.im+b.im)
15-
Base.:-{T}(a::Compl{T}, b::Compl{T}) = Compl{T}(a.re-b.re, a.im-b.im)
14+
Base.:+(a::Compl{T}, b::Compl{T}) where {T} = Compl{T}(a.re+b.re, a.im+b.im)
15+
Base.:-(a::Compl{T}, b::Compl{T}) where {T} = Compl{T}(a.re-b.re, a.im-b.im)
1616

17-
Base.:*{T}(a::Compl{T}, b::Compl{T}) = Compl{T}(a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re)
17+
Base.:*(a::Compl{T}, b::Compl{T}) where {T} = Compl{T}(a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re)
1818

1919

20-
Base.:*{T}(α::Number, z::Compl{T}) = Compl{T}(α*z.re, α*z.im)
20+
Base.:*(α::Number, z::Compl{T}) where {T} = Compl{T}(α*z.re, α*z.im)
2121

22-
Base.one{T}(z::Compl{T}) = Compl{T}(one(T), zero(T))
22+
Base.one(z::Compl{T}) where {T} = Compl{T}(one(T), zero(T))
2323

24-
Base.copy{T}(z::Compl{T}) = Compl{T}(z.re, z.im)
24+
Base.copy(z::Compl{T}) where {T} = Compl{T}(z.re, z.im)
2525

2626
"""
2727
Takes a complex (polynomial) function f and returns a function g:R^2 -> R^2

‎src/contractors.jl

Lines changed: 9 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,5 @@
11
export Bisection, Newton, Krawczyk
22

3-
Base.isinf(X::IntervalBox) = any(isinf.(X))
4-
IntervalArithmetic.mid(X::IntervalBox, α) = mid.(X, α)
53

64
doc"""
75
Contractor{F}
@@ -81,19 +79,19 @@ end
8179
doc"""
8280
Single-variable Newton operator
8381
"""
84-
function 𝒩{T}(f, X::Interval{T})
82+
function 𝒩(f, X::Interval{T}) where {T}
8583
m = Interval(mid(X, where_bisect))
8684

8785
m - (f(m) / ForwardDiff.derivative(f, X))
8886
end
8987

90-
function 𝒩{T}(f, f′, X::Interval{T})
88+
function 𝒩(f, f′, X::Interval{T}) where {T}
9189
m = Interval(mid(X, where_bisect))
9290

9391
m - (f(m) / f′(X))
9492
end
9593

96-
function 𝒩{T}(f, X::Interval{T}, dX::Interval{T})
94+
function 𝒩(f, X::Interval{T}, dX::Interval{T}) where {T}
9795
m = Interval(mid(X, where_bisect))
9896

9997
m - (f(m) / dX)
@@ -104,9 +102,9 @@ Multi-variable Newton operator.
104102
"""
105103
function 𝒩(f::Function, jacobian::Function, X::IntervalBox) # multidimensional Newton operator
106104
m = IntervalBox(Interval.(mid(X, where_bisect)))
107-
J = jacobian(SVector(X))
105+
J = jacobian(X)
108106

109-
return IntervalBox(m - (J \ f(m)))
107+
return IntervalBox(m .- (J \ f(m)))
110108
end
111109

112110

@@ -134,9 +132,9 @@ Single-variable Krawczyk operator
134132
"""
135133
function 𝒦(f, f′, X::Interval{T}) where {T}
136134
m = Interval(mid(X))
137-
Y = 1/f′(m)
135+
Y = 1 / f′(m)
138136

139-
m - Y*f(m) + (1 - Y*f′(X))*(X - m)
137+
m - Y*f(m) + (1 - Y*f′(X)) * (X - m)
140138
end
141139

142140
doc"""
@@ -146,9 +144,9 @@ function 𝒦(f, jacobian, X::IntervalBox{T}) where {T}
146144
m = mid(X)
147145
J = jacobian(X)
148146
Y = inv(jacobian(m))
149-
m = IntervalBox(Interval.(m))
147+
mm = IntervalBox(m)
150148

151-
IntervalBox(m - Y*f(m) + (I - Y*J)*(X - m))
149+
res = m - Y*f(mm) + (I - Y*J) * (X.v - m) # IntervalBox(res)
152150
end
153151

154152
"""

‎src/krawczyk.jl

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
doc"""Returns two intervals, the first being a point within the
88
interval x such that the interval corresponding to the derivative of f there
99
does not contain zero, and the second is the inverse of its derivative"""
10-
function guarded_derivative_midpoint{T}(f::Function, f_prime::Function, x::Interval{T})
10+
function guarded_derivative_midpoint(f::Function, f_prime::Function, x::Interval{T}) where {T}
1111

1212
α = convert(T, 0.46875) # close to 0.5, but exactly representable as a floating point
1313

@@ -31,16 +31,16 @@ function guarded_derivative_midpoint{T}(f::Function, f_prime::Function, x::Inter
3131
end
3232

3333

34-
function K{T}(f::Function, f_prime::Function, x::Interval{T})
34+
function K(f::Function, f_prime::Function, x::Interval{T}) where {T}
3535
m, C = guarded_derivative_midpoint(f, f_prime, x)
3636
deriv = f_prime(x)
3737
Kx = m - C*f(m) + (one(T) - C*deriv) * (x - m)
3838
Kx
3939
end
4040

4141

42-
function krawczyk_refine{T}(f::Function, f_prime::Function, x::Interval{T};
43-
tolerance=eps(one(T)), debug=false)
42+
function krawczyk_refine(f::Function, f_prime::Function, x::Interval{T};
43+
tolerance=eps(one(T)), debug=false) where {T}
4444

4545
debug && (print("Entering krawczyk_refine:"); @show x)
4646

@@ -58,8 +58,8 @@ end
5858

5959

6060

61-
function krawczyk{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
62-
tolerance=eps(one(T)), debug=false, maxlevel=30)
61+
function krawczyk(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
62+
tolerance=eps(one(T)), debug=false, maxlevel=30) where {T}
6363

6464
debug && (print("Entering krawczyk:"); @show(level); @show(x))
6565

@@ -98,15 +98,15 @@ end
9898

9999

100100
# use automatic differentiation if no derivative function given
101-
krawczyk{T}(f::Function,x::Interval{T}; args...) =
101+
krawczyk(f::Function,x::Interval{T}; args...) where {T} =
102102
krawczyk(f, x->D(f,x), x; args...)
103103

104104
# krawczyk for vector of intervals:
105-
krawczyk{T}(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) =
105+
krawczyk(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) where {T} =
106106
vcat([krawczyk(f, f_prime, @interval(x); args...) for x in xx]...)
107107

108-
krawczyk{T}(f::Function, xx::Vector{Interval{T}}, level; args...) =
108+
krawczyk(f::Function, xx::Vector{Interval{T}}, level; args...) where {T} =
109109
krawczyk(f, x->D(f,x), xx; args...)
110110

111-
krawczyk{T}(f::Function, xx::Vector{Root{T}}; args...) =
111+
krawczyk(f::Function, xx::Vector{Root{T}}; args...) where {T} =
112112
krawczyk(f, x->D(f,x), [x.interval for x in xx]; args...)

‎src/newton.jl

Lines changed: 10 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,8 @@
66

77
doc"""If a root is known to be inside an interval,
88
`newton_refine` iterates the interval Newton method until that root is found."""
9-
function newton_refine{N,T}(f::Function, f_prime::Function, X::Union{Interval{T}, IntervalBox{N,T}};
10-
tolerance=eps(T), debug=false)
9+
function newton_refine(f::Function, f_prime::Function, X::Union{Interval{T}, IntervalBox{N,T}};
10+
tolerance=eps(T), debug=false) where {N,T}
1111

1212
debug && (print("Entering newton_refine:"); @show X)
1313

@@ -32,8 +32,8 @@ end
3232

3333
doc"""If a root is known to be inside an interval,
3434
`newton_refine` iterates the interval Newton method until that root is found."""
35-
function newton_refine{T}(f::Function, f_prime::Function, X::Interval{T};
36-
tolerance=eps(T), debug=false)
35+
function newton_refine(f::Function, f_prime::Function, X::Interval{T};
36+
tolerance=eps(T), debug=false) where {T}
3737

3838
debug && (print("Entering newton_refine:"); @show X)
3939

@@ -65,8 +65,8 @@ with its optional derivative `f_prime` and initial interval `x`.
6565
Optional keyword arguments give the `tolerance`, `maxlevel` at which to stop
6666
subdividing, and a `debug` boolean argument that prints out diagnostic information."""
6767

68-
function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
69-
tolerance=eps(T), debug=false, maxlevel=30)
68+
function newton(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
69+
tolerance=eps(T), debug=false, maxlevel=30) where {T}
7070

7171
debug && (print("Entering newton:"); @show(level); @show(x))
7272

@@ -138,20 +138,19 @@ function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
138138
roots = clean_roots(f, roots)
139139

140140
return roots
141-
142141
end
143142

144143

145144
# use automatic differentiation if no derivative function given:
146-
newton{T}(f::Function, x::Interval{T}; args...) =
145+
newton(f::Function, x::Interval{T}; args...) where {T} =
147146
newton(f, x->D(f,x), x; args...)
148147

149148
# newton for vector of intervals:
150-
newton{T}(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) =
149+
newton(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) where {T} =
151150
vcat([newton(f, f_prime, @interval(x); args...) for x in xx]...)
152151

153-
newton{T}(f::Function, xx::Vector{Interval{T}}, level; args...) =
152+
newton(f::Function, xx::Vector{Interval{T}}, level; args...) where {T} =
154153
newton(f, x->D(f,x), xx, 0, args...)
155154

156-
newton{T}(f::Function, xx::Vector{Root{T}}; args...) =
155+
newton(f::Function, xx::Vector{Root{T}}; args...) where {T} =
157156
newton(f, x->D(f,x), [x.interval for x in xx], args...)

‎src/roots.jl

Lines changed: 5 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -7,10 +7,6 @@ export start, next, done, copy, step!, eltype, iteratorsize
77

88
diam(x::Root) = diam(x.interval)
99

10-
Base.size(x::Interval) = (1,)
11-
12-
isinterior{N}(X::IntervalBox{N}, Y::IntervalBox{N}) = all(isinterior.(X, Y))
13-
1410
struct RootSearchState{T <: Union{Interval,IntervalBox}}
1511
working::Vector{T}
1612
outputs::Vector{Root{T}}
@@ -110,13 +106,8 @@ function recursively_branch_and_prune(h, X, contractor=BisectionContractor, fina
110106
end
111107

112108

113-
contains_zero(X::Interval{T}) where {T} = zero(T) ∈ X
114-
contains_zero(X::SVector) = all(contains_zero.(X))
115-
contains_zero(X::IntervalBox) = all(contains_zero.(X))
116-
117-
118-
IntervalLike{T} = Union{Interval{T}, IntervalBox{T}}
119-
NewtonLike = Union{Type{Newton}, Type{Krawczyk}}
109+
const IntervalLike{T} = Union{Interval{T}, IntervalBox{T}}
110+
const NewtonLike = Union{Type{Newton}, Type{Krawczyk}}
120111

121112
"""
122113
roots(f, X, contractor, tol=1e-3)
@@ -183,7 +174,7 @@ function roots(f, Xc::Complex{Interval{T}}, contractor::Type{C},
183174
tol::Float64=1e-3) where {T, C<:Contractor}
184175

185176
g = realify(f)
186-
Y = IntervalBox(reim(Xc))
177+
Y = IntervalBox(reim(Xc)...)
187178
rts = roots(g, Y, contractor, tol)
188179

189180
return [Root(Complex(root.interval...), root.status) for root in rts]
@@ -195,7 +186,7 @@ function roots(f, Xc::Complex{Interval{T}}, C::NewtonLike, tol::Float64=1e-3) wh
195186

196187
g_prime = x -> ForwardDiff.jacobian(g, x)
197188

198-
Y = IntervalBox(reim(Xc))
189+
Y = IntervalBox(reim(Xc)...)
199190
rts = roots(g, g_prime, Y, C, tol)
200191

201192
return [Root(Complex(root.interval...), root.status) for root in rts]
@@ -207,7 +198,7 @@ function roots(f, deriv, Xc::Complex{Interval{T}}, C::NewtonLike, tol::Float64=1
207198

208199
g_prime = realify_derivative(deriv)
209200

210-
Y = IntervalBox(reim(Xc))
201+
Y = IntervalBox(reim(Xc)...)
211202
rts = roots(g, g_prime, Y, C, tol)
212203

213204
return [Root(Complex(root.interval...), root.status) for root in rts]

0 commit comments

Comments
 (0)