[RFC/WIP] MOI wrapper

[joehuchette]

This is a very rough and incomplete draft of an MOI wrapper. I'm requesting feedback particularly on the postsolve interface (the incomplete part), but welcome feedback on it all. Once we converge on how we want this to look, I'll clean it up and add more tests.

The main entrypoint is presolve!(dest::MOI.ModelLike, src::ModelLike, T::Type). You pass in a model src to be reduced, empty model dest for the reduced problem to live, and a coefficient type T. The modification happens in-place to dest, and the return value is a function which (somehow, TBD) returns the mapping between the original problem and the reduced one.

Rather than return a struct containing the termination results of the algorithm (e.g. reduced problem is infeasible, or optimal with cost X), instead I encode this in dest through trivial problems. For "optimal" reductions the problem has no variables/constraints, and a constant term in the objective equal to the optimal cost. I kind of like this since it means 1) you don't need to understand MOP-specific data structures to understand what happened, and 2) you can "solve" dest and immediately infer what's going on anyway. But, I recognize this might be a little cute.

[codecov-io]

Codecov Report

Merging #16 (064e9e1) into master (e29763f) will increase coverage by 17.70%.
The diff coverage is 83.57%.

@@             Coverage Diff             @@
##           master      #16       +/-   ##
===========================================
+ Coverage   52.19%   69.90%   +17.70%     
===========================================
  Files          15       15               
  Lines        1002     1123      +121     
===========================================
+ Hits          523      785      +262     
+ Misses        479      338      -141     

Impacted Files	Coverage Δ
src/lp/empty_column.jl	100.00% <ø> (+7.69%)	⬆️
src/lp/free_column_singleton.jl	2.89% <ø> (-1.39%)	⬇️
src/lp/row_singleton.jl	81.63% <ø> (+19.63%)	⬆️
src/presolve.jl	79.41% <33.33%> (+33.81%)	⬆️
src/moi.jl	84.67% <84.67%> (ø)
src/lp/empty_row.jl	91.66% <0.00%> (+1.66%)	⬆️
src/util.jl	100.00% <0.00%> (+8.33%)	⬆️
src/lp/dominated_column.jl	14.13% <0.00%> (+10.86%)	⬆️
... and 4 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e29763f...064e9e1. Read the comment docs.

[mtanneau]

Great! I’ll have a look during the week

[mtanneau]

OK, so I like the idea of passing a dest argument to presolve!, as it avoids additional copies of the problem.

I'm more skeptical about the Infeasible/Unbounded case.
In the NOT_INFERRED case, then a reasonable workflow could look like

model =  # Your MOI model
presolved_model =  # Your favorite Model/Optimizer

postsolve = MOP.presolve!(presolved_model, model, T)  # run presolve

MOI.optimize!(presolved_model)

# assume no error so far
x_ = MOI.get.(presolved_model, MOI.VariablePrimal(), x)  # Exact definition of x is TBD

x_original = MOP.postsolve(x_)  # or something like this.

The only difference between NOT_INFERRED and OPTIMAL/PRIMAL_INFEASIBLE/DUAL_INFEASIBLE is that, in the latter 3 cases, we already know a solution and don't need to call optimize!; postsolve operations should handle all four cases (they do need to know whether they're manipulating a ray or a point).
What I don't like with with handing out dummy infeasible/unbounded problems is the absence of correspondence to the original problem: if you optimize!(dest) and get an infeasibility certificate, how do you post-solve it?

I looked at what the main solvers are doing:

• Gurobi
  I could not find any useful documentation. There is a GRBpresolvemodel in the C API, but it's not documented anywhere. For C++ and Python, the corresponding docs are just "perform presolve on a model".
  In the C API, the function is called as GRBpresolvemodel(model[], presolved) where model and presolved are Ref{Ptr{Cvoid}}.
  From what I've tried:
  • If presolve detects unboundness or infeasibility, then presolved[] is C_NULL.
  • If all rows and columns are removed (OPTIMAL), then a problem with 0 variables and 0 constraints is returned.
• CPLEX: see CPXgetredlp docs

  It returns NULL if the problem is not presolved or if all the columns and rows are removed by presolve

I would rather see a workflow like this (with one extra step for the user):

1. User creates MOI objects model and dest
2. User calls MOP.presolve!(dest, model, T). In addition to a postsolve function, this would return some information about the problem status. For instance, the inferred status.
3. User checks presolve return status. Two cases:
   1. Status is NOT_INFERRED or OPTIMAL --> user calls optimize!(dest), gets a solution and calls postsolve.
      In the OPTIMAL case, we could give the possibility of querying a solution right-away.
   2. Status is PRIMAL_INFEASIBLE or DUAL_INFEASIBLE --> dest is useless but user can query a certificate if available
      One possibility (but I don't know whether it would be supported by solvers) would be for dest to be an empty problem with objective constant at ±Inf. We'd have to handle the postsolve accordingly though.
   3. There was an error somewhere

[joehuchette]

Interesting. Does that mean that there's no way to do postsolve through the Gurobi interface?

I think on the primal side you could have a nice interface like

postsolve, args... = MOP.presolve!(presolved_model, model, T)
MOI.optimize!(presolved_model)
x_ = MOI.get.(presolved_model, MOI.get(presolved_model, MOI.ListOfVariableIndices()))
x_original = postsolve(x_)

(At least I think this would work). I was less clear on what to do for the duals.

Do you want to return a MOP.ModelStatus that the user can check against to determine what happened? In that case, maybe the signature would be

status::MOP.ModelStatus, postsolve_primal::Function, postsolve_dual::Function = MOP.presolve!(presolved_model, model, T)

(Ideally I would also like to get rid of the T argument for presolve!, but I'll have to think about how best to do that.)

[mtanneau]

Does that mean that there's no way to do postsolve through the Gurobi interface?

Not that I know of 🤷‍♂️ Given the absence of documentation, I don't think Gurobi intended many people to use it...
I'd guess the main intended use case is debugging or factoring out presolve when doing comparisons.

Do you want to return a MOP.ModelStatus that the user can check against to determine what happened?

We can even say that NOT_INFERRED is converted to OPTIMIZE_NOT_CALLED, and return an MOI.TerminationStatusCode.
If the returned status is

• OPTIMIZE_NOT_CALLED --> proceed as expected
• OPTIMAL, PRIMAL_INFEASIBLE or DUAL_INFEASIBLE --> don't optimize!, call postsolve! directly
• Something else --> should not happen for now

Since in the second case, the resulting problem is essentially "empty", we could take the convention that postsolve! will expect a vector of size zero, and return a solution/certificate in the original space.

Then the interface could look like

status, postsolve, args... = MOP.presolve!(presolved_model, model, T)
if status == OPTIMIZE_NOT_CALLED
    # Proceed normally
    MOI.optimize!(presolved_model)
    x_ = MOI.get.(
        presolved_model,
        MOI.VariablePrimal(),
        MOI.get(presolved_model, MOI.ListOfVariableIndices())
    )
elseif status == OPTIMAL || INFEASIBLE || UNBOUNDED
    # Don't need to call optimize!, postsolve expects an empty vector
    x_ = T[]
else
    # Something went wrong
    error("Presolved exited with status $status")
end

x_original = postsolve(x_)

I'm purposefully ignoring the dual part for now.

Oh, and stating the obvious: the ordering of x_original should match the ordering of

MOI.get(model, MOI.ListOfVariableIndices())  # <--- this is the original model

after appropriate "scalarization" of any non-scalar variable indices, if any.

[joehuchette]

That sounds good to me. Using MOI.OPTIMIZE_NOT_CALLED is a bit puny, but I'll get over it :)

[mtanneau]

I won't fight for using MOI status codes over MOP ones :)

I thought not having to deal with both would be more user-friendly.

[joehuchette]

I agree that it's easier, and think it outweighs the puniness.

(I sketched out the postsolve functions. I haven't tested the code though, so they might not even compile...)

[mtanneau]

We should be able (at least internally) to use the same postsolve function for all cases.

All the postsolve does is take an initial point or ray in the presolved space, and lift it to a point or ray in the original space.
Those initial points just happen to be optimal solutions or infeasibility certificates, and they just happen to be generated by the user (after they optimize!) or the presolve directly.
The internal postsolve! functions should work whatever the input is, provided we know, primal- and dual-wise, whether we're given a point or a ray: that's the reason for having these attributes in Solution

MathOptPresolve.jl/src/solution.jl

Lines 5 to 8 in c317b25

	primal_status::SolutionStatus
	dual_status::SolutionStatus
	is_primal_ray::Bool
	is_dual_ray::Bool

As an illustration, the snippet below is an extract from these lines and these lines in Tulip:

st = presolve!(model.presolve_data)
if st == Trm_Optimal || st == Trm_PrimalInfeasible || st == Trm_DualInfeasible || st == Trm_PrimalDualInfeasible
    # Perform post-solve
    sol0 = Solution{T}(model.pbdata.ncon, model.pbdata.nvar)
    postsolve!(sol0, model.presolve_data.solution, model.presolve_data)
    model.solution = sol0
else
    # A call to optimize! will have produce a solution sol_inner
    sol_inner = ...  # <--- extract solution from presolved model
    sol_outer = Solution{T}(model.pbdata.ncon, model.pbdata.nvar)
    postsolve!(sol_outer, sol_inner, model.presolve_data)
    model.solution = sol_outer
end

In the first case, a solution is generated by presolve!, which is stored in presolve_data.solution, and is used as a seed solution in the postsolve. In the second case, the seed solution sol_inner is the one found after solving the presolved model.

[joehuchette]

@mtanneau see if you like this better

[joehuchette]

@mtanneau: Added unit tests for each possible return case.

I'm a bit confused by what's going on with unbounded models: I would think that you should be able to query an unbounded ray by passing in an empty Solution to postsolve!. Just want to make sure that my expectation is correct before trying to debug further.

[mtanneau]

I would think that you should be able to query an unbounded ray by passing in an empty Solution to postsolve!

I checked my code in Tulip: I actually intercept the unbounded/infeasible case just after presolve (here), extract the solution from pbdata, perform post-crush, and terminate the optimization.

I guess that, when detecting unboundedness/infeasibility, we just don't reduce the problem further: there is no need to eliminate x = 0 if we have 1 <= y <= 0. I just happens that we have a satisfactory solution (i.e. a certificate of infeasibility).

The tricky part would be that, in the infeasible (primal or dual) case, we have a non-trivial solution for the reduced problem (e.g. here), which is stored in pbdata.solution.

Coming back to an earlier comment, maybe a more generic approach would be to have postsolve! always expect a Solution of appropriate dimension (i.e. the dimension of the reduced problem).
In the OPTIMAL case, the reduced problem will always be empty, and the solution (in reduced space) will have dimension zero.
In the PRIMAL_INFEASIBLE/DUAL_INFEASIBLE cases, the user could either (i) use the certificate that we computed, or (ii) solve the model themselves and get a solution.

That way we eliminate the need to edge-cases in postsolve!.

Would that make more sense?

[joehuchette]

I'm a bit confused. Let's say that presolve! proves unboundedness leaves dest as a nonempty model. How should the user get the unbounded ray in reduced space? Are they expected to solve the reduced model, or should we return from presolve!? The first seems wasteful, while the second isn't generic (different return type than a model that isn't inferred).

[mtanneau]

How should the user get the unbounded ray in reduced space?

Sorry, I was not clear about that.
We would have to find a nice way to access the certificate for an infeasible model. Or one can always solve the reduced model, yes.

---

Here's my train of thought.

The last approach I mentioned essentially allows to make the following workflow work all the time

status, postsolve = MOP.presolve!(presolved_model, model, T)
MOI.optimize!(presolved_model)
x_ = MOI.get.(
    presolved_model,
    MOI.VariablePrimal(),
    MOI.get(presolved_model, MOI.ListOfVariableIndices())
)
x_original = postsolve(x_)

Of course, if the model was unbounded/infeasible at presolve, then the call to optimize! is wasteful. What I originally suggested was that, in this case, postsolve! would expect an empty vector, and return a certificate in the original space (internally obtained by post-crushing the certificate found at presolve).

For the infeasible/unbounded case, this has the drawback of (i) having a different convention, and (ii) not being able to post-crush anything other than the certificate found at presolve. Since the post-crush is just a lifting from reduced to original space, it should work for any input vertex/ray of appropriate dimension.

This is why I floated the idea of having a generic postsolve!, which however requires a way to query the infeasibility certificate for an infeasible model. For instance:

• Return a third object, which is a certificate if the model is infeasible, and garbage (e.g. nothing or an empty vector) otherwise

  status, postsolve, sol = MOP.presolve!(presolved_model, model, T)

• Since we use closures in postsolve, we could embed a certificate there

  status, postsolve = MOP.presolve!(presolved_model, model, T)
  x_ = postsolve.certificate  # garbage unless model is infeasible/unbounded

Either way, we run into the limitation of not exposing the underlying PresolveData struct, which is where all the info is being held.

[joehuchette]

What about a struct to encapsulate all this? E.g.

struct PresolveResult{T} ... end

presolve!(dest, src, T) --> PresolveResult

get_status(::PresolveResult) --> MOI.TerminationStatusCode
get_optimal_solution(::PresolveResult) --> Vector{T} # fails if get_status != OPTIMAL
get_unbounded_ray(::PresolveResult) --> Vector{T} # fails if get_status != DUAL_INFEASIBLE
get_infeasible_certificate(::PresolveResult) --> Vector{T} # fails if get_status != INFEASIBLE

post_crush(::PresolveResult, ::Vector{T}) --> Vector{T}

This would still populate dest with the presolved model, so you can still solve it and feed the result into post_crush if you'd like.

[mtanneau]

What about a struct to encapsulate all this

I like that a lot 😄

[joehuchette]

Cool, take a look now then.

[joehuchette]

Bump

[mtanneau]

Aside for the post-crushing of unbounded rays, everything looks good

[joehuchette]

Comments addressed, will merge in a day or two ( I want to turn this on in Cerberus :) )

[mtanneau]

👍 Yes, it's about time this gets merged :)