juliaintervals / intervaloptimisation.jl Goto Github PK

View Code? Open in Web Editor NEW

54.0 9.0 21.0 59 KB

Rigorous global optimisation in pure Julia

License: Other

Julia 100.00%

global-optimization interval-arithmetic julia

intervaloptimisation.jl's People

Contributors

Stargazers

Watchers

intervaloptimisation.jl's Issues

MethodError in maximise

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

Add function for finding all stationary points

Problem with boundaries

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]

TagBot trigger issue

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.

HeapedVector sorting?

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.

Example:

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))

Proposal:

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

OptimizationProblem type

It would be useful to have an OptimizationProblem type that encodes the function and original box, dimension etc.

Use some kind of heap / tree instead of sorted vector

Add benchmark suite

Cf. the one in IntervalRootFinding.jl

Add maximum number of function evaluation ?

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.

Add mean-value form and 2nd-order form including Hessian

Maybe in IntervalRootFinding.jl.

Add Goldstein-Price test case:
#11

Document the cross operator for multivariate intervals

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.

Test suite in JuMP syntax

https://lanl-ansi.github.io/MINLPLibJuMP.jl/latest/

`minimise` error because not minimum found

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

Goldstein-Price test case

minimise struggles with this.

https://github.com/JuliaOpt/SumOfSquares.jl/blob/master/examples/goldsteinprice.jl

Error on windows with julia 1.6

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.

Bug when using BigFloat interval

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})

Need advice on how to accelerate a problem of 5 variables

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.

Survey of global optimization test functions

https://arxiv.org/pdf/1308.4008.pdf

Doc link is broken

I am not quite sure what's wrong with it: the url looks correct.

Problem with Broadcasting

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]])

Which gives

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

Unify boxes in result by default?

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

add `minmax` function

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?

Problem with maximize

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

Discontinuous function optimisation

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.

SortedVector.jl is not giving correct result.

Fails to optimize expressions containing a sup/inf on an interval valued function.

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

American english aliases for methods

For us uncivilized folk in the north, minimise and maximise might prove to be confusing

juliaintervals / intervaloptimisation.jl Goto Github PK

intervaloptimisation.jl's People

Contributors

Stargazers

Watchers

Forkers

intervaloptimisation.jl's Issues

Example:

Proposal:

Recommend Projects

Recommend Topics

Recommend Org