Switch to Project.toml, require Julia 1, test agreement with CUTEst #36

timholy · 2019-10-08T09:23:06Z

This is a very handy package, but it's not compatible with newer versions of JuMP because of its current registration: https://github.com/JuliaRegistries/General/blob/master/O/OptimizationProblems/Compat.toml. Tagging a new version will fix that.

One thing I'd like to do is extract the various objectives and try evaluating them at different points. How, for example, would you do this?

julia> nlp = srosenbr(2)
A JuMP Model
Minimization problem with:
Variables: 2
Objective function type: Nonlinear
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: x

julia> objective_value(nlp, [0.0, 0.1])
ERROR: MethodError: no method matching objective_value(::Model, ::Array{Float64,1})
Closest candidates are:
  objective_value(::Model) at /home/tim/.julia/packages/JuMP/iGamg/src/objective.jl:25
Stacktrace:
 [1] top-level scope at none:0

julia> aex = objective_function(nlp);

julia> value(aex, var->error(var, " not recognized"))
ERROR: MethodError: convert(::Type{Union{}}, ::Float64) is ambiguous. Candidates:
  convert(::Type{T}, arg) where T<:VecElement in Base at baseext.jl:7
  convert(::Type{T}, x::Number) where T<:AbstractChar in Base at char.jl:179
  convert(::Type{T}, x::Number) where T<:Number in Base at number.jl:7
Possible fix, define
  convert(::Type{Union{}}, ::Number)

abelsiqueira · 2019-10-08T10:20:06Z

Thanks for the update.

The problem we have here is that we use these JuMP models through NLPModelsJuMP, which uses MathProgBase. So, to update to JuMP 0.20, we'd need to Update NLPModelsJuMP to use MathOptInterface. For now we need to use JuMP 0.18 then.

To compute the objective at various points, you can use NLPModelsJuMP:

julia> model = srosenbr(2)
julia> nlp = MathProgNLPModel(model)
julia> obj(nlp, [0.0; 0.1])
2.0

timholy · 2019-10-08T10:50:30Z

Ah, that's fine then, this can wait until that's ready.

Thanks for the very helpful tip about NLPModels (I should have figured that out on my own). I have an application where I'd like to evaluate the objective for intervals as a way of computing rigorous bounds. But I get this:

julia> obj(nlp, [0.0; 0.1])
2.0

julia> using IntervalArithmetic

julia> iv(x) = prevfloat(x)..nextfloat(x)
iv (generic function with 1 method)

julia> obj(nlp, iv.([0.0; 0.1]))
ERROR: MethodError: no method matching forward_eval(::Array{Float64,1}, ::Array{Float64,1}, ::Array{ReverseDiffSparse.NodeData,1}, ::SparseArrays.SparseMatrixCSC{Bool,Int64}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Interval{Float64},1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}; user_operators=ReverseDiffSparse.UserOperatorRegistry(Dict{Symbol,Int64}(), MathProgBase.SolverInterface.AbstractNLPEvaluator[], Dict{Symbol,Int64}(), Any[], Any[], Any[]))
Closest candidates are:
  forward_eval(::AbstractArray{T,1}, ::AbstractArray{T,1}, ::AbstractArray{ReverseDiffSparse.NodeData,1}, ::Any, ::Any, ::Any, ::AbstractArray{T,1}, ::Any, ::Any, ::Any; user_operators) where T at /home/tim/.julia/packages/ReverseDiffSparse/BQuct/src/forward.jl:20
  forward_eval(::AbstractArray{T,1}, ::AbstractArray{T,1}, ::AbstractArray{ReverseDiffSparse.NodeData,1}, ::Any, ::Any, ::Any, ::AbstractArray{T,1}, ::Any, ::Any) where T at /home/tim/.julia/packages/ReverseDiffSparse/BQuct/src/forward.jl:20 got unsupported keyword argument "user_operators"
  forward_eval(::AbstractArray{T,1}, ::AbstractArray{T,1}, ::AbstractArray{ReverseDiffSparse.NodeData,1}, ::Any, ::Any, ::Any, ::AbstractArray{T,1}, ::Any) where T at /home/tim/.julia/packages/ReverseDiffSparse/BQuct/src/forward.jl:20 got unsupported keyword argument "user_operators"
Stacktrace:
 [1] forward_eval_all(::JuMP.NLPEvaluator, ::Array{Interval{Float64},1}) at /home/tim/.julia/packages/JuMP/I7whV/src/nlp.jl:445
 [2] eval_f(::JuMP.NLPEvaluator, ::Array{Interval{Float64},1}) at /home/tim/.julia/packages/JuMP/I7whV/src/nlp.jl:477
 [3] obj(::MathProgNLPModel, ::Array{Interval{Float64},1}) at /home/tim/.julia/packages/NLPModelsJuMP/Ao59I/src/mpb_model.jl:135
 [4] top-level scope at REPL[11]:1

Any thoughts?

abelsiqueira · 2019-10-08T11:16:54Z

Unfortunately, that issue is in a very large backlog because it involves understanding MOI internals. If you change the compat version of JuMP here he can merge this.

On the interval thing. Maybe you can define NLPModels.obj(nlp, x :: Vector{Interval})?

timholy · 2019-10-08T11:29:05Z

If you change the compat version of JuMP here he can merge this.

It is already set to JuMP = "~0.20".

abelsiqueira · 2019-10-08T12:32:34Z

Sorry, I meant fix at 0.18, so it uses the MPB version of JuMP.

Also fixes deprecated syntax

timholy · 2019-10-12T12:18:26Z

OK, I reverted to JuMP 0.18.

More importantly, I discovered a bug in one of the functions (nondia, fix is not in this PR), which made me wonder how carefully this has been tested. So I decided to check against CUTEst. See the addition to test/runtests.jl for info. Results:

julia> length(problem_status)
168

julia> sum(kv[2].samexval for kv in problem_status)
144

julia> sum(kv[2].sameobjval for kv in problem_status)
137

julia> sum(kv[2].samegradval for kv in problem_status)
124

Looks like there is probably stuff to fix. Is it an explicit goal to default to the same problem size as CUTEst, or is that one area where disagreement is allowed? (Note that the remaining fields redefine the problem so it has the same number of vars as CUTEst.)

dpo · 2019-10-12T14:50:11Z

Thanks for the tests. I'm wary of a comparison with CUTEst because bugs in the SIF files are very hard to detect (and there must be some). The problems here have been implemented directly from reference papers or report. I would rather compare against AMPL problems (with AmplNLReader) but not all of them have been coded up in AMPL though.

timholy · 2019-10-12T15:11:27Z

~~How does one compare with AMPL? Is there a way to automate that?~~ NVM, see below.

Here are my fixes so far, I am getting info mostly by reference to https://github.com/mpf/Optimization-Test-Problems, and with one exception those have agreed with CUTEst:

diff --git a/src/arglina.jl b/src/arglina.jl
index d0beada..c910649 100644
--- a/src/arglina.jl
+++ b/src/arglina.jl
@@ -12,7 +12,7 @@
 export arglina
 
 "Linear function with `n` parameters and `m` observations  - full rank"
-function arglina(n::Int=100, m::Int=200)
+function arglina(n::Int=100, m::Int=2n)
 
   m < n && @warn("arglina: must have m ≥ n")
   m = max(m, n)
diff --git a/src/arglinb.jl b/src/arglinb.jl
index 2b207cc..59edf0e 100644
--- a/src/arglinb.jl
+++ b/src/arglinb.jl
@@ -12,7 +12,7 @@
 export arglinb
 
 "Linear function with `n` parameters and `m` observations - rank 1"
-function arglinb(n::Int=10, m::Int=20)
+function arglinb(n::Int=10, m::Int=2n)
 
   m < n && @warn("arglinb: must have m ≥ n")
   m = max(m, n)
diff --git a/src/arglinc.jl b/src/arglinc.jl
index defb46a..5b8d531 100644
--- a/src/arglinc.jl
+++ b/src/arglinc.jl
@@ -12,7 +12,7 @@
 export arglinc
 
 "Linear function with `n` parameters and `m` observations - rank 1, zero columns and rows"
-function arglinc(n::Int=10, m::Int=20)
+function arglinc(n::Int=10, m::Int=2n)
 
   m < n && @warn("arglinc: must have m ≥ n")
   m = max(m, n)
diff --git a/src/beale.jl b/src/beale.jl
index 1141436..5e68118 100644
--- a/src/beale.jl
+++ b/src/beale.jl
@@ -22,7 +22,7 @@ function beale(args...)
     @NLobjective(
       nlp,
       Min,
-      (1.5 + x[1] * (1.0 - x[2]))^2 + (2.25 + x[1] * (1.0 - x[2]^2))^2 + (2.625 + x[1] * (1.0 - x[2]^3))^2
+      (-1.5 + x[1] * (1.0 - x[2]))^2 + (-2.25 + x[1] * (1.0 - x[2]^2))^2 + (-2.625 + x[1] * (1.0 - x[2]^3))^2
     )
 
     return nlp
diff --git a/src/broydn7d.jl b/src/broydn7d.jl
index 2bdcee4..86b3866 100644
--- a/src/broydn7d.jl
+++ b/src/broydn7d.jl
@@ -50,14 +50,14 @@ function broydn7d(n :: Int=100, p :: Float64=7/3)
 
   nlp = Model()
 
-  @variable(nlp, x[i=1:n], start=(-1.0))
+  @variable(nlp, x[i=1:n], start=1.0)
 
   @NLobjective(
     nlp,
     Min,
-    abs(1 - 2 * x[2] + (3 - x[1] / 2) * x[1])^p +
-    sum(abs(1 - x[i-1] - 2 * x[i+1] + (3 - x[i] / 2) * x[i])^p for i=2:n-1) +
-    abs(1 - x[n-1] + (3 - x[n] / 2) * x[n])^p +
+    abs(1 - 2 * x[2] + (3 - 2x[1]) * x[1])^p +
+    sum(abs(1 - x[i-1] - 2 * x[i+1] + (3 - 2x[i]) * x[i])^p for i=2:n-1) +
+    abs(1 - x[n-1] + (3 - 2x[n]) * x[n])^p +
     sum(abs(x[i] + x[i + n2])^p for i=1:n2)
   )
 
diff --git a/src/brybnd.jl b/src/brybnd.jl
index a4f1ebe..8d5ae88 100644
--- a/src/brybnd.jl
+++ b/src/brybnd.jl
@@ -44,7 +44,7 @@ function brybnd(n :: Int=100; ml :: Int=5, mu :: Int=1)
 
   nlp = Model()
 
-  @variable(nlp, x[i=1:n], start=(-1.0))
+  @variable(nlp, x[i=1:n], start=1.0)   # Note: some start at -1, e.g., https://github.com/mpf/Optimization-Test-Problems/blob/master/cute/brybnd.mod
 
   @NLobjective(
     nlp,
@@ -54,11 +54,12 @@ function brybnd(n :: Int=100; ml :: Int=5, mu :: Int=1)
         x[i] * (2 + 5 * x[i]^2) + 1 -
         sum(
           x[j] * (1 + x[j])
-          for j = max(1, i-ml) : min(n, i+mu) if j != i
+          for j = (max(1, i-ml) : min(n, i+mu)) if j != i  # does this != syntax actually work?
         )
       )^2 for i=1:n
     )
   )
+  error("this is broken")
 
   return nlp
 end
diff --git a/src/nondia.jl b/src/nondia.jl
index 0e44632..cd2a859 100644
--- a/src/nondia.jl
+++ b/src/nondia.jl
@@ -33,7 +33,7 @@ function nondia(n :: Int=100)
     @NLobjective(
       nlp,
       Min,
-           (x[1] - 1.0)^2 + sum((100.0*x[1] - x[i-1]^2)^2 for i=2:n)
+      (x[1] - 1.0)^2 + 100*sum((x[1] - x[i-1]^2)^2 for i=2:n)
     )
 
     return nlp
diff --git a/src/scosine.jl b/src/scosine.jl
index a387734..c2ccb41 100644
--- a/src/scosine.jl
+++ b/src/scosine.jl
@@ -26,10 +26,7 @@ function scosine(n :: Int=100)
     n = max(2, n)
 
     nlp = Model()
-    p = zeros(n)
-    for i=1:n
-      p[i] = exp(6.0 * (i-1) / (n-1))
-    end
+    p = [exp(12.0 * (i-1) / (n-1)) for i = 1:n]
 
     @variable(nlp, x[i=1:n], start=1.0/p[i])

timholy · 2019-10-12T15:13:14Z

Ah, now I see there is https://github.com/JuliaSmoothOptimizers/AmplNLReader.jl. I don't have ampl, unfortunately.

dpo · 2019-10-12T17:21:03Z

Many thanks for checking!

beale is definitely buggy. For consistency with the Moré-Garbow-Hillstrom paper, the objective should really be written (1.5 - x[1] * (1.0 - x[2]))^2 + (2.25 - x[1] * (1.0 - x[2]^2))^2 + (2.625 - x[1] * (1.0 - x[2]^3))^2
broydn7d is correct (the AMPL model is erroneous) based on Toint's paper
brybnd is also correct based on the MGH paper
nondia is buggy but should be (1.0 - x[1])^2 + 100*sum((x[1] - x[i]^2)^2 for i=2:n) based on the Buckley report
scosine is also correct based on the Luksan-Vlcek report

This is all quite difficult to check because there are variants of all these problems in the literature and existing models. Checking is mostly done by visual inspection. That's why we tried to give extensive references in the comments.

timholy · 2019-10-13T09:29:05Z

Wow, really interesting! Is there any kind of central place for reporting errors in AMPL models? That seems pretty important, given that people use what are nominally the same test problems across many papers. (And, do they come with version numbers so one can report which set of bugs one was working with? 😄 )

Is it safe to rule out the possibility of a typo in a paper, and that test suites have incorporated corrections? Seems like that should be clearly documented if it happened, so I suppose it's safe to take the paper as the ultimate authority unless there is a clear comment to the contrary.

For the record I provided the list in #37. I've reduced this to the fixes. What do you want to do about the part of this that is a comparison with CUTEst? Keep it or discard it?

timholy · 2019-10-13T11:00:55Z

Found more errors (morebv, powellsg, and sinquad), fixed those too.

dpo · 2019-10-13T17:04:53Z

src/morebv.jl

@@ -41,7 +41,9 @@ function morebv(n :: Int=100)
    @NLobjective(
      nlp,
      Min,
-	    sum((2.0 * x[i] - x[i-1] - x[i+1] + (h^2 / 2.0) * (x[i] + (i - 1) * h + 1)^3)^2 for i=2:n-1)
+        sum((2.0 * x[i] - x[i-1] - x[i+1] + (h^2 / 2.0) * (x[i] + (i - 1) * h + 1)^3)^2 for i=2:n-1) +


Thanks for catching this! Based on Buckley's report, I think this term should have i * h instead of (i - 1) * h, and the last end-point term should have n * h instead of (n - 1) * h. Also h = 1 / (n+1). Lots of errors here!

Thanks for the further improvements. I now think the initialization setting x0[1] = x0[n] = 0 was a misreading of the problem, and that http://www.cs.cas.cz/matonoha/download/V1081.pdf suggests starting all points at 0.5. I made that further change. But note that https://github.com/mpf/Optimization-Test-Problems/blob/master/cute/morebv.mod seems to use a very different initialization.

Ok. Let's add a note about the initial point that you choose, if you don't mind.

dpo · 2019-10-13T17:16:47Z

test/runtests.jl

+    x0 = getvalue.(omodel[:x])
+    d = JuMP.NLPEvaluator(omodel)
+    MathProgBase.initialize(d, [:Grad])
+    x0 = getvalue.(omodel[:x])


We created NLPModelsJuMP for this.

Interestingly, I get an OutOfMemoryError() with NLPModelsJuMP but not with JuMP.NLPEvaluator.

dpo · 2019-10-13T17:18:58Z

test/runtests.jl

+    samenvars, samexval, sameobjval, samegradval = true, true, true, true
+    isok = false
+    if isa(x0, AbstractVector)
+        if length(x0) != length(cmodel.meta.x0)


Here you can query cmodel.meta.nvar.

dpo · 2019-10-13T17:19:29Z

test/runtests.jl

+            d, x0 = optpr_problem(probname, length(cmodel.meta.x0))
+        end
+        if (x0 ≈ cmodel.meta.x0)
+            samexval = false


Not sure I understand here. Is this back to front?

dpo · 2019-10-13T17:20:52Z

test/runtests.jl

+        gc, go = sort!(grad(cmodel, x0)), sort!(grad(d, x0))
+        if !(gc ≈ go)
+            samegradval = false
+            c = gc \ go


What's happening here?

Trying to see if it's just linearly rescaled, e.g., if f(x) -> c*f(x) then the gradient components would be multiplied by a constant. I'll get rid of that.

dpo · 2019-10-13T17:28:38Z

Thanks for combing through the problems so carefully! Here are some answers:

Is there any kind of central place for reporting errors in AMPL models? That seems pretty important, given that people use what are nominally the same test problems across many papers. (And, do they come with version numbers so one can report which set of bugs one was working with? 😄 )

The only place I can think of is mpf's repository. Unfortunately, there wasn't any good bug-tracking facility in place when those problems were created, so no version number I'm afraid. In CUTEst, we decided to leave the buggy problems untouched for historical heritage reasons (so as not to invalidate numerical results produced with them in the literature), and to create new problems for the corrected version. Not great but hey. You can find them here: https://bitbucket.org/optrove/sif/src/master/

Is it safe to rule out the possibility of a typo in a paper, and that test suites have incorporated corrections? Seems like that should be clearly documented if it happened, so I suppose it's safe to take the paper as the ultimate authority unless there is a clear comment to the contrary.

Yes, we take the paper as the reference problem formulation. If a typo were to be detected, we would probably create a new problem with the corrected version.

For the record I provided the list in #37. I've reduced this to the fixes. What do you want to do about the part of this that is a comparison with CUTEst? Keep it or discard it?

I like it, but I would suggest to use NLPModelsJuMP in there so that the comparisons still work when (one day) we upgrade to MOI. Also I would suggest to move it to a subfolder and not make it a required test because some CUTEst problems have been corrected, but the corrected problem has a different name.

…for others

timholy · 2019-10-13T21:36:19Z

OK, I have made the suggested changes. I dropped the fancy "try to detect a permutation of the coordinates" because that seems to be too difficult to do generically. This just reports on any difference.

timholy · 2019-10-13T21:45:03Z

Also worth noting that while this may correct all bugs here in the unconstrained problems, I did not check any problems with constraints. My current application is for unconstrained methods so this is all I need for now. But #37 at least lists them. Some of the disagreements are just in precision so #37 applies a bit of subjective interpretation to report the ones that I think really count.

dpo · 2019-10-13T22:08:57Z

Many thanks @timholy !

Switch to Project.toml, require Julia 1

52d92ac

Also fixes deprecated syntax

timholy force-pushed the teh/JuMP20 branch from d7956af to 7dc32fc Compare October 12, 2019 12:13

timholy changed the title ~~Switch to Project.toml, require Julia 1, JuMP 0.20~~ Switch to Project.toml, require Julia 1, test agreement with CUTEst Oct 12, 2019

timholy mentioned this pull request Oct 13, 2019

Discrepancies with other test suites #37

Open

22 tasks

timholy force-pushed the teh/JuMP20 branch from 7dc32fc to bdc67fb Compare October 13, 2019 09:28

dpo reviewed Oct 13, 2019

View reviewed changes

Modernize travis script

9fe00d3

timholy force-pushed the teh/JuMP20 branch from 8992b26 to 5dcca51 Compare October 13, 2019 21:33

timholy added 3 commits October 13, 2019 16:34

Check whether these implementations agree with CUTEst

a4d5d0f

Generalize the m argument of arglin[abc]

e489470

Fix mistakes in a pair of problems and document CUTEst discrepancies …

f40cc6f

…for others

timholy force-pushed the teh/JuMP20 branch from 5dcca51 to c634a75 Compare October 13, 2019 21:34

timholy added 4 commits October 13, 2019 16:38

Add endpoints and other fixes to morebv

ff8b099

Fix iteration in powellsg

0e36c5d

Add notes on sbrybnd

9802662

Fix a constant in sinquad

6ed51a9

timholy force-pushed the teh/JuMP20 branch from c634a75 to 6ed51a9 Compare October 13, 2019 21:38

dpo mentioned this pull request Oct 13, 2019

OufOfMemory error on power #38

Open

dpo merged commit 3079e71 into JuliaSmoothOptimizers:master Oct 14, 2019

timholy deleted the teh/JuMP20 branch October 31, 2019 14:10

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Switch to Project.toml, require Julia 1, test agreement with CUTEst #36

Switch to Project.toml, require Julia 1, test agreement with CUTEst #36

timholy commented Oct 8, 2019

abelsiqueira commented Oct 8, 2019

timholy commented Oct 8, 2019

abelsiqueira commented Oct 8, 2019

timholy commented Oct 8, 2019

abelsiqueira commented Oct 8, 2019 via email •

edited

Loading

timholy commented Oct 12, 2019

dpo commented Oct 12, 2019

timholy commented Oct 12, 2019 •

edited

Loading

timholy commented Oct 12, 2019

dpo commented Oct 12, 2019

timholy commented Oct 13, 2019

timholy commented Oct 13, 2019

dpo Oct 13, 2019

timholy Oct 13, 2019

dpo Oct 13, 2019

dpo Oct 13, 2019

timholy Oct 13, 2019

dpo Oct 13, 2019

dpo Oct 13, 2019

dpo Oct 13, 2019

timholy Oct 13, 2019

dpo commented Oct 13, 2019

timholy commented Oct 13, 2019

timholy commented Oct 13, 2019 •

edited

Loading

dpo commented Oct 13, 2019

Switch to Project.toml, require Julia 1, test agreement with CUTEst #36

Switch to Project.toml, require Julia 1, test agreement with CUTEst #36

Conversation

timholy commented Oct 8, 2019

abelsiqueira commented Oct 8, 2019

timholy commented Oct 8, 2019

abelsiqueira commented Oct 8, 2019

timholy commented Oct 8, 2019

abelsiqueira commented Oct 8, 2019 via email • edited Loading

timholy commented Oct 12, 2019

dpo commented Oct 12, 2019

timholy commented Oct 12, 2019 • edited Loading

timholy commented Oct 12, 2019

dpo commented Oct 12, 2019

timholy commented Oct 13, 2019

timholy commented Oct 13, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpo commented Oct 13, 2019

timholy commented Oct 13, 2019

timholy commented Oct 13, 2019 • edited Loading

dpo commented Oct 13, 2019

abelsiqueira commented Oct 8, 2019 via email •

edited

Loading

timholy commented Oct 12, 2019 •

edited

Loading

timholy commented Oct 13, 2019 •

edited

Loading