juliaintervals / intervaloptimisation.jl Goto Github PK
View Code? Open in Web Editor NEWRigorous global optimisation in pure Julia
License: Other
Rigorous global optimisation in pure Julia
License: Other
julia> using IntervalOptimisation, IntervalArithmetic
julia> f(x) = 3x^2 - 1
f (generic function with 1 method)
julia> dom = Interval(2..6)
[2, 6]
julia> minimise(f, dom)
([11, 11.0054], Interval{Float64}[[2, 2.00089]])
julia> maximise(f, dom)
ERROR: MethodError: no method matching minimise(::getfield(IntervalOptimisation, Symbol("##5#6")){typeof(f)}, ::Interval{Float64}, ::Type{HeapedV
ector}, ::Float64)
Closest candidates are:
minimise(::Any, ::T; structure, tol) where T at /Users/forets/.julia/dev/IntervalOptimisation/src/optimise.jl:21
Stacktrace:
[1] #maximise#4(::Type, ::Float64, ::Function, ::typeof(f), ::Interval{Float64}) at /Users/forets/.julia/dev/IntervalOptimisation/src/optimise.j
l:73
[2] maximise(::Function, ::Interval{Float64}) at /Users/forets/.julia/dev/IntervalOptimisation/src/optimise.jl:73
[3] top-level scope at none:0
In the following problem, the actual minimal points are simply (0,1) and (1,0). but the package outputs a lot of boxes. It does not seem that decreasing tolerance helps.
julia> global_min, minimisers = minimise(X->((x,y)=X; -(x/(y+1)+y/(x+1))), IntervalBox(0..1, 2), tol=1e-15);
julia> minimisers
20-element Array{IntervalBox{2,Float64},1}:
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0, 6.00051e-16] × [0.999999, 1]
[6.0005e-16, 1.20956e-15] × [0.999999, 1]
[0.999999, 1] × [0, 6.00051e-16]
[0.999999, 1] × [6.0005e-16, 1.20956e-15]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [1.20955e-15, 1.81906e-15]
[1.20955e-15, 1.81906e-15] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
[0.999999, 1] × [0.999999, 1]
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
It seems to me like IntervalOptimization would be well suited optimization problems where the variables are constrained to be integers.
AIUI, nothing fundermental about the algorithm means it should not work on integers.
And Integer Optimization is already very expensive so the fact that Interval Optimization is expensive is less of a problem.
one problem is the package is currently hardcoded to use mid
and bisect
with the default α
to perform the bisection.
The choice of α
doesn't AIUI effect the correctness of the algorithm.
But the current method does give non-integer values.
Not sure if this is a real issue.
Is a HeapedVector meant to be sorted? I gathered so from the code, specifically the call to floatup!(). If so, then there is a bug in floatup!(). If not, then please disregard the rest of this.
julia> v = collect(1:2:19); insert!(v, length(v)+1, 4)
11-element Array{Int64,1}:
1
3
5
7
9
11
13
15
17
19
4
julia> floatup!(v, length(v), identity)
11-element Array{Int64,1}:
1
3
5
7
4
11
13
15
17
19
9
The example is based on my understanding that when a HeapedVector is created, the data is added one value at a time (obviously intervals, but this is easier to demo with floats) and floatup!() will then sort the newly added value into the existing set.
The error (if this is an error - see above) stems from the line
par = convert(Int, floor(index/2))
I would propose a correction, which works, but is not especially efficient:
function floatup2!(ar, index, by)
if index <= 1
return ar
else
par = index - 1
end
if by(ar[index]) < by(ar[par])
ar[par], ar[index] = ar[index], ar[par]
end
floatup2!(ar, par, by)
end
This sorts as expected:
julia> v = collect(1:2:19); insert!(v, length(v)+1, 4)
11-element Array{Int64,1}:
1
3
5
7
9
11
13
15
17
19
4
julia> floatup2!(v, length(v), identity)
11-element Array{Int64,1}:
1
3
4
5
7
9
11
13
15
17
19
It would be useful to have an OptimizationProblem
type that encodes the function and original box, dimension etc.
Cf. the one in IntervalRootFinding.jl
Sometimes evaluating the error function is quite expensive, so it's good if one can provide a maximum number of function evaluation as an option to minimise
. This would also allow to benchmark the method for a fixed function call budget.
@dpsanders I can try to do a PR if you give me some pointers.
Maybe in IntervalRootFinding.jl.
Add Goldstein-Price test case:
#11
Hi, as a new user, I wasn't sure how to create the x cross
symbol that can be used in place of IntervalBox
. It would be nice to specify the \tab command to create it in the docs. Thanks.
While working on a project I found that minimise
throwed an error because I didn't find a minimum. The MWE is the next
julia> using IntervalArithmetic, IntervalOptimisation
julia> cost(s, a) = 1 + (s^2 + a^2) * abs(s-a) + abs(a)
cost (generic function with 2 methods)
julia> x = -1..1
[-1, 1]
julia> minimise(X -> cost(X[1], X[2]), x×x)
ArgumentError: reducing over an empty collection is not allowed
IntervalOptimisation.jl
version is v0.4.3
, versioninfo()
is
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
JULIA_NUM_THREADS = 4
minimise struggles with this.
https://github.com/JuliaOpt/SumOfSquares.jl/blob/master/examples/goldsteinprice.jl
Using julia 1.6 on windows, the test Discontinuous function in 1D
fails with the following stack trace
julia> H(x) = (sign(x) + 1) / 2 # heaviside function
H (generic function with 1 method)
julia> global_min, minimisers = minimise(x -> abs(x) + H(x) - 1, -10..11, tol=1e-5)
ERROR: could not load symbol "fesetround":
The specified procedure could not be found.
Stacktrace:
[1] setrounding_raw
@ ~\.julia\packages\SetRounding\MjtxK\src\SetRounding.jl:45 [inlined]
[2] setrounding(#unused#::Type{Float64}, r::RoundingMode{:Down})
@ SetRounding ~\.julia\packages\SetRounding\MjtxK\src\SetRounding.jl:47
[3] setrounding(f::IntervalArithmetic.var"#127#128"{Interval{Float64}}, #unused#::Type{Float64}, rounding::RoundingMode{:Down})
@ Base.Rounding .\rounding.jl:174
[4] mig
@ ~\.julia\packages\IntervalArithmetic\odPjL\src\intervals\arithmetic.jl:297 [inlined]
[5] abs
@ ~\.julia\packages\IntervalArithmetic\odPjL\src\intervals\arithmetic.jl:314 [inlined]
[6] (::var"#7#8")(x::Interval{Float64})
@ Main .\REPL[6]:1
[7] minimise(f::var"#7#8", X::Interval{Float64}; structure::Type{HeapedVector}, tol::Float64)
@ IntervalOptimisation ~\.julia\packages\IntervalOptimisation\PjZLd\src\optimise.jl:40
[8] top-level scope
@ REPL[6]:1
with julia 1.5 (on windows) it works fine.
julia> minimise(x->(x^2 - 2)^2, big(-10..10))
ERROR: MethodError: no method matching filter_elements!(::HeapedVector{Tuple{Interval{BigFloat},Float64},IntervalOptimisation.var"#2#3"}, ::Tuple{Interval{BigFloat},BigFloat})
Hi, I am currently using IntervalOptimisation
to find the global minimum (or at least a tight bound) for a nonlinear least-squares problem that includes 5 parameters to be optimized.
The dataset in a CSV file that contains 26 data points can be downloaded here.
Julia code
using CSV, Tables
using IntervalArithmetic, IntervalOptimisation
# read data into a 26×2 Array
const data = CSV.File("data.csv") |> Tables.matrix
const Vt = 0.02638199348809556
function sse_sdm(θ)
Ir, I0, n, Rs, Rp = θ
return sum((Ir - I0*1e-6*(exp((V + I*Rs)/(n*Vt)) - 1) - (V + I*Rs)/Rp - I)^2 for (V, I) in eachrow(data))
end
intervals = (0..1)×(0..1)×(1..2)×(0..0.5)×(0..100)
@time minimum, minimisers = minimise(sse_sdm, intervals)
The above code has run more than one and a half hours (still running now). I am wondering whether there are any suggestions that can help accelerate the optimization. I saw in the Julia discourse that the nonrecursive_powers branch may be faster; however, I cannot find it.
I am not quite sure what's wrong with it: the url looks correct.
I'm having some trouble minimizing a function that includes complex numbers and broadcasting. I'm not sure which is causing the problem. On slack @dpsanders mentioned being a little concerned by complex numbers, but looking at the stack trace I'm actually wondering if its a broadcasting problem (it seems to be trying to assign an interval into the result array from a broadcasting operation, and for some reason the result array has float elements.
I'm coming from a DSP background, and the application is delay estimation between two signals where the delay is not an integer number of samples. The main insight here is that a delay in the time domain is the same as multiplying by a complex exponential in the frequency domain. The frequency of the exponential determines the delay, so we can estimate the delay by finding the frequency that minimizes the distance between the two spectra. Here's a MWE:
using Plots
using IntervalOptimisation
using IntervalArithmetic
using FFTW
# first we generate some test data
D = 0.68 # just picking a random target delay (in samples)
x1 = randn(1000)
N = length(x1) ÷ 2 + 1
ω = range(0, π, length=N)
# create a target signal by delaying our original signal and adding some noise
x2 = irfft(rfft(x1) .* exp.(-im .* ω .* D), length(x1)) .+ randn.() .* 0.1
plot([x1[1:40] x2[1:40]])
We're trying to estimate D
. Here's our error function (which we see has a minimum at 0.68):
X1 = rfft(x1)
X2 = rfft(x2)
err(D) = sum(@.(abs2(X2 - X1 * exp(-im*ω*D))))
x = range(-5,5,length=100)
plot(x, err.(x))
Trying to minimise it gives:
globmin, minxs = IntervalOptimisation.minimise(err, -2..2)
MethodError: no method matching Float64(::Interval{Float64})
Closest candidates are:
Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194
Float64(::T<:Number) where T<:Number at boot.jl:741
Float64(!Matched::Int8) at float.jl:60
...
Stacktrace:
[1] convert(::Type{Float64}, ::Interval{Float64}) at ./number.jl:7
[2] setindex!(::Array{Float64,1}, ::Interval{Float64}, ::Int64) at ./array.jl:767
[3] macro expansion at ./broadcast.jl:843 [inlined]
[4] macro expansion at ./simdloop.jl:73 [inlined]
[5] copyto! at ./broadcast.jl:842 [inlined]
[6] copyto! at ./broadcast.jl:797 [inlined]
[7] copy at ./broadcast.jl:773 [inlined]
[8] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(abs2),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(-),Tuple{Array{Complex{Float64},1},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(*),Tuple{Array{Complex{Float64},1},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(exp),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(*),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(-),Tuple{Complex{Bool}}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Interval{Float64}}}}}}}}}}}) at ./broadcast.jl:753
[9] err(::Interval{Float64}) at ./In[54]:3
[10] #minimise#1(::Type, ::Float64, ::Function, ::typeof(err), ::Interval{Float64}) at /home/sfr/.julia/packages/IntervalOptimisation/A3NrJ/src/optimise.jl:34
[11] minimise(::Function, ::Interval{Float64}) at /home/sfr/.julia/packages/IntervalOptimisation/A3NrJ/src/optimise.jl:19
[12] top-level scope at In[55]:1
There's an implementation using LightGraphs in #5
It doesn't seem so easy to write an implementation by hand. It's basically a union-find algorithm
It would be nice to have a minmax
function that computes both the minimum and the maximum function.
A start could be
function minmax(f::Function, X; kwargs...)
global_min, minimiser = minimise(f, X; kwargs...)
global_max, maximiser = maximise(f, X; kwargs...)
return minimiser, global_min, maximiser, global_max
end
on the longer run, it could be nice to try to be smarted and solve both simultaneously and more efficiently.
@dpsanders what do you think?
When I tried, I got
julia> maximise(X->((x,y)=X; x/(y+1)+y/(x+1)), IntervalBox(0..1, 2), tol=1e-5)
ERROR: MethodError: no method matching minimise(::IntervalOptimisation.var"#5#6"{var"#31#32"}, ::IntervalBox{2,Float64}, ::Type{HeapedVector}, ::Float64)
Closest candidates are:
minimise(::Any, ::T; structure, tol) where T at /home/xinca341/.julia/packages/IntervalOptimisation/A3NrJ/src/optimise.jl:19
Stacktrace:
[1] #maximise#4(::Type, ::Float64, ::typeof(maximise), ::var"#31#32", ::IntervalBox{2,Float64}) at /home/xinca341/.julia/packages/IntervalOptimisation/A3NrJ/src/optimise.jl:61
[2] (::IntervalOptimisation.var"#kw##maximise")(::NamedTuple{(:tol,),Tuple{Float64}}, ::typeof(maximise), ::Function, ::IntervalBox{2,Float64}) at ./none:0
[3] top-level scope at REPL[21]:1
I am using IntervalOptimisation v0.4.0
Hi! Let me start by saying thanks for this great package! I tried the following discontinuous function optimisation where the optimal solution is on the edge of the discontinuous function. However, the output of minimise
seems wrong:
using IntervalArithmetic, IntervalOptimisation
obj = x -> log(x) + (x >= 2 ? 0 : Inf)
minimise(obj, 1..10)
# ([1.69833, 1.69839], Interval{Float64}[[5.46484, 5.46535]])
Any idea if there is a way to fix this? If this works, constraints can be trivially modeled with this discontinuous approach.
y = Interval(-0.1,0.1)
k(x) = ( ((x+y)- 2)^6 + 0.2 ) * (log(1+(x+y)^2))
global_min, minimisers = minimise(x -> sup(k(x)), -0.2..3)
results in:
ERROR: MethodError: no method matching sup(::Float64)
Closest candidates are:
sup(::Interval) at C:\Users\freemint\.julia\packages\IntervalArithmetic\UR6Qe\src\intervals\arithmetic.jl:306
sup(::DecoratedInterval{T}) where T at C:\Users\freemint\.julia\packages\IntervalArithmetic\UR6Qe\src\decorations\functions.jl:54
Stacktrace:
[1] minimise(f::var"#1#2", X::Interval{Float64}; structure::Type{HeapedVector}, tol::Float64)
@ IntervalOptimisation C:\Users\freemint\.julia\packages\IntervalOptimisation\deTUP\src\optimise.jl:40
[2] minimise(f::Function, X::Interval{Float64})
@ IntervalOptimisation C:\Users\freemint\.julia\packages\IntervalOptimisation\deTUP\src\optimise.jl:22
[3] top-level scope
@ REPL[5]:1
However:
julia> sup(x::Real) = x
sup (generic function with 3 methods)
julia> inf(x::Real) = x
inf (generic function with 3 methods)
does not fix it and raises a new issue which suggests that there might be something else going on.
ERROR: ArgumentError: reducing over an empty collection is not allowed
Stacktrace:
[1] _empty_reduce_error()
@ Base .\reduce.jl:299
[2] reduce_empty(op::Function, #unused#::Type{Real})
@ Base .\reduce.jl:309
[3] mapreduce_empty(#unused#::typeof(identity), op::Function, T::Type)
@ Base .\reduce.jl:343
[4] reduce_empty(op::Base.MappingRF{typeof(identity), typeof(min)}, #unused#::Type{Real})
@ Base .\reduce.jl:329
[5] reduce_empty_iter
@ .\reduce.jl:355 [inlined]
[6] mapreduce_empty_iter(f::Function, op::Function, itr::Vector{Real}, ItrEltype::Base.HasEltype)
@ Base .\reduce.jl:351
[7] _mapreduce(f::typeof(identity), op::typeof(min), #unused#::IndexLinear, A::Vector{Real})
@ Base .\reduce.jl:400
[8] _mapreduce_dim
@ .\reducedim.jl:318 [inlined]
[9] #mapreduce#672
@ .\reducedim.jl:310 [inlined]
[10] mapreduce
@ .\reducedim.jl:310 [inlined]
[11] #_minimum#694
@ .\reducedim.jl:878 [inlined]
[12] _minimum
@ .\reducedim.jl:878 [inlined]
[13] #_minimum#693
@ .\reducedim.jl:877 [inlined]
[14] _minimum
@ .\reducedim.jl:877 [inlined]
[15] #minimum#691
@ .\reducedim.jl:873 [inlined]
[16] minimum(a::Vector{Real})
@ Base .\reducedim.jl:873
[17] minimise(f::var"#3#4", X::Interval{Float64}; structure::Type{HeapedVector}, tol::Float64)
@ IntervalOptimisation C:\Users\freemint\.julia\packages\IntervalOptimisation\deTUP\src\optimise.jl:60
[18] minimise(f::Function, X::Interval{Float64})
@ IntervalOptimisation C:\Users\freemint\.julia\packages\IntervalOptimisation\deTUP\src\optimise.jl:22
[19] top-level scope
@ REPL[11]:1
For us uncivilized folk in the north, minimise
and maximise
might prove to be confusing
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.