Support Gurobi

[mlubin]

Add optional dependency and support for @lindahua's interface for Gurobi.

[lindahua]

Thanks. Let me know if need anything on my part.

[mlubin]

@carlobaldassi: Is it a good time to standardize linprog and mixintprog?

[IainNZ]

Very cool @lindahua !!

[mlubin]

@carlobaldassi, @lindahua, @IainNZ, @timholy:

What do you think of the following standard for linprog:

solution = linprog(c, A, b, sense, lb, ub[, kwargs])

where
c is the objective vector, always in the sense of minimization
A is the constraint matrix (inequalities and equalities combined. the matlab way of splitting them really annoys me)
b is the right-hand side vector. Tuple elements (l,u) are allowed, in which case the constraint is double sided, i.e. l <= a'x <= u.
sense is the vector of constraint sense characters: '<', '=', and '>'. If the corresponding element of b is a tuple, the sense is ignored.
lb is the vector of lower bounds on the variables
ub is the vector of upper bounds on the variables

A scalar is accepted for the b, sense, lb, and ub arguments, in which case its value is replicated. -Inf and Inf should be treated correctly when given for lb or ub.

A shortened version should be available as

linprog(c, A, b, sense[, kwargs]) = linprog(c, A, b, sense, 0, Inf[, kwargs])

kwargs contains solution parameters as keyword arguments. We'll have to wait and see what the syntax will look like, but this will be in Julia soon enough.

solution should be an instance of some type with the following members

type LinprogSolution
    status
    objval
    sol
    redcost
    lambda
end

where
status is a termination status string, one of "Optimal", "Primal Infeasible", "Dual Infeasible", "Stopped due to user limit" (iteration limit or timeout), "Stopped due to errors", maybe others.
If status is "Optimal", the other members have the following values
objval - optimal objective value
sol - primal solution vector
redcost - dual multipliers for active variable bounds (zero if inactive)
lambda - dual multipliers for active linear constraints (equalities are always active)
colbasis - optimal simplex basis statuses for the variables (columns) if available. Possible values are "NonbasicAtLower", "NonbasicAtUpper", "Basic", "Superbasic"
rowbasis - optimal simplex basis statuses for the constraints (rows) if available. Same statuses.

By convention the dual multipliers should have the sign following the interpretation of marginal change in the objective when the corresponding active right-hand side is increased. This corresponds to the standard that reduced costs should be nonnegative when a variable is at a lower bound and nonpositive when a variable is at an upper bound. Different solvers might have different conventions for the lambda vector, so tranformations might be needed.

It will be easy to extend this to mixintprog once we've settled down on linprog. We may also have to iterate on things like status strings and basis strings depending on what solvers provide.

@lindahua: Gurobi doesn't directly support double-sided constraints (afaik), but most other solvers do, and I think it's an important feature. To make things easy for now, the gurobi interface can throw an error in this case, but eventually it should be wrapped to work.

[timholy]

My only concern about this interface is that b and sense encode overlapping information, and that some of this information is ignored in some cases. I'm wondering whether it would be better to make b an array of type AbstractBound, something like this:

abstract AbstractBound{T}
immutable LowerBound{T} <: AbstractBound{T}
    val::T
end
immutable UpperBound{T} <: AbstractBound{T}
    val::T
end
immutable DoubleBound{T} <: AbstractBound{T}
    lower::T
    upper::T
end

[lindahua]

@mlubin Thanks for outlining an interface.

For Gurobi, as of the latest version (v5.0+), it supports double sided constraint directly. In particular, they call it range constraint. The C functions for this are GRBaddrangeconstr and GRBaddrangeconstrs. I will add it to the Julia port soon.

[lindahua]

Some minor suggestion of the interface:

1. I guess we need more thoughts about how to express constraints consistently.

for example, we have this set of linear constraints:

a <= A x <= b
Cx <= c
Aeq x == beq

It is possible to merge these into one constraint as

[a; -inf; beq] <= [A; C; Aeq] x <= [b; c; beq] 

But let the user/caller do this may not be the optimal interface.

What about a two level interface? one for the solver implementers, and the other for the caller/user.

On the implementer side, an interface can be

linprog(f, constraint[, params])
linprog(f, (constraint1, constraint2, ...)[, params])

Here, constraints can be instances of the following types

type UpperBoundConstraint   // the internal definition is omitted
type LowerBoundConstraint
type RangeConstraint  // or DoubleSided ?
type EqualityConstraint

The implementer can all relevant functions to add these constraints to their models (e.g. Gurobi have different functions to add different types of constraints).

One the caller side, an interface can be

linprog(f, (a, '<', A, '<', b), (C, '<', c), (Aeq, '=', aeq), ...)

You can also provide macros and other higher level interface to make this even simpler and more comfortable to write.

Anyway, of the first priority is to standardize the implementer side's interface, such that package maintainer can do their job to make their packages conform to this interface.

1. For the solution type, this seems to be tailored to the simplex method. In Gurobi, the user can choose to use other ways to solve the problem (e.g. Barrier method).

I think only objv, sol, and status should be a must for all solvers. What about a modification as follows?

type LinprogSolution
     objv
     sol
     status
     attrs  # a dictionary to contain other relevant attributes such as lambda and colbasis
end

In this way, users can enquire whether certain attributes are available and get them when they do. Of course, we can standardize the name for commonly used attributes.

Also, I think the type of status should be symbols (Julia use symbols for enumerators).

[lindahua]

Some of these types can be immutable. But we can defer this later.

[mlubin]

@timholy: That's a good point, but the AbstractBound approach seems like it could be confusing to new users and also a bit ugly to construct. I tried to make the interface accessible to someone with a matlab-level understanding of Julia. It would be ok though if using ranges required more julian syntax, because this is an advanced feature that matlab can't do, but I'd like the most basic syntax to be completely straightforward.

What if we add a 'r' sense character for ranges. So if sense[i] == 'r', the lower bound is b[1] and the upper bound is b[2]. If sense[i] != 'r', then we can @assert isa(b[i],Real). This way there's no ignored information and we throw an error if conflicting input is given.

@lindahua: Interesting idea. Given the interface for solver implementers, we could make a unified linprog package that provides the caller/user interface. I.e. the caller/user interface can be implemented on top of the solver-implementer interface in a package-independent way, instead of requiring two different package-dependent interfaces. Given that, let's focus on the solver-implementer interface as you suggested. I still think the solver-implementer interface should be matlab-user friendly because it's the most direct replacement. I agree with your changes to LinprogSolution.

[mlubin]

Alternatively we can make the solver-implementer interface very precise and Julian and push off the user friendly magic to a linprog package. This could make things simpler on the solver side.

[mlubin]

Let's continue this discussion on the mailing list: https://groups.google.com/forum/#!topic/julia-users/23PzvW_2nvs

[mlubin]

We can now switch to Gurobi with MathProgBase.setlpsolver(Gurobi) and MathProgBase.setmipsolver(Gurobi) (syntax subject to change in the future).