Merge pull request #7 from JuliaPolyhedra/mpb

Fix MPB support

src/lp.jl

@@ -87,18 +87,32 @@ function getsolution(sol::CDDLPSolution{Cdouble})
   solutiontmp[2:end]
 end
 
-function getreducedcosts(sol::CDDLPSolution{GMPRational})
+function getconstrduals(sol::CDDLPSolution{GMPRational})
   soldata = unsafe_load(sol.sol)
   # -1 because there is the objective
-  dualsolutiontmp = myconvert(Array, soldata.dsol, soldata.m-1)
-  dualsolution = Array{Rational{BigInt}}(dualsolutiontmp)
-  myfree(dualsolutiontmp)
-  dualsolution
+  nbindex = myconvert(Array, soldata.nbindex, soldata.d+1)
+  dsol = myconvert(Array, soldata.dsol, soldata.d+1)
+  dual = zeros(Rational{BigInt}, soldata.m-1)
+  for j in 2:soldata.d
+      if nbindex[j+1] > 0
+          dual[nbindex[j+1]] = dsol[j]
+      end
+  end
+  myfree(dsol)
+  dual
 end
-function getreducedcosts(sol::CDDLPSolution{Cdouble})
+function getconstrduals(sol::CDDLPSolution{Cdouble})
   soldata = unsafe_load(sol.sol)
   # -1 because there is the objective
-  myconvert(Array, soldata.dsol, soldata.m-1)
+  nbindex = myconvert(Array, soldata.nbindex, soldata.d+1)
+  dsol = myconvert(Array, soldata.dsol, soldata.d+1)
+  dual = zeros(Cdouble, soldata.m-1)
+  for j in 2:soldata.d
+      if nbindex[j+1] > 0
+          dual[nbindex[j+1]] = dsol[j]
+      end
+  end
+  dual
 end
 
 type Cdd_LPData{T<:MyType}

src/mathprogbase.jl

@@ -10,7 +10,6 @@ type CDDPolyhedraModel <: AbstractPolyhedraModel
   objval
   solution
   constrsolution
-  reducedcosts
   constrduals
   infeasibilityray
 end
@@ -28,7 +27,7 @@ type CDDSolver <: AbstractMathProgSolver
 end
 
 function PolyhedraModel(s::CDDSolver)
-  CDDPolyhedraModel(s.solver_type, s.exact, nothing, :Undefined, 0, [], [], [], [], [])
+  CDDPolyhedraModel(s.solver_type, s.exact, nothing, :Undefined, 0, [], [], [], [])
 end
 LinearQuadraticModel(s::CDDSolver) = PolyhedraToLPQPBridge(PolyhedraModel(s))
 
@@ -56,21 +55,24 @@ function optimize!(lpm::CDDPolyhedraModel)
   lpm.objval = getobjval(sol)
   lpm.solution = getsolution(sol)
 
-  lpm.reducedcosts = getreducedcosts(sol)
-  # FIXME if A has equalities, cddlib splits them as 2 inequalities
-  lpm.reducedcosts = lpm.reducedcosts[1:nhreps(prob)]
+  lpm.constrduals = getconstrduals(sol)
+  # if A has equalities, cddlib splits them as 2 inequalities
+  m = nhreps(prob)
+  if length(lpm.constrduals) > m
+      secondeqduals = lpm.constrduals[m+1:end]
+      lpm.constrduals = lpm.constrduals[1:m]
+      lpm.constrduals[collect(linset(prob))] -= secondeqduals
+  end
   # FIXME if A is GMPRational, check that no creation/leak
 
   T = eltype(prob)
 
   lpm.constrsolution = Vector{T}(nhreps(prob))
-  lpm.constrduals = zeros(T, fulldim(prob))
   lpm.infeasibilityray = zeros(T, nhreps(prob))
 
   eps = 1e-7
   for (i,h) in enumerate(hreps(prob))
     lpm.constrsolution[i] = dot(coord(h), lpm.solution)
-    lpm.constrduals += lpm.reducedcosts[i] * coord(h)
     if Polyhedra.mygt(lpm.constrsolution[i], h.β)
       lpm.infeasibilityray[i] = -1
     end
@@ -96,7 +98,8 @@ function getconstrsolution(lpm::CDDPolyhedraModel)
 end
 
 function getreducedcosts(lpm::CDDPolyhedraModel)
-  copy(lpm.reducedcosts)
+  prob = get(lpm.prob)
+  spzeros(eltype(prob), fulldim(prob))
 end
 
 function getconstrduals(lpm::CDDPolyhedraModel)

src/matrix.jl

@@ -273,7 +273,6 @@ function setobjective{N, T}(matrix::CDDInequalityMatrix{N, T}, c, sense)
   dd_setmatrixobjective(matrix.matrix, sense == :Max ? dd_LPmax : dd_LPmin)
   obj = [zero(T); Vector{T}(c)]
   dd_copyArow(unsafe_load(matrix.matrix).rowvec, obj)
-  myfree(obj)
 end
 
 # V-representation

src/mytype.jl

@@ -153,12 +153,8 @@ polytypefor{T <: AbstractFloat}(::Type{T}) = Cdouble
 polytypefor{T <: PolyType}(::Type{T})        = T
 
 # Used by mathprogbase.jl
-function myconvert(::Type{Array}, x::Ptr{Cdouble}, n)
-  y = Array{Cdouble, 1}(n)
-  for i = 1:n
-    y[i] = unsafe_load(x, i)
-  end
-  y
+function myconvert{T<:Union{Cdouble, Clong}}(::Type{Array}, x::Ptr{T}, n)
+  copy(unsafe_wrap(Array, x, n))
 end
 function myconvert(::Type{Array}, x::Ptr{GMPRational}, n)
   y = Array{GMPRationalMut, 1}(n)

test/mathprogbase.jl

@@ -1,9 +1,12 @@
+@test_throws ErrorException CDDSolver(;solver_type=:Simplex)
+
 const mathprogbase_test = joinpath(Pkg.dir("MathProgBase"), "test")
 
-include(joinpath(mathprogbase_test, "linprog.jl"))
-linprogtest(CDDSolver())
-linprogtest(CDDSolver(;solver_type=:DualSimplex,exact=true))
-# Failing for
+include(joinpath(mathprogbase_test, "linproginterface.jl"))
+linprogsolvertest(CDDSolver())
+linprogsolvertest(CDDSolver(; solver_type=:DualSimplex, exact=true))
+
+# CrissCross is failing for the following problem of linprog.jl
 #
 #  # test unbounded problem:
 #  # min -x-y
@@ -14,6 +17,5 @@ linprogtest(CDDSolver(;solver_type=:DualSimplex,exact=true))
 #
 # CrossCross says dd_Inconsistent while it should say dd_DualInconsistent or dd_Unbounded (Unbounded)
 
-#linprogtest(CDDSolver(;solver_type=:CrissCross))
-#linprogtest(CDDSolver(;solver_type=:CrissCross,exact=true))
-@test_throws ErrorException CDDSolver(;solver_type=:Simplex)
+linprogsolvertest(CDDSolver(; solver_type=:CrissCross))
+linprogsolvertest(CDDSolver(; solver_type=:CrissCross, exact=true))

test/runtests.jl

@@ -5,5 +5,5 @@ using Base.Test
 include("simplex.jl")
 include("permutahedron.jl")
 include("board.jl")
-#include("mathprogbase.jl")
+include("mathprogbase.jl")
 include("polyhedron.jl")