More nlp tests [WIP] (#5)

andrewrosemberg · web-flow · commit 1c9caddbf62a · 2025-12-18T17:19:20.000+01:00
* hs10 oracle tests

* add hs14_oracle

* normalize constraints test

* add more tests nlp

* start checking jacobian and hessian

* fix bug hessian

* fix bg consistency
diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl
@@ -312,7 +312,7 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals:
         s.x[i] = x[f.variables[i].value]
       end
       nnz_oracle = length(s.set.jacobian_structure)
-      s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x)
+      s.set.eval_jacobian(view(vals, (offset + 1):(offset + nnz_oracle)), s.x)
       offset += nnz_oracle
     end
   end
@@ -350,8 +350,8 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs
         index += qcon.nnzg
       end
     end
+    offset = nlp.lincon.nnzj + nlp.quadcon.nnzj
     if nlp.meta.nnln > nlp.quadcon.nquad
-      offset = nlp.lincon.nnzj + nlp.quadcon.nnzj
       ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj)
       MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x)
       offset += nlp.nlcon.nnzj
@@ -360,7 +360,7 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs
           s.x[i] = x[f.variables[i].value]
         end
         nnz_oracle = length(s.set.jacobian_structure)
-        s.set.eval_jacobian(view(values, (offset + 1):(offset + nnz_oracle)), s.x)
+        s.set.eval_jacobian(view(vals, (offset + 1):(offset + nnz_oracle)), s.x)
         offset += nnz_oracle
       end
     end
@@ -568,42 +568,55 @@ function NLPModels.hess_coord!(
 )
   increment!(nlp, :neval_hess)
 
-  # Running index over Hessian nonzeros we've filled so far.
-  index = 0
-
-  # Objective Hessian block
+  # 1. Quadratic objective block (if any)
   if nlp.obj.type == "QUADRATIC"
-    view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals
-    index += nlp.obj.nnzh
+    view(vals, 1:nlp.obj.nnzh) .= obj_weight .* nlp.obj.hessian.vals
   end
-  # Nonlinear objective Hessian block
+
+  # 2. Nonlinear block (objective + JuMP @NLconstraint)
   if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
-    λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
+    # Multipliers for the JuMP nonlinear constraints (not the oracles)
+    λ = view(
+      y,
+      (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon),
+    )
+
+    first_nl = nlp.obj.nnzh + nlp.quadcon.nnzh + 1
+    last_nl  = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh
+
     MOI.eval_hessian_lagrangian(
       nlp.eval,
-      view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)),
+      view(vals, first_nl:last_nl),
       x,
       obj_weight,
       λ,
     )
   end
-  # Quadratic constraint Hessian blocks
+
+  # 3. Quadratic constraint Hessian blocks
   if nlp.quadcon.nquad > 0
+    index = nlp.obj.nnzh
     for i = 1:(nlp.quadcon.nquad)
       qcon = nlp.quadcon.constraints[i]
-      view(vals, (index + 1):(index + qcon.nnzh)) .= y[nlp.meta.nlin + i] .* qcon.A.vals
+      ind = (index + 1):(index + qcon.nnzh)
+      view(vals, ind) .= y[nlp.meta.nlin + i] .* qcon.A.vals
       index += qcon.nnzh
     end
   end
-  # Oracle Hessian blocks
+
+  # 4. Oracle Hessian blocks are appended at the very end
   if !isempty(nlp.oracles_data.oracles)
+    # Multipliers for oracle constraints only
     λ_oracle_all = view(
       y,
       (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1):
       (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + nlp.oracles_data.noracle),
     )
 
     λ_offset = 0
+    # Start after obj + quadcon + nonlinear block
+    index = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh
+
     for (f, s) in nlp.oracles_data.oracles
       # build local x for this oracle
       for i in 1:s.set.input_dimension
@@ -622,7 +635,6 @@ function NLPModels.hess_coord!(
     end
   end
 
-  @assert index == nlp.meta.nnzh
   return vals
 end
 
diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl
@@ -1,4 +1,4 @@
-for problem in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_oracle_problems))
+for problem in [nlp_problems; extra_nlp_oracle_problems]
   @testset "Problem $problem" begin
     nlp_manual = eval(Symbol(problem))()
     problem_f = eval(Symbol(lowercase(problem)))
diff --git a/test/nlp_problems/hs100.jl b/test/nlp_problems/hs100.jl
@@ -22,3 +22,93 @@ function hs100(args...; kwargs...)
 
   return nlp
 end
+
+"HS100 model with the 2 middle constraints as VectorNonlinearOracle"
+function hs100_oracle(args...; kwargs...)
+    model = Model()
+    x0 = [1, 2, 0, 4, 0, 1, 1]
+    @variable(model, x[i = 1:7], start = x0[i])
+
+    # 1st constraint: keep as NLconstraint
+    @NLconstraint(model,
+        127 - 2 * x[1]^2 - 3 * x[2]^4 - x[3] - 4 * x[4]^2 - 5 * x[5] ≥ 0
+    )
+
+    # 2nd constraint as oracle:
+    # Original: 282 - 7x1 - 3x2 - 10x3^2 - x4 + x5 ≥ 0
+    # Canonical: f2(x) = -7x1 - 3x2 - 10x3^2 - x4 + x5, l2 = -282, u2 = +∞
+    set2 = MOI.VectorNonlinearOracle(;
+        dimension = 5,               # inputs: x1, x2, x3, x4, x5
+        l = [-282.0],
+        u = [Inf],
+        eval_f = (ret, xv) -> begin
+            # xv = [x1, x2, x3, x4, x5]
+            ret[1] = -7.0 * xv[1] - 3.0 * xv[2] - 10.0 * xv[3]^2 - xv[4] + xv[5]
+        end,
+        # ∇f2 = [-7, -3, -20*x3, -1, 1]
+        jacobian_structure = [(1, 1), (1, 2), (1, 3), (1, 4), (1, 5)],
+        eval_jacobian = (ret, xv) -> begin
+            ret[1] = -7.0              # ∂f2/∂x1
+            ret[2] = -3.0              # ∂f2/∂x2
+            ret[3] = -20.0 * xv[3]     # ∂f2/∂x3
+            ret[4] = -1.0              # ∂f2/∂x4
+            ret[5] =  1.0              # ∂f2/∂x5
+        end,
+        # Hessian of f2: only (3,3) = -20
+        hessian_lagrangian_structure = [(3, 3)],
+        eval_hessian_lagrangian = (ret, xv, μ) -> begin
+            # Hessian of μ[1] * f2(x)
+            ret[1] = μ[1] * (-20.0)    # (3,3)
+        end,
+    )
+    @constraint(model, [x[1], x[2], x[3], x[4], x[5]] in set2)
+
+    # 3rd constraint as oracle:
+    # Original: -196 + 23x1 + x2^2 + 6x6^2 - 8x7 ≤ 0
+    # Canonical: f3(x) = 23x1 + x2^2 + 6x6^2 - 8x7, l3 = -∞, u3 = 196
+    set3 = MOI.VectorNonlinearOracle(;
+        dimension = 4,               # inputs: x1, x2, x6, x7
+        l = [-Inf],
+        u = [196.0],
+        eval_f = (ret, xv) -> begin
+            # xv = [x1, x2, x6, x7]
+            ret[1] = 23.0 * xv[1] + xv[2]^2 + 6.0 * xv[3]^2 - 8.0 * xv[4]
+        end,
+        # ∇f3 = [23, 2*x2, 12*x6, -8]
+        jacobian_structure = [(1, 1), (1, 2), (1, 3), (1, 4)],
+        eval_jacobian = (ret, xv) -> begin
+            ret[1] = 23.0             # ∂f3/∂x1
+            ret[2] =  2.0 * xv[2]     # ∂f3/∂x2
+            ret[3] = 12.0 * xv[3]     # ∂f3/∂x6
+            ret[4] = -8.0             # ∂f3/∂x7
+        end,
+        # Hessian of f3: (2,2) = 2, (3,3) = 12
+        hessian_lagrangian_structure = [(2, 2), (3, 3)],
+        eval_hessian_lagrangian = (ret, xv, μ) -> begin
+            # Hessian of μ[1] * f3(x)
+            ret[1] = μ[1] *  2.0      # (2,2)
+            ret[2] = μ[1] * 12.0      # (3,3)
+        end,
+    )
+    @constraint(model, [x[1], x[2], x[6], x[7]] in set3)
+
+    # 4th constraint: keep as standard nonlinear
+    @constraint(model,
+        -4 * x[1]^2 - x[2]^2 + 3 * x[1] * x[2] - 2 * x[3]^2 - 5 * x[6] + 11 * x[7] ≥ 0
+    )
+
+    # Objective: same as original
+    @NLobjective(
+        model,
+        Min,
+        (x[1] - 10)^2 +
+        5 * (x[2] - 12)^2 +
+        x[3]^4 +
+        3 * (x[4] - 11)^2 +
+        10 * x[5]^6 +
+        7 * x[6]^2 +
+        x[7]^4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]
+    )
+
+    return model
+end
diff --git a/test/nlp_problems/hs14.jl b/test/nlp_problems/hs14.jl
@@ -66,4 +66,4 @@ function hs14_oracle()
     @constraint(model, x[1] - 2 * x[2] + 1 == 0)
 
     return model
-end
+end
diff --git a/test/nlp_problems/hs61.jl b/test/nlp_problems/hs61.jl
@@ -10,3 +10,76 @@ function hs61(args...; kwargs...)
 
   return nlp
 end
+
+"HS61 model with both constraints as VectorNonlinearOracle"
+function hs61_oracle(args...; kwargs...)
+    model = Model()
+    @variable(model, x[i = 1:3], start = 0)
+
+    @NLobjective(model, Min,
+        4 * x[1]^2 + 2 * x[2]^2 + 2 * x[3]^2 -
+        33 * x[1] + 16 * x[2] - 24 * x[3]
+    )
+
+    # First equality: 3*x1 - 2*x2^2 - 7 == 0
+    # Canonical form for MathOptNLPModel:
+    #   f1(x) = 3*x1 - 2*x2^2
+    #   l1 = u1 = 7
+    set1 = MOI.VectorNonlinearOracle(;
+        dimension = 2,           # inputs: x1, x2
+        l = [7.0],
+        u = [7.0],
+        eval_f = (ret, xv) -> begin
+            # xv = [x1, x2]
+            ret[1] = 3.0 * xv[1] - 2.0 * xv[2]^2
+        end,
+        # ∇f1 = [3, -4*x2]
+        jacobian_structure = [(1, 1), (1, 2)],
+        eval_jacobian = (ret, xv) -> begin
+            ret[1] = 3.0             # d f1 / d x1
+            ret[2] = -4.0 * xv[2]    # d f1 / d x2
+        end,
+        # Hessian of f1:
+        # ∂²f1/∂x1² = 0
+        # ∂²f1/∂x2² = -4
+        # (mixed derivatives are 0)
+        hessian_lagrangian_structure = [(2, 2)],
+        eval_hessian_lagrangian = (ret, xv, μ) -> begin
+            # Hessian of μ[1] * f1(x)
+            ret[1] = μ[1] * (-4.0)   # (2,2)
+        end,
+    )
+
+    # Second equality: 4*x1 - x3^2 - 11 == 0
+    # Canonical form:
+    #   f2(x) = 4*x1 - x3^2
+    #   l2 = u2 = 11
+    set2 = MOI.VectorNonlinearOracle(;
+        dimension = 2,           # inputs: x1, x3
+        l = [11.0],
+        u = [11.0],
+        eval_f = (ret, xv) -> begin
+            # xv = [x1, x3]
+            ret[1] = 4.0 * xv[1] - xv[2]^2
+        end,
+        # ∇f2 = [4, -2*x3]
+        jacobian_structure = [(1, 1), (1, 2)],
+        eval_jacobian = (ret, xv) -> begin
+            ret[1] = 4.0             # d f2 / d x1
+            ret[2] = -2.0 * xv[2]    # d f2 / d x3
+        end,
+        # Hessian of f2:
+        # ∂²f2/∂x1² = 0
+        # ∂²f2/∂x3² = -2
+        hessian_lagrangian_structure = [(2, 2)],
+        eval_hessian_lagrangian = (ret, xv, μ) -> begin
+            ret[1] = μ[1] * (-2.0)   # (2,2)
+        end,
+    )
+
+    # Equality constraints as oracles
+    @constraint(model, [x[1], x[2]] in set1)
+    @constraint(model, [x[1], x[3]] in set2)
+
+    return model
+end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -8,7 +8,7 @@ nls_problems = NLPModelsTest.nls_problems
 extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon", "operatorspb", "nf"]
 extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"]
 
-extra_nlp_oracle_problems = ["hs10_oracle", "hs14_oracle"]
+extra_nlp_oracle_problems = ["hs10_oracle", "hs14_oracle", "hs61_oracle", "hs100_oracle"]
 
 for problem in lowercase.(nlp_problems ∪ extra_nlp_problems)
   include(joinpath("nlp_problems", "$problem.jl"))
diff --git a/test/test_moi_nlp_oracle.jl b/test/test_moi_nlp_oracle.jl
@@ -2,20 +2,43 @@
     prob_no_oracle = replace(prob, "_oracle" => "")
     prob_fn_no_oracle = eval(Symbol(prob_no_oracle))
     prob_fn = eval(Symbol(prob))
-    nlp_no_oracle = MathOptNLPModel(prob_fn_no_oracle(), hessian = true)
-    nlp_with_oracle = MathOptNLPModel(prob_fn(), hessian = true)
+
+    nlp_no_oracle   = MathOptNLPModel(prob_fn_no_oracle(), hessian = true)
+    nlp_with_oracle = MathOptNLPModel(prob_fn(),          hessian = true)
+
     n = nlp_no_oracle.meta.nvar
     m = nlp_no_oracle.meta.ncon
     x = nlp_no_oracle.meta.x0
-    fx_no_oracle = obj(nlp_no_oracle, x)
+
+    # Objective value
+    fx_no_oracle   = obj(nlp_no_oracle, x)
     fx_with_oracle = obj(nlp_with_oracle, x)
     @test isapprox(fx_no_oracle, fx_with_oracle; atol = 1e-8, rtol = 1e-8)
-    ngx_no_oracle = grad(nlp_no_oracle, x)
+
+    # Gradient of objective
+    ngx_no_oracle   = grad(nlp_no_oracle, x)
     ngx_with_oracle = grad(nlp_with_oracle, x)
     @test isapprox(ngx_no_oracle, ngx_with_oracle; atol = 1e-8, rtol = 1e-8)
-    if m > 0
-        ncx_no_oracle = cons(nlp_no_oracle, x)
-        ncx_with_oracle = cons(nlp_with_oracle, x)
-        @test isapprox(ncx_no_oracle, ncx_with_oracle; atol = 1e-8, rtol = 1e-8)
-    end
+
+    # Constraint values (up to ordering)
+    ncx_no_oracle   = cons(nlp_no_oracle, x)
+    ncx_with_oracle = cons(nlp_with_oracle, x)
+    @test isapprox(sort(ncx_no_oracle), sort(ncx_with_oracle);
+                    atol = 1e-8, rtol = 1e-8)
+
+    # Jacobian: compare J'J, which is invariant to row permutations
+    J_no   = jac(nlp_no_oracle, x)
+    J_with = jac(nlp_with_oracle, x)
+
+    G_no   = Matrix(J_no)' * Matrix(J_no)
+    G_with = Matrix(J_with)' * Matrix(J_with)
+
+    @test isapprox(G_no, G_with; atol = 1e-8, rtol = 1e-8)
+
+    # Hessian of the objective: use y = 0 so constraints don't enter
+    λ = zeros(m)
+    H_no   = hess(nlp_no_oracle, x, λ)
+    H_with = hess(nlp_with_oracle, x, λ)
+
+    @test isapprox(Matrix(H_no), Matrix(H_with); atol = 1e-8, rtol = 1e-8)
 end

Original file line number	Diff line number	Diff line change
`@@ -1,4 +1,4 @@`
`1`		`-for problem in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_oracle_problems))`
	`1`	`+for problem in [nlp_problems; extra_nlp_oracle_problems]`
`2`	`2`	`@testset "Problem $problem" begin`
`3`	`3`	`nlp_manual = eval(Symbol(problem))()`
`4`	`4`	`problem_f = eval(Symbol(lowercase(problem)))`