Skip to content

Commit ae55541

Browse files
eeshan9815eeshan9815
committed
Add icp contractors
1 parent a9c1eab commit ae55541

File tree

3 files changed

+195
-66
lines changed

3 files changed

+195
-66
lines changed

‎examples/optimise1.jl

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
using IntervalOptimisation, IntervalArithmetic, IntervalConstraintProgramming
2+
3+
4+
# Unconstrained Optimisation
5+
f(x, y) = (1.5 - x * (1 - y))^2 + (2.25 - x * (1 - y^2))^2 + (2.625 - x * (1 - y^3))^2
6+
# f (generic function with 2 methods)
7+
8+
f(X) = f(X...)
9+
# f (generic function with 2 methods)
10+
11+
X = (-1e6..1e6) × (-1e6..1e6)
12+
# [-1e+06, 1e+06] × [-1e+06, 1e+06]
13+
14+
minimise_icp(f, X)
15+
# ([0, 2.39527e-07], IntervalArithmetic.IntervalBox{2,Float64}[[2.99659, 2.99748] × [0.498906, 0.499523], [2.99747, 2.99834] × [0.499198, 0.500185], [2.99833, 2.99922] × [0.499198, 0.500185], [2.99921, 3.00011] × [0.499621, 0.500415], [3.0001, 3.001] × [0.499621, 0.500415], [3.00099, 3.0017] × [0.500169, 0.500566], [3.00169, 3.00242] × [0.500169, 0.500566]])
16+
17+
18+
19+
f(X) = X[1]^2 + X[2]^2
20+
# f (generic function with 2 methods)
21+
22+
X = (-∞..∞) × (-∞..∞)
23+
# [-∞, ∞] × [-∞, ∞]
24+
25+
minimise_icp(f, X)
26+
# ([0, 0], IntervalArithmetic.IntervalBox{2,Float64}[[-0, 0] × [-0, 0], [0, 0] × [-0, 0]])
27+
28+
29+
30+
31+
# Constrained Optimisation
32+
f(X) = -1 * (X[1] + X[2])
33+
# f (generic function with 1 method)
34+
35+
X = (-∞..∞) × (-∞..∞)
36+
# [-∞, ∞] × [-∞, ∞]
37+
38+
constraints = [IntervalOptimisation.constraint(x->(x[1]), -∞..6), IntervalOptimisation.constraint(x->x[2], -∞..4), IntervalOptimisation.constraint(x->(x[1]*x[2]), -∞..4)]
39+
# 3-element Array{IntervalOptimisation.constraint{Float64},1}:
40+
# IntervalOptimisation.constraint{Float64}(#3, [-∞, 6])
41+
# IntervalOptimisation.constraint{Float64}(#4, [-∞, 4])
42+
# IntervalOptimisation.constraint{Float64}(#5, [-∞, 4])
43+
44+
minimise_icp_constrained(f, X, constraints)
45+
# ([-6.66676, -6.66541], IntervalArithmetic.IntervalBox{2,Float64}[[5.99918, 6] × [0.666233, 0.666758], [5.99918, 6] × [0.665717, 0.666234], [5.99887, 5.99919] × [0.666233, 0.666826], [5.99984, 6] × [0.665415, 0.665718], [5.99856, 5.99888] × [0.666233, 0.666826], [5.99969, 5.99985] × [0.665415, 0.665718]])
46+
47+
48+
49+
# One-dimensional case
50+
minimise1d(x -> (x^2 - 2)^2, -10..11)
51+
# ([0, 1.33476e-08], IntervalArithmetic.Interval{Float64}[[-1.41426, -1.41358], [1.41364, 1.41429]])
52+
53+
minimise1d_deriv(x -> (x^2 - 2)^2, -10..11)
54+
# ([0, 8.76812e-08], IntervalArithmetic.Interval{Float64}[[-1.41471, -1.41393], [1.41367, 1.41444]])

‎src/IntervalOptimisation.jl

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,12 +3,14 @@
33
module IntervalOptimisation
44

55
export minimise, maximise,
6-
minimize, maximize
6+
minimize, maximize,
7+
minimise1d, minimise1d_deriv,
8+
minimise_icp, minimise_icp_constrained
79

810
include("SortedVectors.jl")
911
using .SortedVectors
1012

11-
using IntervalArithmetic, IntervalRootFinding
13+
using IntervalArithmetic, IntervalRootFinding, DataStructures, IntervalConstraintProgramming, StaticArrays, ForwardDiff
1214

1315
if !isdefined(:sup)
1416
const sup = supremum
@@ -19,8 +21,9 @@ if !isdefined(:inf)
1921
end
2022

2123
include("optimise.jl")
24+
include("optimise1.jl")
2225

2326
const minimize = minimise
24-
const maximize = maximise
27+
const maximize = maximise
2528

2629
end

‎src/optimise1.jl

Lines changed: 135 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,25 @@
1-
using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat, IntervalOptimisation, DataStructures
1+
using IntervalRootFinding, IntervalArithmetic, StaticArrays, ForwardDiff, BenchmarkTools, Compat, IntervalOptimisation, DataStructures, IntervalConstraintProgramming
22

33
import Base.isless
44

55
struct IntervalMinima{T<:Real}
66
intval::Interval{T}
7-
minima::T
7+
global_minimum::T
88
end
99

1010
function isless(a::IntervalMinima{T}, b::IntervalMinima{T}) where {T<:Real}
11-
return isless(a.minima, b.minima)
11+
return isless(a.global_minimum, b.global_minimum)
1212
end
1313

1414
function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real}
1515

1616
Q = binary_minheap(IntervalMinima{T})
1717

18-
minima = f(interval(mid(x))).hi
18+
global_minimum = f(interval(mid(x))).hi
1919

2020
arg_minima = Interval{T}[]
2121

22-
push!(Q, IntervalMinima(x, minima))
22+
push!(Q, IntervalMinima(x, global_minimum))
2323

2424
while !isempty(Q)
2525

@@ -29,51 +29,90 @@ function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where
2929
continue
3030
end
3131

32-
if p.minima > minima
32+
if p.global_minimum > global_minimum
3333
continue
3434
end
3535

36-
# deriv = ForwardDiff.derivative(f, p.intval)
37-
# if 0 ∉ deriv
38-
# continue
39-
# end
36+
current_minimum = f(interval(mid(p.intval))).hi
37+
38+
if current_minimum < global_minimum
39+
global_minimum = current_minimum
40+
end
41+
42+
if diam(p.intval) < abstol
43+
push!(arg_minima, p.intval)
44+
else
45+
x1, x2 = bisect(p.intval)
46+
push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo))
47+
end
48+
end
49+
lb = minimum(inf.(f.(arg_minima)))
50+
51+
return lb..global_minimum, arg_minima
52+
end
53+
54+
function minimise1d_deriv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real}
55+
56+
Q = binary_minheap(IntervalMinima{T})
57+
58+
global_minimum = f(interval(mid(x))).hi
59+
60+
arg_minima = Interval{T}[]
61+
62+
push!(Q, IntervalMinima(x, global_minimum))
63+
64+
while !isempty(Q)
4065

41-
doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.intval)
42-
if doublederiv < 0
66+
p = pop!(Q)
67+
68+
if isempty(p.intval)
69+
continue
70+
end
71+
72+
if p.global_minimum > global_minimum
4373
continue
4474
end
4575

76+
deriv = ForwardDiff.derivative(f, p.intval)
77+
if 0 ∉ deriv
78+
continue
79+
end
80+
# Second derivative contractor
81+
# doublederiv = ForwardDiff.derivative(x->ForwardDiff.derivative(f, x), p.intval)
82+
# if doublederiv < 0
83+
# continue
84+
# end
85+
4686
m = mid(p.intval)
4787

48-
current_minima = f(interval(m)).hi
88+
current_minimum = f(interval(m)).hi
4989

50-
if current_minima < minima
51-
minima = current_minima
90+
if current_minimum < global_minimum
91+
global_minimum = current_minimum
5292
end
5393

5494

55-
# # Contractor 1
56-
# x = m .+ extended_div((interval(-∞, minima) - f(m)), deriv)
57-
#
58-
# x = x .∩ p.intval
59-
60-
# Contractor 2
61-
@show p.intval
62-
@show (doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima))
63-
@show (m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima)))
64-
x = m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, minima))
95+
# Contractor 1
96+
x = m .+ extended_div((interval(-∞, global_minimum) - f(m)), deriv)
6597

6698
x = x .∩ p.intval
6799

100+
# # Contractor 2 (Second derivative, expanding more on the Taylor Series, not beneficial in practice)
101+
# x = m .+ quadratic_roots(doublederiv/2, ForwardDiff.derivative(f, interval(m)), f(m) - interval(-∞, global_minimum))
102+
#
103+
# x = x .∩ p.intval
104+
68105
if diam(p.intval) < abstol
69106
push!(arg_minima, p.intval)
70107
else
71-
if length(x) == 1
108+
if isempty(x[2])
72109
x1, x2 = bisect(x[1])
73110
push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo))
74111
else
75112
push!(Q, IntervalMinima.(x, inf.(f.(x)))...)
76113
end
114+
115+
# Second Deriv contractor
77116
# if isempty(x[2])
78117
# x1, x2 = bisect(x[1])
79118
# push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo))
@@ -83,94 +122,127 @@ function minimise1d(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where
83122
# end
84123
end
85124
end
125+
lb = minimum(inf.(f.(arg_minima)))
86126

87-
return minima, arg_minima
127+
return lb..global_minimum, arg_minima
88128
end
89129

90130

91-
function minimise1d_noderiv(f::Function, x::Interval{T}; reltol=1e-3, abstol=1e-3) where {T<:Real}
92131

93-
Q = binary_minheap(IntervalMinima{T})
94132

95-
minima = f(interval(mid(x))).hi
96-
arg_minima = Interval{T}[]
133+
struct IntervalBoxMinima{N, T<:Real}
134+
intval::IntervalBox{N, T}
135+
global_minimum::T
136+
end
137+
138+
struct constraint{T<:Real}
139+
f::Function
140+
c::Interval{T}
141+
end
142+
143+
function isless(a::IntervalBoxMinima{N, T}, b::IntervalBoxMinima{N, T}) where {N, T<:Real}
144+
return isless(a.global_minimum, b.global_minimum)
145+
end
146+
147+
function minimise_icp(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e-3) where {N, T<:Real}
148+
149+
Q = binary_minheap(IntervalBoxMinima{N, T})
150+
151+
global_minimum = ∞
152+
153+
x = icp(f, x, -∞..global_minimum)
97154

98-
push!(Q, IntervalMinima(x, minima))
155+
arg_minima = IntervalBox{N, T}[]
156+
157+
push!(Q, IntervalBoxMinima(x, global_minimum))
99158

100159
while !isempty(Q)
101160

102161
p = pop!(Q)
103162

104-
if p.minima > minima
163+
if isempty(p.intval)
105164
continue
106165
end
107-
#
108-
# if 0 ∉ ForwardDiff.derivative(f, p.intval)
109-
# continue
110-
# end
111166

112-
# current_minima = f(p.intval).lo
167+
if p.global_minimum > global_minimum
168+
continue
169+
end
113170

114-
current_minima = f(interval(mid(p.intval))).hi
171+
current_minimum = f(interval.(mid(p.intval))).hi
115172

116-
if current_minima < minima
117-
minima = current_minima
173+
if current_minimum < global_minimum
174+
global_minimum = current_minimum
118175
end
176+
# if all(0 .∉ ForwardDiff.gradient(f, p.intval.v))
177+
# continue
178+
# end
179+
180+
X = icp(f, p.intval, -∞..global_minimum)
119181

120182
if diam(p.intval) < abstol
121183
push!(arg_minima, p.intval)
184+
122185
else
123-
x1, x2 = bisect(p.intval)
124-
push!(Q, IntervalMinima(x1, f(x1).lo), IntervalMinima(x2, f(x2).lo))
186+
x1, x2 = bisect(X)
187+
push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo))
125188
end
126189
end
127190

128-
return minima, arg_minima
129-
end
191+
lb = minimum(inf.(f.(arg_minima)))
130192

193+
return lb..global_minimum, arg_minima
194+
end
131195

196+
function minimise_icp_constrained(f::Function, x::IntervalBox{N, T}, constraints::Vector{constraint{T}} = Vector{constraint{T}}(); reltol=1e-3, abstol=1e-3) where {N, T<:Real}
132197

133-
struct IntervalBoxMinima{N, T<:Real}
134-
intval::IntervalBox{N, T}
135-
minima::T
136-
end
198+
Q = binary_minheap(IntervalBoxMinima{N, T})
137199

138-
function isless(a::IntervalBoxMinima{N, T}, b::IntervalBoxMinima{N, T}) where {N, T<:Real}
139-
return isless(a.minima, b.minima)
140-
end
200+
global_minimum = ∞
141201

142-
function minimisend(f::Function, x::IntervalBox{N, T}; reltol=1e-3, abstol=1e-3) where {N, T<:Real}
202+
for t in constraints
203+
x = icp(t.f, x, t.c)
204+
end
143205

144-
Q = binary_minheap(IntervalBoxMinima{N, T})
206+
x = icp(f, x, -∞..global_minimum)
145207

146-
minima = f(x).lo
147208
arg_minima = IntervalBox{N, T}[]
148209

149-
push!(Q, IntervalBoxMinima(x, minima))
210+
push!(Q, IntervalBoxMinima(x, global_minimum))
150211

151212
while !isempty(Q)
152213

153214
p = pop!(Q)
154215

155-
if p.minima > minima
216+
if isempty(p.intval)
156217
continue
157218
end
158219

159-
if 0 .∉ ForwardDiff.gradient(f, p.intval)
220+
if p.global_minimum > global_minimum
160221
continue
161222
end
162223

163-
if p.minima < minima
164-
minima = p.minima
224+
# current_minimum = f(interval.(mid(p.intval))).hi
225+
current_minimum = f(p.intval).hi
226+
227+
if current_minimum < global_minimum
228+
global_minimum = current_minimum
229+
end
230+
# if 0 .∉ ForwardDiff.gradient(f, p.intval.v)
231+
# continue
232+
# end
233+
X = icp(f, p.intval, -∞..global_minimum)
234+
235+
for t in constraints
236+
X = icp(t.f, X, t.c)
165237
end
166238

167239
if diam(p.intval) < abstol
168240
push!(arg_minima, p.intval)
169241
else
170-
x1, x2 = bisect(p.intval)
242+
x1, x2 = bisect(X)
171243
push!(Q, IntervalBoxMinima(x1, f(x1).lo), IntervalBoxMinima(x2, f(x2).lo))
172244
end
173245
end
174-
175-
return minima, arg_minima
246+
lb = minimum(inf.(f.(arg_minima)))
247+
return lb..global_minimum, arg_minima
176248
end

0 commit comments

Comments
 (0)