From c44a8876daa2dc08892cd30c66f2af43db704738 Mon Sep 17 00:00:00 2001
From: Twan Koolen <koolen.twan@gmail.com>
Date: Wed, 5 Dec 2018 16:23:07 -0500
Subject: [PATCH] Add test for quadratic constraints.

---
 test/REQUIRE  |  1 +
 test/model.jl | 47 ++++++++++++++++++++++++++++++++---------------
 2 files changed, 33 insertions(+), 15 deletions(-)

diff --git a/test/REQUIRE b/test/REQUIRE
index 859ebf9..ca87ca5 100644
--- a/test/REQUIRE
+++ b/test/REQUIRE
@@ -1,3 +1,4 @@
 OSQP
 StaticArrays 0.7
 GLPK
+Gurobi
diff --git a/test/model.jl b/test/model.jl
index 9265539..0ab8d68 100644
--- a/test/model.jl
+++ b/test/model.jl
@@ -7,6 +7,7 @@ using Parametron
 using OSQP
 using OSQP.MathOptInterfaceOSQP: OSQPSettings
 using GLPK
+using Gurobi
 using StaticArrays: SVector
 
 import MathOptInterface
@@ -233,21 +234,6 @@ end
     @test value(model, x) ≈ 1.0 atol=1e-8
 end
 
-glpk_works = false
-try
-    optimizer = GLPK.Optimizer()
-    x = MOI.add_variable(optimizer)
-    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
-    s = MOI.GreaterThan(0.0)
-    ci = MOI.add_constraint(optimizer, f, s)
-    MOI.set(optimizer, MOI.ConstraintFunction(), ci, f)
-    global glpk_works = true
-catch e
-    if !(e isa MOI.UnsupportedAttribute{MOI.ConstraintFunction})
-        rethrow(e)
-    end
-end
-
 @testset "integer basics" begin
     # https://github.com/JuliaOpt/GLPK.jl/issues/58
     @testset "scalar constraint" begin
@@ -376,4 +362,35 @@ end
     @test gv ≈ ggt rtol=0.01
 end
 
+@testset "Quadratic constraints" begin
+    if !haskey(ENV, "CI")
+        rng = MersenneTwister(1)
+        optimizer = Gurobi.Optimizer(OutputFlag=0)
+        model = Model(optimizer)
+
+        direction = Parameter(x -> normalize!(randn!(rng, x)), zeros(2), model)
+        zmax = Parameter{Float64}(() -> rand(rng), model)
+
+        x = Variable(model)
+        y = Variable(model)
+        z = Variable(model)
+        μ = 0.7
+
+        @constraint model x^2 + y^2 <= μ^2 * z^2
+        @constraint model z >= 0
+        @constraint model z <= zmax
+        @objective model Maximize direction ⋅ [x, y]
+
+        for i = 1 : 5
+            solve!(model)
+            @test terminationstatus(model) == MOI.Success
+            @test primalstatus(model) == MOI.FeasiblePoint
+
+            xval, yval, zval = value.(model, (x, y, z))
+            @test zval ≈ zmax() atol=1e-6
+            @test [xval, yval] ⋅ direction() ≈ zval * μ atol=1e-6
+        end
+    end
+end
+
 end # module