Updates

REQUIRE

@@ -5,3 +5,4 @@ MultivariateMoments 0.1.0
 MathOptInterface 0.6
 JuMP 0.18.2+
 PolyJuMP 0.2 0.3
+Compat 1

src/SumOfSquares.jl

@@ -4,17 +4,21 @@ module SumOfSquares
 
 export SOSModel
 
+using Compat
+
 using MultivariatePolynomials
 const MP = MultivariatePolynomials
-using MultivariateMoments
 using SemialgebraicSets
+export @set
+using MultivariateMoments
 
 include("matpoly.jl")
 include("sosdec.jl")
 
 include("certificate.jl")
 
-using PolyJuMP, JuMP
+using JuMP
+using PolyJuMP
 
 include("variable.jl")
 include("constraint.jl")

src/certificate.jl

@@ -1,23 +1,6 @@
 # Inspired from SOSTools
 export getmonomialsforcertificate, randpsd, randsos
 
-function map_extrema(f::Function, itr::AbstractVector)
-    state = start(itr)
-    it1, state = next(itr, state)
-    mini = maxi = f(it1)
-    while !done(itr, state)
-        it, state = next(itr, state)
-        fit = f(it)
-        if fit < mini
-            mini = fit
-        end
-        if fit > maxi
-            maxi = fit
-        end
-    end
-    mini, maxi
-end
-
 cfld(x::NTuple{2,Int}, n) = (cld(x[1], n), fld(x[2], n))
 
 # TODO sparse with Newton polytope (Polyhedra.jl for convex hull)
@@ -44,11 +27,14 @@ function _getmonomialsforcertificate(X::AbstractVector{<:AbstractMonomial}, spar
         # +---------
         mindeg, maxdeg = cfld(extdegree(X), 2)
         n = nvariables(X)
-        minmultideg, maxmultideg = Vector{Int}(n), Vector{Int}(n)
+        minmultideg, maxmultideg = Vector{Int}(undef, n), Vector{Int}(undef, n)
         for i in 1:n
-            minmultideg[i], maxmultideg[i] = cfld(map_extrema(m -> exponents(m)[i], X), 2)
+            exponent_i(m) = exponents(m)[i]
+            minmultideg[i] = cld(Compat.mapreduce(exponent_i, min, X), 2)
+            maxmultideg[i] = fld(Compat.mapreduce(exponent_i, max, X), 2)
         end
-        monomials(variables(X), mindeg:maxdeg, m -> reduce(&, true, minmultideg .<= exponents(m) .<= maxmultideg))
+        monomials(variables(X), mindeg:maxdeg,
+                  m -> all(minmultideg .<= exponents(m) .<= maxmultideg))
     else
         error("Not supported yet :(")
     end
@@ -58,8 +44,8 @@ getmonomialsforcertificate(X::AbstractVector, sparse=:No) = _getmonomialsforcert
 function randpsd(n; r=n, eps=0.1)
     Q = randn(n,n)
     d = zeros(Float64, n)
-    d[1:r] = eps + abs.(randn(r))
-    Q' * Diagonal(d) * Q
+    d[1:r] = eps .+ abs.(randn(r))
+    Q' * Compat.LinearAlgebra.Diagonal(d) * Q
 end
 
 function _randsos(X::AbstractVector{<:AbstractMonomial}; r=-1, monotype=:Classic, eps=0.1)

src/constraint.jl

@@ -42,8 +42,8 @@ JuMP.result_dual(c::SOSConstraint) = JuMP.result_dual(c.zero_constraint)
 PolyJuMP.getslack(c::SOSConstraint) = JuMP.result_value(c.slack)
 
 function PolyJuMP.addpolyconstraint!(m::JuMP.Model, P::Matrix{PT}, ::SOSMatrixCone, domain::AbstractBasicSemialgebraicSet, basis) where PT <: APL
-    n = Base.LinAlg.checksquare(P)
-    if !issymmetric(P)
+    n = Compat.LinearAlgebra.checksquare(P)
+    if !Compat.LinearAlgebra.issymmetric(P)
         throw(ArgumentError("The polynomial matrix constrained to be SOS must be symmetric"))
     end
     y = [similarvariable(PT, gensym()) for i in 1:n]

src/sosdec.jl

@@ -36,7 +36,11 @@ function SOSDecomposition(p::MatPolynomial)
     # TODO LDL^T factorization for SDP is missing in Julia
     # it would be nice to have though
     A = getmat(p)
-    Q = chol(A)
+    @static if VERSION >= v"0.7-"
+        Q = cholesky(Matrix(A)).U
+    else
+        Q = chol(A)
+    end
     m = size(Q, 1)
     ps = [polynomial(Q[i,:], p.x) for i in 1:m]
     SOSDecomposition(ps)

src/variable.jl

@@ -40,7 +40,7 @@ ArXiv e-prints, 2017
 function constraint_matpoly! end
 
 function constraint_matpoly!(m, p::MatPolynomial, ::SOSPoly)
-    JuMP.addconstraint(m, JuMP.VectorOfVariablesConstraint(p.Q.Q, MOI.PositiveSemidefiniteConeTriangle(length(p.x))))
+    JuMP.add_constraint(m, JuMP.VectorOfVariablesConstraint(p.Q.Q, MOI.PositiveSemidefiniteConeTriangle(length(p.x))))
 end
 function constraint_matpoly!(model, p::MatPolynomial, ::SDSOSPoly)
     # `p.Q` is SDD iff it is the sum of psd matrices Mij that are zero except for
@@ -65,7 +65,7 @@ function constraint_matpoly!(model, p::MatPolynomial, ::SDSOSPoly)
 end
 function constraint_matpoly!(model, p::MatPolynomial, ::DSOSPoly)
     n = length(p.x)
-    Q = Matrix{JuMP.VariableRef}(n, n)
+    Q = Matrix{JuMP.VariableRef}(undef, n, n)
     for i in 1:n
         for j in i:n
             if i == j

test/constraint.jl

@@ -1,6 +1,9 @@
+import SemialgebraicSets
+import PolyJuMP
+
 @testset "Non-symmetric matrix SOS constraint" begin
     @polyvar x
     m = SOSModel()
     # Throws an Argument Error because it should be symmetric to be an SOS matrix
-    @test_throws ArgumentError PolyJuMP.addpolyconstraint!(m, [1 x; -x 0], SOSMatrixCone(), BasicSemialgebraicSet{Int, polynomialtype(x, Int)}(), PolyJuMP.MonomialBasis)
+    @test_throws ArgumentError PolyJuMP.addpolyconstraint!(m, [1 x; -x 0], SOSMatrixCone(), SemialgebraicSets.BasicSemialgebraicSet{Int, polynomialtype(x, Int)}(), PolyJuMP.MonomialBasis)
 end

test/extract_domain_polyopt.jl

@@ -2,7 +2,6 @@
 
 using JuMP
 using SumOfSquares
-using SemialgebraicSets
 using MultivariateMoments
 
 @testset "Polynomial Optimization example with $(factory.constructor)" for factory in sdp_factories
@@ -26,9 +25,9 @@ using MultivariateMoments
         ν = matmeasure(μ, X)
         ranktol = 1e-3
         atoms = extractatoms(ν, ranktol)
-        @test isnull(atoms) == !found
-        if !isnull(atoms)
-            η = get(atoms)
+        @test (atoms === nothing) == !found
+        if atoms !== nothing
+            η = atoms
             @test η.atoms[1].weight ≈ 1/2 atol=1e-2
             @test η.atoms[2].weight ≈ 1/2 atol=1e-2
             @test isapprox(η.atoms[1].center, [0, 1], atol=1e-2) || isapprox(η.atoms[1].center, [1, 0], atol=1e-2)

test/options_pricing.jl

@@ -10,7 +10,7 @@ using MultivariateMoments
     @polyvar x y z
     σ = [184.04, 164.88, 164.88, 184.04, 164.88, 184.04]
     X = [x^2, x*y, x*z, y^2, y*z, z^2, x, y, z, 1]
-    μ = measure([σ + 44.21^2; 44.21 * ones(3); 1],
+    μ = measure([σ .+ 44.21^2; 44.21 * ones(3); 1],
                 X)
     function optionspricing(K, cone)
         m = SOSModel(factory)

test/runtests.jl

@@ -1,12 +1,11 @@
-using Compat, Compat.LinearAlgebra, Compat.Test
+using Compat, Compat.LinearAlgebra, Compat.SparseArrays, Compat.Test
 
 using MultivariatePolynomials
-using SemialgebraicSets
 
 using MathOptInterface
 const MOI = MathOptInterface
+
 using JuMP
-using PolyJuMP
 using SumOfSquares
 
 # Taken from JuMP/test/solvers.jl

test/solvers.jl

@@ -27,8 +27,8 @@ sdp_factories = JuMP.OptimizerFactory[]
 mos && push!(sdp_factories, with_optimizer(MathOptInterfaceMosek.MosekOptimizer, QUIET=true))
 # Currently, need to create a file param.csdp to avoid printing, see https://github.com/JuliaOpt/CSDP.jl/issues/15
 csd && push!(sdp_factories, with_opt(CSDP, printlevel=0))
-iscsdp(factory) = contains(string(factory.constructor), "CSDP")
+iscsdp(factory) = occursin("CSDP", string(factory.constructor))
 # Need 54000 iterations for sosdemo3 to pass on Linux 64 bits
 # With 55000, sosdemo3 passes for every platform except Windows 64 bits on AppVeyor
 scs && push!(sdp_factories, SCS.SCSOptimizer(eps=1e-6, max_iters=60000, verbose=0))
-isscs(factory) = contains(string(factory.constructor), "SCS")
+isscs(factory) = occursin("SCS", string(factory.constructor))

test/sosdemo10.jl

@@ -24,7 +24,7 @@
 
     Sc = [theta^2-s*(gamma-p) g0+g1; g0+g1 1]
 
-    @SDconstraint m eps * eye(2) ⪯ Sc
+    @SDconstraint m Matrix(eps * I, 2, 2) ⪯ Sc
 
     JuMP.optimize!(m)
 

test/sosdemo5.jl

@@ -32,7 +32,8 @@
 
         # -- Q(x)'s -- : sums of squares
         # Monomial vector: [x1; ... x8]
-        Q = Vector{MatPolynomial{JuMP.VariableRef, monomialtype(x[1]), monovectype(x[1])}}(4)
+        Q = Vector{MatPolynomial{JuMP.VariableRef, monomialtype(x[1]),
+                                 monovectype(x[1])}}(undef, 4)
         @variable m Q[1:4] SOSPoly(Z)
 
         # -- r's -- : constant sum of squares

test/sosdemo6.jl

@@ -9,7 +9,7 @@
     f = 2.5 - 0.5*x[1]*x[2] - 0.5*x[2]*x[3] - 0.5*x[3]*x[4] - 0.5*x[4]*x[5] - 0.5*x[5]*x[1]
 
     # Boolean constraints
-    bc = vec(x).^2 - 1
+    bc = vec(x).^2 .- 1
 
     for (gamma, feasible) in [(3.9, false), (4, true)]
 

test/variable.jl

@@ -1,3 +1,5 @@
+import PolyJuMP
+
 @testset "Creating polynomial with empty MonomialVector" begin
     @polyvar x
     X = emptymonovec(typeof(x^2))