Talked to Juan Pablo today, and he expressed interest in adding incumbent and branching callbacks at some point. Are these supported in Gurobi/GLPK?

The following test of the fallback conic solver fails. (It works when using ECOS.) (I'm using the integers branch of Convex here.)

using Base.Test
using Convex
using GLPKMathProgInterface

TOL = 1e-2

x = Variable()
p = minimize(x, x>=4.3)
solve!(p, GLPKSolverLP())
@test_approx_eq_eps p.optval 4.3 TOL

(It thinks the optimal value is 0.)

The assertion here errors out if, e.g. the row and column index vectors are not identical, even if the off-diagonal terms have zero coefficients. This comes up in the vectorized stuff in JuMP, when we do something like x'Q*x. Should we relax this?

I have an issue when using quadprog with both CPLEX and Gurobi, this seems to be due to the fact that they don't support range constraints, i.e. either rowlb or rowub has to be Inf. This would be fine but unfortunately, if I use senses and my problem data are integers then the one-sided constraints are converted to the smallest (or biggest) Int64, not Inf which makes both CPLEX als well as Gurobi trip up.
This should be fixed in quadprog.jl:29 and 33, by replacing typemin(realtype) and typemax(realtype) with -Inf and Inf respectively.
The individual solvers should have to decide how they want to deal with Inf.

In the README the solver parameter to linprog is presented as a keyword argument, but in the code it is a positional argument with a default.

Is there a specific reason for the following definition in
MathProgBase/src/HighLevelinterface/quadrpog.jl

quadprog(c,Q,A,rowlb,rowub,solver::AbstractMathProgSolver=defaultQPsolver) = quadprog(c,Q,A,rowlb,rowub,0,Inf,solver)

I would expect instead something like:

quadprog(c,Q,A,rowlb,rowub,solver::AbstractMathProgSolver=defaultQPsolver) = quadprog(c,Q,A,rowlb,rowub,-Inf,Inf,solver)

described my MOSEK http://docs.mosek.com/modeling-cookbook/cqo.html
and Vandenberghe http://www.seas.ucla.edu/~vandenbe/236C/lectures/conic.pdf

encourage more solvers to implement power cone!

I can ask whether the objective function is linear, but is there any better way of easily getting the objective coefficients (ideally as a sparse vector, or equivalently the index and value vectors) than evaluating the gradient?

And JuMP.prepConstrMatrix is a useful feature of providing the constraint matrix for the linear set of constraints, but it doesn't appear to be exposed through MathProgBase for nonlinear problems. Maybe this should be a JuMP issue of providing an impementation of buildInternalModel for nonlinear that does what prepConstrMatrix does, then MathProgBase.getconstrmatrix could be used for this?

Will need AbstractNLPEvaluator in order to access objective and constraint expressions. It's of course possible to implement a writer as solve(m), but also fairly awkward when writing and solving are separate steps. Could come up with a contrived set of "solver" options for whether or not to actually run the solver, and which sub-solver and other options are selected, but this would really not be the most elegant interface.

Also probably want a version of loadnonlinearproblem! that can populate an AbstractMathProgModel and an AbstractNLPEvaluator from a filename, to load optimization problems that were created and saved in other formats.

I ran Pkg.test on the MathProgBase and ran into following error after a couple of Ipopt tests.

ERROR: LoadError: UndefVarError: model not defined
 in test_conic_fallback at /Users/arora/.julia/v0.4/MathProgBase/test/runtests.jl:20
 in include at ./boot.jl:261
 in include_from_node1 at ./loading.jl:304
 in process_options at ./client.jl:280
 in _start at ./client.jl:378
while loading /Users/arora/.julia/v0.4/MathProgBase/test/runtests.jl, in expression starting on line 35
=============================================================================[ ERROR: MathProgBase ]==============================================================================

failed process: Process(`/Users/arora/src/julia/usr/bin/julia --check-bounds=yes --code-coverage=none --color=yes /Users/arora/.julia/v0.4/MathProgBase/test/runtests.jl`, ProcessExited(1)) [1]

==================================================================================================================================================================================
INFO: Removing Cbc v0.2.0
INFO: Removing Clp v0.2.0
ERROR: MathProgBase had test errors
 in error at ./error.jl:21
 in test at pkg/entry.jl:803
 in anonymous at pkg/dir.jl:31
 in cd at file.jl:22
 in cd at pkg/dir.jl:31
 in test at pkg.jl:71

My Julia setup is:

Julia Version 0.4.3-pre+6
Commit adffe19 (2015-12-11 00:38 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin14.1.1)
  CPU: Intel(R) Core(TM) i7-3840QM CPU @ 2.80GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

and I using MathProgBase 0.4.0

thanks,
Nitin

Current API looks great. But, we still have to document this (the semantics of solver & model, as well as the methods acting upon them).

In Gurobi at a MIP node, cbgetobj returns the objective of the proposed solution (which may not be feasible), not the objective of the "current best integer-feasible solution" as the docs state. We should check this behavior across solvers and update the documentation.

The test hasn't been updated for the new symmetric definition of the SDP cone.

We should add documentation specifying that part of the API contract is that solver interfaces (and solvers, I guess) will not alter user data (except a mathprog model) when passed through an interface function, even one with a bang. Ref jump-dev/CPLEX.jl#31

getconstrmatrix has an entry in the docs, but it doesn't appear at http://mathprogbasejl.readthedocs.org/en/latest/mathprogbase.html

This is just an API idea I was having: Instead of having functions like linprog, mixintprog, setmipsolver, etc, there could also be one method solve(problem, solver) with immutable problem and solver types.

Example:
Instead of writing

setlpsolver(:GLPK)
solution = linprog(c, A, sense, b, lb, ub; options...)

one would write

solution = solve(LP(c, A, sense, b, lb, ub), GLPK(solveroptions))

MathProgBase.jl would then provide the immutable types LP, MIP, LinprogSolution, etc, and GLPK.jl would depend on MathProgBase.jl and provide the immutable type GLPK as well as the functions solve(problem::LP, solver::GLPK) and solve(problem::MIP, solver::GLPK).

I think we should change names of mutating methods such as addvar, addconstr, and setobj to addvar!, addconstr! and setobj!. There are some others that fall in this family.

I hope this is not too much work (at this stage).

I think it would be useful to have an addvar! method that does not require constraint coefficients:

addvar!(model, lb, ub, objcoef)

It is not uncommon in practice that people first add variables, and then constraints.

I am exploring using eval_grad_f to perform coordinate ascent steps in mean field variational models. At each step, I only need a few terms of the gradient, but eval_grad_f always returns the whole thing. Is there a way to save computation by only requesting a few select terms of the gradient to evaluate?

The abstraction is a getting a bit crowded with methods overloaded to mean different things depending on the sequence of operations you perform on a model object. E.g., if you load a nonlinear problem, then getconstrduals returns the multipliers for all constraints, while if you load a linear problem, it returns the multipliers for the linear constraints. A similar issue leads to jump-dev/JuMP.jl#472 where you're only allowed to call setquadobj! after loadproblem! but not after loadnonlinearproblem!. I think we should be a bit more careful about solver features (jump-dev/JuMP.jl#83).

Concrete proposal:
Define a finite list of problem classes which a solver might support, and force users to request the features that they need as arguments to model(). Example:

m = model(IpoptSolver(), Nonlinear)

or

m = model(IpoptSolver(), LP)

Having these as arguments to model lets the solver wrapper potentially return different types depending on the input. The major problem classes might be LP, Conic, and Nonlinear. These are meant to be very broad but capture the different styles of input that we currently support. Further distinctions can be made by using applicable as we do now. (Extra bonus if we can resolve the GLPKSolverMIP/GLPKSolverLP split.) These changes would most directly affect JuMP; Convex.jl shouldn't need more than a couple lines of code tweaked.

CC @carlobaldassi @tkelman @madeleineudell

I'm new to this package, so I tried to start with a simple case, the binary knapsack problem solution shown in the [manual](https://mathprogbasejl.readthedocs.org/en/latest/mixintprog.html#mixed-integer-programming manual). Unfortunately on my installation, this example throws a "no method matching" error:

julia> using MathProgBase

julia> using Cbc

julia> mixintprog(-[5.,3.,2.,7.,4.],[2. 8. 4. 2. 5.],'<',10,'I',0,1)
ERROR: `mixintprog` has no method matching mixintprog(::Array{Float64,1}, ::Array{Float64,2}, ::Char, ::Int64, ::Char, ::Int64, ::Int64)

julia> versioninfo()
Julia Version 0.3.7
Commit cb9bcae* (2015-03-23 21:36 UTC)
Platform Info:
  System: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i7-3610QM CPU @ 2.30GHz
  WORD_SIZE: 64
  BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

julia> Pkg.installed("MathProgBase")
v"0.3.14"

julia> Pkg.installed("Cbc")
v"0.1.7"

I'll study the source to gain some insight and will update this issue with a suggestion for an example that works when I succeed in running mixintprog. Meanwhile please feel free to let me know if I'm making some silly newbie error in this attempt.

I'd like to add the setwarmstart! functionality to SCS.jl. For the purposes of warmstarting SCS, it would be quite useful to have a primal & dual warmstart. But doing this requires a change to the MathProgBase interface, I think. For a minimally breaking change, I propose the following signature:

setwarmstart!(m::MathProgModel, v; kwargs...)

which, in SCS.jl, would be implemented as

function setwarmstart!(m::SCSMathProgModel, primal_sol; dual_sol=Float64[], kwargs...)
    m.primal_sol = primal_sol
    m.dual_sol = dual_sol
    m
end

while all other solvers would ignore all keyword arguments. (For now, at least.)

Are there any solutions that don't require modifying the signature of setwarmstart! in all solvers that have implemented it?

setquadobjterms!(m::AbstractMathProgModel, rowidx, colidx, quadval)
Provides an alternative “terms”-based interface to setquadobj!. A list of quadratic terms is specified instead of the matrix Q. For example, the objective x21+2x1x2 is specified by setquadobjterms!(m,[1,1],[1,2],[1.0,2.0]). Duplicate terms are summed together. Note: this method does not need to be implemented by solvers.

It doesn't need to be implemented, but its what JuMP uses...

Right now, MathProgBase supports warmstarts only for MIP solvers. However, it could be extremely useful for other solvers as well: for example, SCS supports warmstarts. Iain recommends relaxing the text in MathProgBase documentation so it doesn't specify MIP, then just continue not implementing it for solvers that don't support it, as we are now. At that point we can add it to the interface for non-MIP solvers that do support it.

Currently, each solver package is required to export a model method.

If we use two different solver packages at the same time, it is going to cause trouble -- which model method should be invoked ?

I suggest using the first argument for dispatch:

# note that the first argument should be the type of the solver
# instead of an solver instance.
m = model(GurobiSolver; options...)  # create a gurobi model
m = model(GLPKSolver; options...)  # create a GLPK model

# ordinary user can use a default version as follows
# this should translate to model(default_solvertype, options...);
m = model(;options...)

This approach has been successfully used in other packages like Distributions.jl.

There are three approaches I can think of:

1. Check the return value and exit if == :Exit
2. Add a new interface function (abortcallback?) that exits
3. Do nothing and just use errors in Julia

I think I prefer 2), since 1) can easily lead to type instability. Thoughts? @mlubin @IainNZ

Given the discussion in cvxgrp/scs#31, it seems like we should revise the definition of the SDP cone in the conic interface to contain only the lower triangle ( n*(n+1)/2 variables ). This is the format that Mosek supports anyway.
@madeleineudell @karanveerm @ulfworsoe

Can probably just start with Lorentz and SDP cones. Maybe something like

idx = addLorentz(m, dim)
addconicconstr(m, idx, varidx, xnidx)

for \sqrt{\sum_{i \in varidx} x_i^2} \leq x_xnidx

Is there a way to pass user defined immutable type arrays to the functions (as additional function arguments or vector ) ? I need those to evaluate my cost function and gradients. I am trying to use the IPOPT.jl package.

thanks,
Nitin

In test/mixintprog the line
@test_approx_eq norm(sol.sol[1:5] - [1.0, 0.0, 0.0, 1.0, 1.0]) 0.0
tests the result of an integer optimization, but the epsilon value used in the comparison is way too small (~ 5e-28 ?). Gurobi uses a default IntFeasTol=1e-5, which means that the correct epsilon value should be at least
norm([IntFeasTol for i=1:5])
in this case. For continuous cases the relevant parameter is FeasibilityTol (= 1e-6) in Gurobi.

Would be nice to get the number of explored branch and bound at the end of a solve without callbacks.

In the callback docs it's written about cbaddcut! and cbaddlazy!:

Adds cut to model, as specified by a sparse representation identical to addconstr!

but the function signature is different, the last two arguments are sense and rhs, rather then lb and ub as in addconstr!.

Also, it's only when looking in the JuMP or CPLEX code that one sees that sense is one of :==, :<= or :>=, which BTW seems inconsistent with other parts of MathProgBase where chars (e.g. '>') are used. I think that if we keep this signature, we should either be consistent and use chars everywhere (fixing JuMP and solvers accordingly), or allow both chars and symbols. But personally I'd rather actually change the signature to match that of addconstr!.

using JuMP, MathProgBase, Gurobi
m = Model(solver=GurobiSolver())
@defVar(m, x)
@addConstraint(m, norm(x) <= 1)
buildInternalModel(m)
MathProgBase.getrawsolver(getInternalModel(m))

ERROR: MethodError: `getrawsolver` has no method matching getrawsolver(::MathProgBase.SolverInterface.LPQPtoConicBridge)

You can work around it by accessing the fields directly (MathProgBase.getrawsolver(getInternalModel(m).lpqpmodel)), but I think the abstraction should probably do that for you. CC @mlubin

I have noticed this that appears as a type instability in AbstractNonlinearModel. Specifically, using the HS071() problem in the test directory

function nlptest(solver=MathProgBase.defaultNLPsolver)

    m = MathProgBase.NonlinearModel(solver)
    l = [1,1,1,1]
    u = [5,5,5,5]
    lb = [25, 40]
    ub = [Inf, 40]
    MathProgBase.loadproblem!(m, 4, 2, l, u, lb, ub, :Min, HS071())
    MathProgBase.setwarmstart!(m,[1,5,5,1])

    MathProgBase.optimize!(m)
    stat = MathProgBase.status(m)

    @test stat == :Optimal
    x = MathProgBase.getsolution(m)
    @test_approx_eq_eps x[1] 1.0000000000000000 1e-5
    @test_approx_eq_eps x[2] 4.7429996418092970 1e-5
    @test_approx_eq_eps x[3] 3.8211499817883077 1e-5
    @test_approx_eq_eps x[4] 1.3794082897556983 1e-5
    @test_approx_eq_eps MathProgBase.getobjval(m) 17.014017145179164 1e-5

    # Test that a second call to optimize! works
    MathProgBase.setwarmstart!(m,[1,5,5,1])
    MathProgBase.optimize!(m)
    stat = MathProgBase.status(m)
    @test stat == :Optimal
    m
end

mm = nlptest()

[OMITTED...]
EXIT: Optimal Solution Found.
Ipopt.IpoptMathProgModel(Ipopt.IpoptProblem(Ptr{Void} @0x00007ff7fbf551f0,4,2,[1.0,4.742999641809297,3.8211499817883072,1.3794082897556983],[24.999999754436296,40.00000000001771],[-0.5522936588816064,0.16146856313488814],[1.0878712258659022,6.693166200640359e-10,8.887657145295419e-10,6.570872591660427e-9],[6.262653086171725e-10,9.788835007031796e-9,2.1228492520638034e-9,6.925197858854474e-10],17.014017145179164,0,eval_f_cb,eval_g_cb,eval_grad_f_cb,eval_jac_g_cb,eval_h_cb,nothing,:Min),4,0,[1.0,5.0,5.0,1.0],Any[])

u() = mm.inner.x
@code_warntype u()

Variables:

Body:
  begin  # none, line 1:
      return (top(getfield))((top(getfield))(Main.mm,:inner)::Any,:x)::Any
  end::Any

##### Current package: MathProgBase
INFO: Installing BinDeps v0.2.12
INFO: Installing Cbc v0.0.6
INFO: Installing MathProgBase v0.1.6
INFO: Installing URIParser v0.0.1
INFO: Building Cbc
INFO: REQUIRE updated.
INFO: Installing Clp v0.0.6
INFO: Building Cbc
INFO: REQUIRE updated.
INFO: No packages to install, update or remove.
INFO: REQUIRE updated.
ERROR: No QP solver detected. Try installing one of the following packages: "Gurobi",  "CPLEX",  "Mosek",  and restarting Julia
 in error at error.jl:21
 in model at /home/idunning/pkgtest/.julia/v0.2/MathProgBase/src/MathProgSolverInterface.jl:58
 in quadprog at /home/idunning/pkgtest/.julia/v0.2/MathProgBase/src/quadprog.jl:10
 in quadprogtest at /home/idunning/pkgtest/.julia/v0.2/MathProgBase/test/quadprog.jl:8
 in quadprogtest at /home/idunning/pkgtest/.julia/v0.2/MathProgBase/test/quadprog.jl:6
 in include at boot.jl:238
at /home/idunning/pkgtest/.julia/v0.2/MathProgBase/test/runtests.jl:10
INFO: No packages to install, update or remove.
INFO: REQUIRE updated.
INFO: No packages to install, update or remove.
INFO: REQUIRE updated.
INFO: Removing BinDeps v0.2.12
INFO: Removing Cbc v0.0.6
INFO: Removing Clp v0.0.6
INFO: Removing MathProgBase v0.1.6
INFO: Removing URIParser v0.0.1
INFO: REQUIRE updated.

I guess this means we can't automatically test MathProgBase and should except it, right?

Loading default solvers for every problem class at the using MathProgBase invocation is getting a bit burdensome. It seems reasonable to defer this to a call to optimize!, and only load the solver needed for whatever problem class you are calling optimize! on.

Normally I'd just hack this up and submit a PR, but I'm feeling short of time so will first float a suggestion and see what people think.

For nonlinear optimization, the API is good and flexible, but to me it still feels like there's more boilerplate required from the user than would be ideal. I'm wondering whether it would be better to define a traits-based API, e.g.,

abstract NLConstraintsType
immutable NLConstraintsNone <: NLConstraintsType end
immutable NLConstraintsJac <: NLConstraintsType end
immutable NLConstraintsJacProd <: NLConstraintsType end

nlconstraints(d::AbstractNLPEvaluator) = NLConstraintsJac()   # default value

# fallback definition, so the user doesn't have to define this
SolverInterface.eval_g(d, g, x) = SolverInterface.eval_g(nlconstraints(d), d, g, x)
SolverInterface.eval_g(::NLConstraintsNone, d, g, x) = nothing

# do the same thing for jac_structure, eval_jac_g
...

so the user can define

MathProgBase.nlconstraints(::MyEvaluator) = NLConstraintsNone()

and then not have to worry about any of the g-related functions. Same strategy for the various Hessian options, etc. These would essentially replace initialize and features_available, at least as far as users are concerned, and as illustrated the more granular design reduces the need for boilerplate stub functions.

In my own code I do this via inheritance, e.g.,

abstract BoundsOnly <: SolverInterface.AbstractNLPEvaluator
SolverInterface.eval_g(::BoundsOnly, g, x) = nothing

but the traits-based mechanism is far better because it allows you to mix-and-match traits arbitrarily.

cbgetbestbound should mean the best lower bound from the objective so far (when minimizing), not the objective value of the best feasible solution, which is returned by cbgetobj.

@joehuchette

Solver-specific parameters should be prefixed (e.g. GLPKPresolve) so that the same call to linprog works when the solver is changed; each solver can pick out its own parameters and ignore the others. After this is set up, we can start thinking about standardizing a common subset of parameters.
@carlobaldassi

SCS just added power cones {(x,y,z) | x^a * y^(1-a) >= |z|, x>=0, y>=0}! Let's add these to the MathProgBase conic interface. (I think this only requires a change to the docs.)

EXIT: Optimal Solution Found.
ERROR: assertion failed: |:(sol.objval) - -0.75| <= 1.1102230246251565e-12
  :(sol.objval) = -0.75000001248
  -0.75 = -0.75
  difference = 1.2480000033399108e-8 > 1.1102230246251565e-12
 in error at error.jl:22

IPOPT can't be that precise...

The methods eval_grad_f(), eval_jac_g(), eval_hesslag() and eval_hesslag_prod() all fill in one of their input arguments. By convention, it might be better to rename them eval_grad_f!(), eval_jac_g!(), eval_hesslag!() and eval_hesslag_prod!(). Hopefully, this doesn't mean too much breakage...

It would be nice to have a unified interface to grab solve statistics (e.g. solve time, iteration count, basic B&B stuff, etc). Whatever is useful and available across solvers, really. Then this info could be accessible directly through JuMP.

The documentation is ambiguous about whether getconstrLB/setconstrLB! etc.. apply to just the linear constraints or to all constraints that have been added to the model so far (including quadratic constraints). Mosek appears to count all constraints, while Gurobi counts only linear constraints. JuMP implicitly assumes the behavior of Gurobi. Unless there are objections, I propose to clarify that these functions only apply to the linear constraints.

@IainNZ @joehuchette @wflu @ulfworsoe

In test/conicinterface.jl, there are tests for getduals, but the tests are not run by default.

I am getting errors when calling getduals and using Mosek and SCS, and a starting point for fixing these would be to actually run the tests.

E.g. addvar!(m, lb, ub, coeff, vtype::Char). This will be useful for SCIP where variables are created one-by-one and the variable type is asked for during construction.

I am using MathProgBase.jl in a publication. Is there a standard citation that I may cite? If not, then should I cite JuMP?

The method getreducedcosts() returns dual variable information, but according to @mlubin is not part of the conic interface. I propose a getconicvarduals() (for lack of a better name), that returns the dual solution values. This can in principle be computed as

s = A'y+c

but as solvers work in epsilon-precision, the solver's dual variable may only be approximately equal to A'y+c. In some cases this may make a difference when validating eps-optimality/-feasibility.