Use square_getindex

Author: blegat

src/Sets/polynomials.jl

@@ -63,9 +63,7 @@ function zero_eliminate(set::ConvexPolySet, I)
     vars = space_variables(set)[I]
     K = findall(mono -> all(var -> iszero(degree(mono, var)), vars),
                 set.p.basis.monomials)
-    M = set.p.Q[K, K]
-    Q = SumOfSquares.SymMatrix([M[i, j] for j in 1:length(K) for i in 1:j],
-                               length(K))
+    Q = SumOfSquares.square_getindex(set.p.Q, K)
     monos = set.p.basis.monomials[K]
     J = setdiff(1:dimension(set), I)
     monos = DynamicPolynomials.MonomialVector(monos.vars[J],

test/Sets.jl

@@ -0,0 +1,18 @@
+using Test
+using DynamicPolynomials
+using SetProg, SetProg.Sets
+
+@testset "Sets" begin
+    @testset "zero_eliminate" begin
+        @polyvar x y z
+        p = SetProg.GramMatrix{Float64}((i, j) -> convert(Float64, i + j),
+                                           monovec([x, y, z]))
+        set = Sets.ConvexPolySet(2, p, nothing)
+        el = Sets.zero_eliminate(set, 1:2)
+        @test el.p.Q == 6ones(1, 1)
+        el = Sets.zero_eliminate(set, 3:3)
+        @test el.p.Q == [2 3; 3 4]
+        el = Sets.zero_eliminate(set, 2:2)
+        @test el.p.Q == [2 4; 4 6]
+    end
+end

test/Tests/Tests.jl

@@ -2,7 +2,7 @@ module Tests
 
 using JuMP
 import GLPK
-const lp_solver = optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.ON)
+const lp_solver = optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.GLP_ON)
 import Polyhedra
 const lib = Polyhedra.DefaultLibrary{Float64}(lp_solver)
 

test/recipes.jl

@@ -58,7 +58,8 @@ end
             recipe_test(circle,
                         [1.0, 0.0, -1.0, 0.0], [0.0, 1.0, 0.0, -1.0])
             recipe_test(Sets.polar(circle),
-                        [-1.0, 1.0, 1.0, -1.0, -1.0], [1.0, 1.0, -1.0, -1.0, 1.0])
+                        [-1.0, -1.0, 1.0,  1.0, -1.0],
+                        [-1.0,  1.0, 1.0, -1.0, -1.0])
         end
         @testset "Scaled circle" begin
             p = SetProg.GramMatrix{Float64}((i, j) -> 2convert(Float64, i == j),
@@ -67,8 +68,8 @@ end
             recipe_test(circle,
                         [1/√2, 0.0, -1/√2, 0.0], [0.0, 1/√2, 0.0, -1/√2])
             recipe_test(Sets.polar(circle),
-                        [-√2, √2,  √2, -√2, -√2],
-                        [ √2, √2, -√2, -√2,  √2])
+                        [-√2, -√2, √2,  √2, -√2],
+                        [-√2,  √2, √2, -√2, -√2])
         end
         @testset "Non-homogeneous Circle" begin
             @testset "Basic" begin
@@ -77,17 +78,17 @@ end
                                                0 0 1], monovec([z, x, y]))
                 shifted_circle = SetProg.perspective_dual_polyset(2, q, SetProg.InteriorPoint(zeros(2)), z, [x, y])
                 recipe_test(shifted_circle,
-                            [-1.0, 1.0, 1.0, -1.0, -1.0],
-                            [1.0, 1.0, -1.0, -1.0, 1.0])
+                            [-1.0, -1.0, 1.0,  1.0, -1.0],
+                            [-1.0,  1.0, 1.0, -1.0, -1.0])
             end
             @testset "Scaled" begin
                 q = SetProg.GramMatrix(Float64[0 0 0
                                                0 2 0
                                                0 0 2], monovec([z, x, y]))
                 shifted_circle = SetProg.perspective_dual_polyset(2, q, SetProg.InteriorPoint(zeros(2)), z, [x, y])
                 recipe_test(shifted_circle,
-                           [-√2, √2,  √2, -√2, -√2],
-                           [ √2, √2, -√2, -√2,  √2])
+							[-√2, -√2, √2,  √2, -√2],
+							[-√2,  √2, √2, -√2, -√2])
             end
             @testset "z-Scaled" begin
                 # z: -1/2 + 1 = 1/2
@@ -96,8 +97,8 @@ end
                                         0   0 1], monovec([z, x, y]))
                 shifted_circle = SetProg.perspective_dual_polyset(2, q, SetProg.InteriorPoint(zeros(2)), z, [x, y])
 				recipe_test(shifted_circle,
-							[-√2, √2,  √2, -√2, -√2],
-							[ √2, √2, -√2, -√2,  √2])
+							[-√2, -√2, √2,  √2, -√2],
+							[-√2,  √2, √2, -√2, -√2])
             end
         end
     end
@@ -129,7 +130,8 @@ end
         β = 0.732050807568877
         recipe_test(
             Sets.polar(set),
-            [-1.0, -α, 0.0, 1.0, 1.0, β, -1.0, -1.0],
-            [α, 1.0, 1.0, 0.0, -β, -1.0, -1.0, α], 8)
+            [-1.0, -1.0, -α, 0.0, 1.0, 1.0, β, -1.0],
+            [-1.0, α, 1.0, 1.0, 0.0, -β, -1.0, -1.0],
+            8)
     end
 end

test/runtests.jl

@@ -3,6 +3,7 @@ const lp_solver = GLPK.Optimizer
 import Polyhedra
 const lib = Polyhedra.DefaultLibrary{Float64}(lp_solver)
 
+include("Sets.jl")
 include("apply.jl")
 include("variables.jl")
 include("L1_heuristic.jl")