From 2ee190ee392b46172e24e5e030e965b60ace3ffb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Fri, 25 Nov 2016 14:32:16 +0000
Subject: [PATCH 1/2] Add support for exact newton polytope #2

---
 REQUIRE              |  1 +
 src/certificate.jl   | 27 +++++++++++++++++----------
 test/runtests.jl     |  5 ++++-
 test/simplesparse.jl | 15 +++++++++++++++
 test/sparse.jl       |  1 +
 5 files changed, 38 insertions(+), 11 deletions(-)
 create mode 100644 test/simplesparse.jl
 create mode 100644 test/sparse.jl

diff --git a/REQUIRE b/REQUIRE
index d081a8e9f..9eae640a7 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,3 +1,4 @@
 julia 0.4
 JuMP
 MultivariatePolynomials
+Polyhedra
diff --git a/src/certificate.jl b/src/certificate.jl
index a8b8a1682..8a2b12fbc 100644
--- a/src/certificate.jl
+++ b/src/certificate.jl
@@ -1,3 +1,5 @@
+using Polyhedra
+
 # Inspired from SOSTools
 import Base.extrema
 export getmonomialsforcertificate, randpsd, randsos
@@ -21,9 +23,10 @@ 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, lib=nothing)
+    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 +46,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)
diff --git a/test/runtests.jl b/test/runtests.jl
index ab89f5d47..1242f2faf 100644
--- a/test/runtests.jl
+++ b/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()
diff --git a/test/simplesparse.jl b/test/simplesparse.jl
new file mode 100644
index 000000000..a98f1a4da
--- /dev/null
+++ b/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
diff --git a/test/sparse.jl b/test/sparse.jl
new file mode 100644
index 000000000..33cd7a82c
--- /dev/null
+++ b/test/sparse.jl
@@ -0,0 +1 @@
+include("simplesparse.jl")

From e130fd58ca50d338d85164324dbdfd1f924fca0f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sun, 8 Jan 2017 18:47:51 +0100
Subject: [PATCH 2/2] wip

---
 src/certificate.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/certificate.jl b/src/certificate.jl
index 8a2b12fbc..9fa8d43c2 100644
--- a/src/certificate.jl
+++ b/src/certificate.jl
@@ -23,7 +23,9 @@ end
 
 cfld(x::NTuple{2,Int}, n) = (cld(x[1], 2), fld(x[2], 2))
 
-function getmonomialsforcertificate(x::MonomialVector, lib=nothing)
+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