[WIP] Add support for exact newton polytope #2

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.4
 JuMP
 MultivariatePolynomials
+Polyhedra

src/certificate.jl

@@ -1,3 +1,5 @@
+using Polyhedra
+
 # Inspired from SOSTools
 import Base.extrema
 export getmonomialsforcertificate, randpsd, randsos
@@ -21,9 +23,12 @@ end
 
 cfld(x::NTuple{2,Int}, n) = (cld(x[1], 2), fld(x[2], 2))
 
-# TODO sparse with Newton polytope (Polyhedra.jl for convex hull)
-function getmonomialsforcertificate(Z::MonomialVector, sparse=:No)
-    if sparse == :No
+function getmonomialsforcertificate(x::MonomialVector, parts::Vector, libs::Vector)
+    @assert length(parts) == length(libs)
+    for i in 1:length(parts)
+    n = length(x.Z[1])
+    mindeg, maxdeg = cfld(extrema(map(sum, x.Z)), 2)
+    if lib === nothing
         # Cheap approximation of the convex hull as the approximation of:
         #
         # z such that mindeg < sum(z) < maxdeg
@@ -43,19 +48,23 @@ function getmonomialsforcertificate(Z::MonomialVector, sparse=:No)
         # | +----+
         # | ^---------- minmultideg
         # +---------
-        mindeg, maxdeg = cfld(extrema(map(sum, Z.Z)), 2)
-        n = length(Z.Z[1])
         minmultideg, maxmultideg = Vector{Int}(n), Vector{Int}(n)
         for i in 1:n
-            a, b = extrema(z->z[i], Z.Z)
-            minmultideg[i], maxmultideg[i] = cfld(extrema(z->z[i], Z.Z), 2)
+            a, b = extrema(z->z[i], x.Z)
+            minmultideg[i], maxmultideg[i] = cfld(extrema(z->z[i], x.Z), 2)
         end
-        MonomialVector(vars(Z), mindeg:maxdeg, z -> reduce(&, true, minmultideg .<= z .<= maxmultideg))
+        MonomialVector(vars(x), mindeg:maxdeg, z -> reduce(&, true, minmultideg .<= z .<= maxmultideg))
     else
-        error("Not supported yet :(")
+        Zm = Matrix{Int}(length(x), n)
+        for (i, z) in enumerate(x.Z)
+            Zm[i,:] = z
+        end
+        vrep = SimpleVRepresentation(Zm)
+        newtonpoly = polyhedron(vrep, lib)
+        MonomialVector(vars(x), mindeg:maxdeg, z -> 2*z in newtonpoly)
     end
 end
-getmonomialsforcertificate(Z::Vector, sparse=:No) = getmonomialsforcertificate(MonomialVector(Z), sparse)
+getmonomialsforcertificate(Z::Vector, lib=nothing) = getmonomialsforcertificate(MonomialVector(Z), lib)
 
 function randpsd(n; r=n, eps=0.1)
     Q = randn(n,n)

test/runtests.jl

@@ -5,6 +5,8 @@ using PolyJuMP
 using FactCheck
 
 include("solvers.jl")
+include("sparse.jl")
+error()
 
 include("certificate.jl")
 include("constraint.jl")
@@ -13,11 +15,12 @@ include("motzkin.jl")
 
 # SOSTools demos
 include("sospoly.jl")
+include("sosmatrix.jl")
+
 include("sosdemo2.jl")
 include("sosdemo3.jl")
 include("sosdemo4.jl")
 #include("sosdemo5.jl")
 include("sosdemo6.jl")
-include("sosmatrix.jl")
 
 FactCheck.exitstatus()

test/simplesparse.jl

@@ -0,0 +1,15 @@
+# Example 3.95 of
+# Blekherman, G., Parrilo, P. A., & Thomas, R. R. (Eds.).
+# Semidefinite optimization and convex algebraic geometry SIAM 2013
+facts("Example 3.25") do
+for solver in sdp_solvers
+context("With solver $(typeof(solver))") do
+    @polyvar w x y z
+    p = (w^4+1)*(x^4+1)*(y^4+1)*(z^4+1)+2w+3x+4y+5z
+    m = Model(solver=solver)
+    c = @polyconstraint m p >= 0
+    status = solve(m)
+    @fact status --> :Optimal
+    #s = getslack(c)
+    #@show length(s.x)
+end; end; end

test/sparse.jl

@@ -0,0 +1 @@
+include("simplesparse.jl")