diff --git a/Project.toml b/Project.toml
index 446fff9..5566c2d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,12 +14,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 JuMP = "1.25"
 LinearAlgebra = "1.10"
-MathOptInterface = "1.40"
-NLPModels = "0.21.5"
-NLPModelsTest = "0.10.3"
+MathOptInterface = "1.46"
+NLPModels = "0.21.6"
+NLPModelsTest = "0.10.6"
 Percival = "0.7.3"
 Printf = "1.10"
-SolverCore = "0.3"
+SolverCore = "0.3.9"
 SparseArrays = "1.10"
 Test = "1.10"
 julia = "1.10"
diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl
index c85a7c0..ed53b66 100644
--- a/src/moi_nlp_model.jl
+++ b/src/moi_nlp_model.jl
@@ -6,6 +6,7 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}}
   lincon::LinearConstraints
   quadcon::QuadraticConstraints
   nlcon::NonLinearStructure
+  oracles::Oracles
   λ::Vector{Float64}
   hv::Vector{Float64}
   obj::Objective
@@ -34,8 +35,10 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
     parser_MOI(moimodel, index_map, nvar)
 
   nlp_data = _nlp_block(moimodel)
-  nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian)
-  λ = zeros(Float64, nnln)  # Lagrange multipliers for hess_coord! and hprod! without y
+  nlcon = parser_NL(nlp_data, hessian = hessian)
+  oracles = parser_oracles(moimodel)
+  counters = Counters()
+  λ = zeros(Float64, nlcon.nnln)  # Lagrange multipliers for hess_coord! and hprod! without y
   hv = zeros(Float64, nvar)  # workspace for ghjvprod!
 
   if nlp_data.has_objective
@@ -44,12 +47,12 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
     obj = parser_objective_MOI(moimodel, nvar, index_map)
   end
 
-  ncon = nlin + quadcon.nquad + nnln
-  lcon = vcat(lin_lcon, quad_lcon, nl_lcon)
-  ucon = vcat(lin_ucon, quad_ucon, nl_ucon)
-  nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj
-  nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh
-
+  # Total counts
+  ncon = nlin + quadcon.nquad + nlcon.nnln + oracles.ncon
+  lcon = vcat(lin_lcon, quad_lcon, nlcon.nl_lcon, oracles.lcon)
+  ucon = vcat(lin_ucon, quad_ucon, nlcon.nl_ucon, oracles.ucon)
+  nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj + oracles.nnzj
+  nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh + oracles.nnzh
   meta = NLPModelMeta(
     nvar,
     x0 = x0,
@@ -63,13 +66,15 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
     nnzh = nnzh,
     lin = collect(1:nlin),
     lin_nnzj = lincon.nnzj,
-    nln_nnzj = quadcon.nnzj + nlcon.nnzj,
+    nln_nnzj = quadcon.nnzj + nlcon.nnzj + oracles.nnzj,
     minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE,
-    islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0),
+    islp = (obj.type == "LINEAR") && (nlin == ncon),
     name = name,
+    hprod_available = hessian && oracles.hessian_oracles_supported,
+    hess_available = hessian && oracles.hessian_oracles_supported,
   )
 
-  return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, λ, hv, obj, Counters()),
+  return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, oracles, λ, hv, obj, counters),
   index_map
 end
 
@@ -113,36 +118,65 @@ end
 
 function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector)
   increment!(nlp, :neval_cons_nln)
+  offset = 0
   if nlp.quadcon.nquad > 0
+    offset += nlp.quadcon.nquad
     for i = 1:(nlp.quadcon.nquad)
       qcon = nlp.quadcon.constraints[i]
       c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
     end
   end
-  if nlp.meta.nnln > nlp.quadcon.nquad
-    MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(nlp.meta.nnln)), x)
+  if nlp.nlcon.nnln > 0
+    offset += nlp.nlcon.nnln
+    index_nlcon = (nlp.quadcon.nquad + 1):(offset)
+    MOI.eval_constraint(nlp.eval, view(c, index_nlcon), x)
+  end
+  if nlp.oracles.ncon > 0
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      index_oracle = (offset + 1):(offset + s.set.output_dimension)
+      c_oracle = view(c, index_oracle)
+      s.set.eval_f(c_oracle, s.x)
+      offset += s.set.output_dimension
+    end
   end
+  @assert offset == nlp.meta.nnln
   return c
 end
 
 function NLPModels.cons!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector)
   increment!(nlp, :neval_cons)
+  offset = 0
   if nlp.meta.nlin > 0
+    offset += nlp.meta.nlin
     coo_prod!(nlp.lincon.jacobian.rows, nlp.lincon.jacobian.cols, nlp.lincon.jacobian.vals, x, c)
   end
-  if nlp.meta.nnln > 0
-    if nlp.quadcon.nquad > 0
-      for i = 1:(nlp.quadcon.nquad)
-        qcon = nlp.quadcon.constraints[i]
-        c[nlp.meta.nlin + i] =
-          0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
-      end
+  if nlp.quadcon.nquad > 0
+    offset += nlp.quadcon.nquad
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon.constraints[i]
+      c[nlp.meta.nlin + i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
     end
-    if nlp.meta.nnln > nlp.quadcon.nquad
-      index_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)
-      MOI.eval_constraint(nlp.eval, view(c, index_nnln), x)
+  end
+  if nlp.nlcon.nnln > 0
+    offset += nlp.nlcon.nnln
+    index_nlcon = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(offset)
+    MOI.eval_constraint(nlp.eval, view(c, index_nlcon), x)
+  end
+  if nlp.oracles.ncon > 0
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      index_oracle = (offset + 1):(offset + s.set.output_dimension)
+      c_oracle = view(c, index_oracle)
+      s.set.eval_f(c_oracle, s.x)
+      offset += s.set.output_dimension
     end
   end
+  @assert offset == nlp.meta.ncon
   return c
 end
 
@@ -162,22 +196,39 @@ function NLPModels.jac_nln_structure!(
   rows::AbstractVector{<:Integer},
   cols::AbstractVector{<:Integer},
 )
+  offset = 0
   if nlp.quadcon.nquad > 0
-    index = 0
     for i = 1:(nlp.quadcon.nquad)
       # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon
       qcon = nlp.quadcon.constraints[i]
-      ind_quad = (index + 1):(index + qcon.nnzg)
+      ind_quad = (offset + 1):(offset + qcon.nnzg)
       view(rows, ind_quad) .= i
       view(cols, ind_quad) .= qcon.g
-      index += qcon.nnzg
+      offset += qcon.nnzg
     end
   end
-  if nlp.meta.nnln > nlp.quadcon.nquad
-    ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj)
+  @assert offset == nlp.quadcon.nnzj
+  if nlp.nlcon.nnln > 0
+     # non-oracle nonlinear constraints
+    ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj)
     view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows
     view(cols, ind_nnln) .= nlp.nlcon.jac_cols
+    offset += nlp.nlcon.nnzj
+  end
+  if nlp.oracles.ncon > 0
+    # oracle nonlinear constraints
+    row_offset = nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      for (r, c) in s.set.jacobian_structure
+        offset += 1
+        rows[offset] = row_offset + r
+        cols[offset] = f.variables[c].value
+      end
+      row_offset += s.set.output_dimension
+    end
+    @assert row_offset == nlp.meta.nnln
   end
+  @assert offset == nlp.meta.nln_nnzj
   return rows, cols
 end
 
@@ -186,29 +237,46 @@ function NLPModels.jac_structure!(
   rows::AbstractVector{<:Integer},
   cols::AbstractVector{<:Integer},
 )
+  offset = 0
   if nlp.meta.nlin > 0
     index_lin = 1:(nlp.lincon.nnzj)
     view(rows, index_lin) .= nlp.lincon.jacobian.rows
     view(cols, index_lin) .= nlp.lincon.jacobian.cols
+    offset += nlp.lincon.nnzj
   end
-  if nlp.meta.nnln > 0
-    if nlp.quadcon.nquad > 0
-      index = 0
-      for i = 1:(nlp.quadcon.nquad)
-        # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon
-        qcon = nlp.quadcon.constraints[i]
-        ind_quad = (nlp.lincon.nnzj + index + 1):(nlp.lincon.nnzj + index + qcon.nnzg)
-        view(rows, ind_quad) .= nlp.meta.nlin .+ i
-        view(cols, ind_quad) .= qcon.g
-        index += qcon.nnzg
-      end
+  if nlp.quadcon.nquad > 0
+    for i = 1:(nlp.quadcon.nquad)
+      # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon
+      qcon = nlp.quadcon.constraints[i]
+      ind_quad = (offset + 1):(offset + qcon.nnzg)
+      view(rows, ind_quad) .= nlp.meta.nlin .+ i
+      view(cols, ind_quad) .= qcon.g
+      offset += qcon.nnzg
     end
-    if nlp.meta.nnln > nlp.quadcon.nquad
-      ind_nnln = (nlp.lincon.nnzj + nlp.quadcon.nnzj + 1):(nlp.meta.nnzj)
-      view(rows, ind_nnln) .= nlp.meta.nlin .+ nlp.quadcon.nquad .+ nlp.nlcon.jac_rows
-      view(cols, ind_nnln) .= nlp.nlcon.jac_cols
+  end
+  @assert offset == nlp.lincon.nnzj + nlp.quadcon.nnzj
+  if nlp.nlcon.nnln > 0
+    # non-oracle nonlinear constraints
+    offset = nlp.lincon.nnzj + nlp.quadcon.nnzj
+    ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj)
+    view(rows, ind_nnln) .= nlp.meta.nlin .+ nlp.quadcon.nquad .+ nlp.nlcon.jac_rows
+    view(cols, ind_nnln) .= nlp.nlcon.jac_cols
+    offset += nlp.nlcon.nnzj
+  end
+  if nlp.oracles.ncon > 0
+    # oracle nonlinear constraints
+    row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      for (r, c) in s.set.jacobian_structure
+        offset += 1
+        rows[offset] = row_offset + r
+        cols[offset] = f.variables[c].value
+      end
+      row_offset += s.set.output_dimension
     end
+    @assert row_offset == nlp.meta.ncon
   end
+  @assert offset == nlp.meta.nnzj
   return rows, cols
 end
 
@@ -221,73 +289,98 @@ end
 
 function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector)
   increment!(nlp, :neval_jac_nln)
+  offset = 0
   if nlp.quadcon.nquad > 0
-    index = 0
     ind_quad = 1:(nlp.quadcon.nnzj)
     view(vals, ind_quad) .= 0.0
     for i = 1:(nlp.quadcon.nquad)
       qcon = nlp.quadcon.constraints[i]
       for (j, ind) in enumerate(qcon.b.nzind)
         k = qcon.dg[ind]
-        vals[index + k] += qcon.b.nzval[j]
+        vals[offset + k] += qcon.b.nzval[j]
       end
       for j = 1:(qcon.nnzh)
         row = qcon.A.rows[j]
         col = qcon.A.cols[j]
         val = qcon.A.vals[j]
         k1 = qcon.dg[row]
-        vals[index + k1] += val * x[col]
+        vals[offset + k1] += val * x[col]
         if row != col
           k2 = qcon.dg[col]
-          vals[index + k2] += val * x[row]
+          vals[offset + k2] += val * x[row]
         end
       end
-      index += qcon.nnzg
+      offset += qcon.nnzg
     end
   end
-  if nlp.meta.nnln > nlp.quadcon.nquad
-    ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj)
+  if nlp.nlcon.nnln > 0
+    ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj)
     MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x)
+    offset += nlp.nlcon.nnzj
+  end
+  if nlp.oracles.ncon > 0
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      nnzj_oracle = length(s.set.jacobian_structure)
+      ind_oracle = (offset + 1):(offset + nnzj_oracle)
+      s.set.eval_jacobian(view(vals, ind_oracle), s.x)
+      offset += nnzj_oracle
+    end
   end
+  @assert offset == nlp.meta.nln_nnzj
   return vals
 end
 
 function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector)
   increment!(nlp, :neval_jac)
+  offset = 0
   if nlp.meta.nlin > 0
     index_lin = 1:(nlp.lincon.nnzj)
     view(vals, index_lin) .= nlp.lincon.jacobian.vals
+    offset += nlp.lincon.nnzj
   end
-  if nlp.meta.nnln > 0
-    if nlp.quadcon.nquad > 0
-      ind_quad = (nlp.lincon.nnzj + 1):(nlp.lincon.nnzj + nlp.quadcon.nnzj)
-      view(vals, ind_quad) .= 0.0
-      index = nlp.lincon.nnzj
-      for i = 1:(nlp.quadcon.nquad)
-        qcon = nlp.quadcon.constraints[i]
-        for (j, ind) in enumerate(qcon.b.nzind)
-          k = qcon.dg[ind]
-          vals[index + k] += qcon.b.nzval[j]
-        end
-        for j = 1:(qcon.nnzh)
-          row = qcon.A.rows[j]
-          col = qcon.A.cols[j]
-          val = qcon.A.vals[j]
-          k1 = qcon.dg[row]
-          vals[index + k1] += val * x[col]
-          if row != col
-            k2 = qcon.dg[col]
-            vals[index + k2] += val * x[row]
-          end
+  if nlp.quadcon.nquad > 0
+    ind_quad = (nlp.lincon.nnzj + 1):(nlp.lincon.nnzj + nlp.quadcon.nnzj)
+    view(vals, ind_quad) .= 0.0
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon.constraints[i]
+      for (j, ind) in enumerate(qcon.b.nzind)
+        k = qcon.dg[ind]
+        vals[offset + k] += qcon.b.nzval[j]
+      end
+      for j = 1:(qcon.nnzh)
+        row = qcon.A.rows[j]
+        col = qcon.A.cols[j]
+        val = qcon.A.vals[j]
+        k1 = qcon.dg[row]
+        vals[offset + k1] += val * x[col]
+        if row != col
+          k2 = qcon.dg[col]
+          vals[offset + k2] += val * x[row]
         end
-        index += qcon.nnzg
       end
+      offset += qcon.nnzg
     end
-    if nlp.meta.nnln > nlp.quadcon.nquad
-      ind_nnln = (nlp.lincon.nnzj + nlp.quadcon.nnzj + 1):(nlp.meta.nnzj)
-      MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x)
+  end
+  if nlp.nlcon.nnln > 0
+    ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj)
+    MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x)
+    offset += nlp.nlcon.nnzj
+  end
+  if nlp.oracles.ncon > 0
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      nnzj_oracle = length(s.set.jacobian_structure)
+      ind_oracle = (offset + 1):(offset + nnzj_oracle)
+      s.set.eval_jacobian(view(vals, ind_oracle), s.x)
+      offset += nnzj_oracle
     end
   end
+  @assert offset == nlp.meta.nnzj
   return vals
 end
 
@@ -323,10 +416,30 @@ function NLPModels.jprod_nln!(
       Jv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
     end
   end
-  if nlp.meta.nnln > nlp.quadcon.nquad
-    ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
+  if nlp.nlcon.nnln > 0
+    ind_nnln = (nlp.quadcon.nquad + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln)
     MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v)
   end
+  if nlp.oracles.ncon > 0
+    for i in (nlp.quadcon.nquad + nlp.nlcon.nnln + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln + nlp.oracles.ncon)
+      Jv[i] = 0
+    end
+    row_offset = nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      s.set.eval_jacobian(s.nzJ, s.x)
+      k = 0
+      for (r, c) in s.set.jacobian_structure
+        k += 1
+        row = row_offset + r
+        col = f.variables[c].value
+        Jv[row] += s.nzJ[k] * v[col]
+      end
+      row_offset += s.set.output_dimension
+    end
+  end
   return Jv
 end
 
@@ -358,10 +471,30 @@ function NLPModels.jprod!(
         coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
     end
   end
-  if nlp.meta.nnln > nlp.quadcon.nquad
-    ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)
+  if nlp.nlcon.nnln > 0
+    ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln)
     MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v)
   end
+  if nlp.oracles.ncon > 0
+    for i in (nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1):(nlp.meta.ncon)
+      Jv[i] = 0
+    end
+    row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      s.set.eval_jacobian(s.nzJ, s.x)
+      k = 0
+      for (r, c) in s.set.jacobian_structure
+        k += 1
+        row = row_offset + r
+        col = f.variables[c].value
+        Jv[row] += s.nzJ[k] * v[col]
+      end
+      row_offset += s.set.output_dimension
+    end
+  end
   return Jv
 end
 
@@ -390,11 +523,11 @@ function NLPModels.jtprod_nln!(
   Jtv::AbstractVector,
 )
   increment!(nlp, :neval_jtprod_nln)
-  if nlp.meta.nnln > nlp.quadcon.nquad
-    ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
+  if nlp.nlcon.nnln > 0
+    ind_nnln = (nlp.quadcon.nquad + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln)
     MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln))
   end
-  (nlp.quadcon.nquad == nlp.meta.nnln) && (Jtv .= 0.0)
+  (nlp.nlcon.nnln == 0) && (Jtv .= 0.0)
   if nlp.quadcon.nquad > 0
     for i = 1:(nlp.quadcon.nquad)
       # Jtv += v[i] * (Aᵢ * x + bᵢ)
@@ -403,6 +536,23 @@ function NLPModels.jtprod_nln!(
       Jtv .+= v[i] .* qcon.b
     end
   end
+  if nlp.oracles.ncon > 0
+    row_offset = nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      s.set.eval_jacobian(s.nzJ, s.x)
+      k = 0
+      for (r, c) in s.set.jacobian_structure
+        k += 1
+        row = row_offset + r
+        col = f.variables[c].value
+        Jtv[col] += s.nzJ[k] * v[row]
+      end
+      row_offset += s.set.output_dimension
+    end
+  end
   return Jtv
 end
 
@@ -413,11 +563,11 @@ function NLPModels.jtprod!(
   Jtv::AbstractVector,
 )
   increment!(nlp, :neval_jtprod)
-  if nlp.meta.nnln > nlp.quadcon.nquad
-    ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)
+  if nlp.nlcon.nnln > 0
+    ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln)
     MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln))
   end
-  (nlp.quadcon.nquad == nlp.meta.nnln) && (Jtv .= 0.0)
+  (nlp.nlcon.nnln == 0) && (Jtv .= 0.0)
   if nlp.meta.nlin > 0
     transpose = true
     coo_unsym_add_mul!(
@@ -437,6 +587,23 @@ function NLPModels.jtprod!(
       Jtv .+= v[nlp.meta.nlin + i] .* qcon.b
     end
   end
+  if nlp.oracles.ncon > 0
+    row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      s.set.eval_jacobian(s.nzJ, s.x)
+      k = 0
+      for (r, c) in s.set.jacobian_structure
+        k += 1
+        row = row_offset + r
+        col = f.variables[c].value
+        Jtv[col] += s.nzJ[k] * v[row]
+      end
+      row_offset += s.set.output_dimension
+    end
+  end
   return Jtv
 end
 
@@ -445,16 +612,13 @@ function NLPModels.hess_structure!(
   rows::AbstractVector{<:Integer},
   cols::AbstractVector{<:Integer},
 )
+  nlp.meta.hess_available || error("The function hess_structure! is not supported by this MathOptNLPModel.")
   if nlp.obj.type == "QUADRATIC"
     view(rows, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.rows
     view(cols, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.cols
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
-    view(rows, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_rows
-    view(cols, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_cols
-  end
+  index = nlp.obj.nnzh
   if nlp.quadcon.nquad > 0
-    index = nlp.obj.nnzh
     for i = 1:(nlp.quadcon.nquad)
       qcon = nlp.quadcon.constraints[i]
       view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows
@@ -462,6 +626,20 @@ function NLPModels.hess_structure!(
       index += qcon.nnzh
     end
   end
+  if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0)
+    view(rows, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_rows
+    view(cols, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_cols
+    index += nlp.nlcon.nnzh
+  end
+  if nlp.oracles.ncon > 0
+    for (f, s) in nlp.oracles.oracles
+      for (i, j) in s.set.hessian_lagrangian_structure
+        index += 1
+        rows[index] = f.variables[i].value
+        cols[index] = f.variables[j].value
+      end
+    end
+  end
   return rows, cols
 end
 
@@ -472,28 +650,66 @@ function NLPModels.hess_coord!(
   vals::AbstractVector;
   obj_weight::Float64 = 1.0,
 )
+  nlp.meta.hess_available || error("The function hess_coord! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_hess)
+
+  # 1. Quadratic objective block (if any)
   if nlp.obj.type == "QUADRATIC"
-    view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals
+    view(vals, 1:nlp.obj.nnzh) .= obj_weight .* nlp.obj.hessian.vals
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
-    λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
+
+  # 2. Nonlinear block (objective + JuMP @NLconstraint)
+  if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0)
+    # Multipliers for the JuMP nonlinear constraints (not the oracles)
+    ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln)
+    λ_nnln = view(y, ind_nnln)
+    ind_nnzh = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(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, ind_nnzh),
       x,
       obj_weight,
-      λ,
+      λ_nnln,
     )
   end
+
+  # 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
+
+  # 4. Oracle Hessian blocks are appended at the very end
+  if nlp.oracles.ncon > 0
+    offset_λ = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
+    offset_nnzh = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh
+
+    for (f, s) in nlp.oracles.oracles
+      # build local x for this oracle
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+
+      ncon_oracle = s.set.output_dimension
+      index_λ = (offset_λ + 1):(offset_λ + ncon_oracle)
+      λ_oracles = view(y, index_λ)
+
+      nnzh_oracle = length(s.set.hessian_lagrangian_structure)
+      index_nnzh = (offset_nnzh + 1):(offset_nnzh + nnzh_oracle)
+      s.set.eval_hessian_lagrangian(view(vals, index_nnzh), s.x, λ_oracles)
+
+      offset_λ += ncon_oracle
+      offset_nnzh += nnzh_oracle
+    end
+    @assert offset_λ == nlp.meta.ncon
+    @assert offset_nnzh == nlp.meta.nnzh
+  end
+
   return vals
 end
 
@@ -503,6 +719,7 @@ function NLPModels.hess_coord!(
   vals::AbstractVector;
   obj_weight::Float64 = 1.0,
 )
+  nlp.meta.hess_available || error("The function hess_coord! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_hess)
   if nlp.obj.type == "LINEAR"
     vals .= 0.0
@@ -513,7 +730,8 @@ function NLPModels.hess_coord!(
   end
   if nlp.obj.type == "NONLINEAR"
     view(vals, 1:(nlp.obj.nnzh + nlp.quadcon.nnzh)) .= 0.0
-    ind_nnln = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)
+    view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh + 1):(nlp.meta.nnzh)) .= 0.0
+    ind_nnln = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh)
     MOI.eval_hessian_lagrangian(nlp.eval, view(vals, ind_nnln), x, obj_weight, nlp.λ)
   end
   return vals
@@ -525,9 +743,12 @@ function NLPModels.jth_hess_coord!(
   j::Integer,
   vals::AbstractVector,
 )
+  nlp.meta.hess_available || error("The function jth_hess_coord! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_jhess)
   @rangecheck 1 nlp.meta.ncon j
   vals .= 0.0
+
+  # Quadratic constraints
   if nlp.meta.nlin + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad
     index = nlp.obj.nnzh
     for i = 1:(nlp.quadcon.nquad)
@@ -537,17 +758,53 @@ function NLPModels.jth_hess_coord!(
       end
       index += qcon.nnzh
     end
-  elseif nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.ncon
+  end
+
+  # Non-oracle nonlinear constraints
+  if nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
     nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 1.0
+    index_nnzh = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)
     MOI.eval_hessian_lagrangian(
       nlp.eval,
-      view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)),
+      view(vals, index_nnzh),
       x,
       0.0,
       nlp.λ,
     )
     nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 0.0
   end
+
+  # Oracle nonlinear constraints
+  if nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1 ≤ j ≤ nlp.meta.ncon
+    offset_λ = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
+    offset_nnzh = nlp.obj.nnzh + nlp.quadcon.nnzh + nlp.nlcon.nnzh
+    find_oracle = false
+    for (f, s) in nlp.oracles.oracles
+      if !find_oracle
+        ncon_oracle = s.set.output_dimension
+        nnzh_oracle = length(s.set.hessian_lagrangian_structure)
+        if offset_λ + 1 ≤ j ≤ offset_λ + s.set.output_dimension
+          find_oracle = true
+          # build local x
+          for i in 1:s.set.input_dimension
+            s.x[i] = x[f.variables[i].value]
+          end
+          # We could preallocate the vector in the structure nlp
+          λ_oracle = zeros(Float64, s.set.output_dimension)
+          λ_oracle[j - offset_λ] = 1.0
+          index_nnzh = (offset_nnzh + 1):(offset_nnzh + nnzh_oracle)
+          s.set.eval_hessian_lagrangian(
+            view(vals, index_nnzh),
+            s.x,
+            λ_oracle,
+          )
+        else
+          offset_λ += ncon_oracle
+          offset_nnzh += nnzh_oracle
+        end
+      end
+    end
+  end
   return vals
 end
 
@@ -559,16 +816,18 @@ function NLPModels.hprod!(
   hv::AbstractVector;
   obj_weight::Float64 = 1.0,
 )
+  nlp.meta.hprod_available || error("The function hprod! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_hprod)
   if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0)
     hv .= 0.0
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
-    λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
+  if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0)
+    ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln)
+    λ = view(y, ind_nnln)
     MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, λ)
   end
   if nlp.obj.type == "QUADRATIC"
-    (nlp.meta.nnln == nlp.quadcon.nquad) && (hv .= 0.0)
+    (nlp.nlcon.nnln == 0) && (hv .= 0.0)
     coo_sym_add_mul!(
       nlp.obj.hessian.rows,
       nlp.obj.hessian.cols,
@@ -579,12 +838,35 @@ function NLPModels.hprod!(
     )
   end
   if nlp.quadcon.nquad > 0
-    (nlp.obj.type == "LINEAR") && (hv .= 0.0)
+    (nlp.obj.type == "LINEAR") && (nlp.nlcon.nnln == 0) && (hv .= 0.0)
     for i = 1:(nlp.quadcon.nquad)
       qcon = nlp.quadcon.constraints[i]
       coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nlp.meta.nlin + i])
     end
   end
+  if nlp.oracles.ncon > 0
+    (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == nlp.oracles.ncon) && (hv .= 0.0)
+    offset_y = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      for i in 1:s.set.input_dimension
+        s.x[i] = x[f.variables[i].value]
+      end
+      ind_y = (offset_y + 1):(offset_y + s.set.output_dimension)
+      y_oracle = view(y, ind_y)
+      s.set.eval_hessian_lagrangian(s.nzH, s.x, y_oracle)
+      k = 0
+      for (r, c) in s.set.hessian_lagrangian_structure
+        k += 1
+        row = f.variables[r].value
+        col = f.variables[c].value
+        hv[row] += s.nzH[k] * v[col]
+        if row != col
+          hv[col] += s.nzH[k] * v[row]
+        end
+      end
+      offset_y += s.set.output_dimension
+    end
+  end
   return hv
 end
 
@@ -595,6 +877,7 @@ function NLPModels.hprod!(
   hv::AbstractVector;
   obj_weight::Float64 = 1.0,
 )
+  nlp.meta.hprod_available || error("The function hprod! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_hprod)
   if nlp.obj.type == "LINEAR"
     hv .= 0.0
@@ -623,16 +906,51 @@ function NLPModels.jth_hprod!(
   j::Integer,
   hv::AbstractVector,
 )
+  nlp.meta.hprod_available || error("The function jth_hprod! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_jhprod)
   @rangecheck 1 nlp.meta.ncon j
   hv .= 0.0
   if nlp.meta.nlin + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad
     qcon = nlp.quadcon.constraints[j - nlp.meta.nlin]
     coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, 1.0)
-  elseif nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.ncon
+  elseif nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
     nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 1.0
     MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, 0.0, nlp.λ)
     nlp.λ[j - nlp.meta.nlin - nlp.quadcon.nquad] = 0.0
+  elseif nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln + 1 ≤ j ≤ nlp.meta.ncon
+    find_oracle = false
+    offset_λ = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
+    for (f, s) in nlp.oracles.oracles
+      if !find_oracle
+        if offset_λ + 1 ≤ j ≤ offset_λ + s.set.output_dimension
+          find_oracle = true
+          # build local x
+          for i in 1:s.set.input_dimension
+            s.x[i] = x[f.variables[i].value]
+          end
+          # We could preallocate the vector in the structure nlp
+          λ_oracle = zeros(Float64, s.set.output_dimension)
+          λ_oracle[j - offset_λ] = 1.0
+          s.set.eval_hessian_lagrangian(
+            s.nzH,
+            s.x,
+            λ_oracle,
+          )
+          k = 0
+          for (r, c) in s.set.hessian_lagrangian_structure
+            k += 1
+            row = f.variables[r].value
+            col = f.variables[c].value
+            hv[row] += s.nzH[k] * v[col]
+            if row != col
+              hv[col] += s.nzH[k] * v[row]
+            end
+          end
+        else
+          offset_λ += ncon_oracle
+        end
+      end
+    end
   end
   return hv
 end
@@ -644,6 +962,7 @@ function NLPModels.ghjvprod!(
   v::AbstractVector,
   ghv::AbstractVector,
 )
+  nlp.meta.hprod_available || error("The function ghjvprod! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_hprod)
   ghv .= 0.0
   for i = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad)
diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl
index 8eeccf9..e70d701 100644
--- a/src/moi_nls_model.jl
+++ b/src/moi_nls_model.jl
@@ -36,16 +36,16 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri
     parser_MOI(moimodel, index_map, nvar)
 
   nlp_data = _nlp_block(moimodel)
-  nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian)
-  λ = zeros(nnln)  # Lagrange multipliers for hess_coord! and hprod! without y
+  nlcon = parser_NL(nlp_data, hessian = hessian)
+  λ = zeros(nlcon.nnln)  # Lagrange multipliers for hess_coord! and hprod! without y
 
   nequ = nlinequ + nnlnequ
   Fnnzj = linequ.nnzj + nlequ.nnzj
   Fnnzh = nlequ.nnzh
 
-  ncon = nlin + quadcon.nquad + nnln
-  lcon = vcat(lin_lcon, quad_lcon, nl_lcon)
-  ucon = vcat(lin_ucon, quad_ucon, nl_ucon)
+  ncon = nlin + quadcon.nquad + nlcon.nnln
+  lcon = vcat(lin_lcon, quad_lcon, nlcon.nl_lcon)
+  ucon = vcat(lin_ucon, quad_ucon, nlcon.nl_ucon)
   cnnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj
   cnnzh = lls.nnzh + quadcon.nnzh + nlcon.nnzh
 
@@ -265,7 +265,7 @@ function NLPModels.cons_nln!(nls::MathOptNLSModel, x::AbstractVector, c::Abstrac
       c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
     end
   end
-  if nls.meta.nnln > nls.quadcon.nquad
+  if nls.nlcon.nnln > 0
     MOI.eval_constraint(nls.ceval, view(c, (nls.quadcon.nquad + 1):(nls.meta.nnln)), x)
   end
   return c
@@ -284,7 +284,7 @@ function NLPModels.cons!(nls::MathOptNLSModel, x::AbstractVector, c::AbstractVec
           0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
       end
     end
-    if nls.meta.nnln > nls.quadcon.nquad
+    if nls.nlcon.nnln > 0
       index_nnln = (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon)
       MOI.eval_constraint(nls.ceval, view(c, index_nnln), x)
     end
@@ -319,7 +319,7 @@ function NLPModels.jac_nln_structure!(
       index += qcon.nnzg
     end
   end
-  if nls.meta.nnln > nls.quadcon.nquad
+  if nls.nlcon.nnln > 0
     ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj)
     view(rows, ind_nnln) .= nls.quadcon.nquad .+ nls.nlcon.jac_rows
     view(cols, ind_nnln) .= nls.nlcon.jac_cols
@@ -349,7 +349,7 @@ function NLPModels.jac_structure!(
         index += qcon.nnzg
       end
     end
-    if nls.meta.nnln > nls.quadcon.nquad
+    if nls.nlcon.nnln > 0
       ind_nnln = (nls.lincon.nnzj + nls.quadcon.nnzj + 1):(nls.meta.nnzj)
       view(rows, ind_nnln) .= nls.meta.nlin .+ nls.quadcon.nquad .+ nls.nlcon.jac_rows
       view(cols, ind_nnln) .= nls.nlcon.jac_cols
@@ -391,7 +391,7 @@ function NLPModels.jac_nln_coord!(nls::MathOptNLSModel, x::AbstractVector, vals:
       index += qcon.nnzg
     end
   end
-  if nls.meta.nnln > nls.quadcon.nquad
+  if nls.nlcon.nnln > 0
     ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj)
     MOI.eval_constraint_jacobian(nls.ceval, view(vals, ind_nnln), x)
   end
@@ -429,7 +429,7 @@ function NLPModels.jac_coord!(nls::MathOptNLSModel, x::AbstractVector, vals::Abs
         index += qcon.nnzg
       end
     end
-    if nls.meta.nnln > nls.quadcon.nquad
+    if nls.nlcon.nnln > 0
       ind_nnln = (nls.lincon.nnzj + nls.quadcon.nnzj + 1):(nls.meta.nnzj)
       MOI.eval_constraint_jacobian(nls.ceval, view(vals, ind_nnln), x)
     end
@@ -469,7 +469,7 @@ function NLPModels.jprod_nln!(
       Jv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
     end
   end
-  if nls.meta.nnln > nls.quadcon.nquad
+  if nls.nlcon.nnln > 0
     ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln)
     MOI.eval_constraint_jacobian_product(nls.ceval, view(Jv, ind_nnln), x, v)
   end
@@ -504,7 +504,7 @@ function NLPModels.jprod!(
         coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
     end
   end
-  if nls.meta.nnln > nls.quadcon.nquad
+  if nls.nlcon.nnln > 0
     ind_nnln = (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon)
     MOI.eval_constraint_jacobian_product(nls.ceval, view(Jv, ind_nnln), x, v)
   end
@@ -536,7 +536,7 @@ function NLPModels.jtprod_nln!(
   Jtv::AbstractVector,
 )
   increment!(nls, :neval_jtprod_nln)
-  if nls.meta.nnln > nls.quadcon.nquad
+  if nls.nlcon.nnln > 0
     ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln)
     MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, view(v, ind_nnln))
   end
@@ -559,7 +559,7 @@ function NLPModels.jtprod!(
   Jtv::AbstractVector,
 )
   increment!(nls, :neval_jtprod)
-  if nls.meta.nnln > nls.quadcon.nquad
+  if nls.nlcon.nnln > 0
     ind_nnln = (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon)
     MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, view(v, ind_nnln))
   end
@@ -595,7 +595,7 @@ function NLPModels.hess_structure!(
     view(rows, 1:(nls.lls.nnzh)) .= nls.lls.hessian.rows
     view(cols, 1:(nls.lls.nnzh)) .= nls.lls.hessian.cols
   end
-  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad)
+  if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0)
     view(rows, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_rows
     view(cols, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_cols
   end
@@ -622,7 +622,7 @@ function NLPModels.hess_coord!(
   if nls.nls_meta.nlin > 0
     view(vals, 1:(nls.lls.nnzh)) .= obj_weight .* nls.lls.hessian.vals
   end
-  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad)
+  if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0)
     λ = view(y, (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon))
     MOI.eval_hessian_lagrangian(
       nls.ceval,
@@ -658,7 +658,7 @@ function NLPModels.hess_coord!(
     ind_nnln = (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)
     MOI.eval_hessian_lagrangian(nls.ceval, view(vals, ind_nnln), x, obj_weight, nls.λ)
   else
-    if nls.meta.nnln > nls.quadcon.nquad
+    if nls.nlcon.nnln > 0
       view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= 0.0
     end
   end
@@ -674,7 +674,7 @@ function NLPModels.hprod!(
   obj_weight::Float64 = 1.0,
 )
   increment!(nls, :neval_hprod)
-  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad)
+  if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0)
     λ = view(y, (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon))
     MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, λ)
   end
diff --git a/src/utils.jl b/src/utils.jl
index a98b54b..5ff1ef8 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -24,6 +24,23 @@ const SNF = MOI.ScalarNonlinearFunction
 const VNF = MOI.VectorNonlinearFunction
 const NF = Union{SNF, VNF}
 
+# VectorNonlinearOracle
+const ORACLE = MOI.VectorNonlinearOracle{Float64}  # VectorNonlinearOracle{Float64}(input_dimension, output_dimension, l, u, eval_f, jacobian_structure, eval_jacobian, hessian_lagrangian_structure, eval_hessian_lagrangian)
+
+# Cache of VectorNonlinearOracle
+mutable struct _VectorNonlinearOracleCache
+    set::MOI.VectorNonlinearOracle{Float64}
+    x::Vector{Float64}
+    nzJ::Vector{Float64}
+    nzH::Vector{Float64}
+
+    function _VectorNonlinearOracleCache(set::MOI.VectorNonlinearOracle{Float64})
+        nnzj = length(set.jacobian_structure)
+        nnzh = length(set.hessian_lagrangian_structure)
+        return new(set, zeros(set.input_dimension), zeros(nnzj), zeros(nnzh))
+    end
+end
+
 # AffLinSets and VecLinSets
 const ALS = Union{
   MOI.EqualTo{Float64},
@@ -74,7 +91,24 @@ mutable struct QuadraticConstraints
   nnzh::Int
 end
 
+"""
+    NonLinearStructure
+
+Structure containing Jacobian and Hessian structures of nonlinear constraints:
+- nnln: number of nonlinear constraints
+- nl_lcon: lower bounds of nonlinear constraints
+- nl_ucon: upper bounds of nonlinear constraints
+- jac_rows: row indices of the Jacobian in Coordinate format (COO) format
+- jac_cols: column indices of the Jacobian in COO format
+- nnzj: number of non-zero entries in the Jacobian
+- hess_rows: row indices of the Hessian in COO format
+- hess_cols: column indices of the Hessian in COO format
+- nnzh: number of non-zero entries in the Hessian
+"""
 mutable struct NonLinearStructure
+  nnln::Int
+  nl_lcon::Vector{Float64}
+  nl_ucon::Vector{Float64}
   jac_rows::Vector{Int}
   jac_cols::Vector{Int}
   nnzj::Int
@@ -97,6 +131,30 @@ mutable struct Objective
   nnzh::Int
 end
 
+"""
+    Oracles
+
+Structure containing nonlinear oracles data:
+- oracles: vector of tuples (MOI.VectorOfVariables, _VectorNonlinearOracleCache)
+- ncon: number of scalar constraints represented by all oracles
+- lcon: lower bounds of oracle constraints
+- ucon: upper bounds of oracle constraints
+- nnzj: number of non-zero entries in the Jacobian of all oracles
+- nnzh: number of non-zero entries in the Hessian of all oracles
+- nzJ: buffer to store the nonzeros of the Jacobian for all oracles (needed for the functions jprod and jtprod)
+- nzH: buffer to store the nonzeros of the Hessian for all oracles (needed for the function hprod)
+- hessian_oracles_supported: support of the Hessian for all oracles
+"""
+mutable struct Oracles
+  oracles::Vector{Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}}
+  ncon::Int
+  lcon::Vector{Float64}
+  ucon::Vector{Float64}
+  nnzj::Int
+  nnzh::Int
+  hessian_oracles_supported::Bool
+end
+
 """
     replace!(ex, x)
 
@@ -407,6 +465,10 @@ function parser_MOI(moimodel, index_map, nvar)
 
   contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent())
   for (F, S) in contypes
+    # Ignore VectorNonlinearOracle here, we'll parse it separately
+    if F == MOI.VectorOfVariables && S <: MOI.VectorNonlinearOracle{Float64}
+        continue
+    end
     (F == VNF) && error(
       "The function $F is not supported. Please use `.<=`, `.==`, and `.>=` in your constraints to ensure compatibility with ScalarNonlinearFunction.",
     )
@@ -444,7 +506,6 @@ function parser_MOI(moimodel, index_map, nvar)
   coo = COO(linrows, lincols, linvals)
   lin_nnzj = length(linvals)
   lincon = LinearConstraints(coo, lin_nnzj)
-
   quad_nnzj = 0
   quad_nnzh = 0
   for i = 1:nquad
@@ -529,6 +590,9 @@ end
     parser_NL(nlp_data; hessian)
 
 Parse nonlinear constraints of an `nlp_data`.
+
+Returns:
+- nlcon: NonLinearStructure containing Jacobian and Hessian structures
 """
 function parser_NL(nlp_data; hessian::Bool = true)
   nnln = length(nlp_data.constraint_bounds)
@@ -546,9 +610,54 @@ function parser_NL(nlp_data; hessian::Bool = true)
   hess_rows = hessian ? getindex.(hess, 1) : Int[]
   hess_cols = hessian ? getindex.(hess, 2) : Int[]
   nnzh = length(hess)
-  nlcon = NonLinearStructure(jac_rows, jac_cols, nnzj, hess_rows, hess_cols, nnzh)
+  nlcon = NonLinearStructure(nnln, nl_lcon, nl_ucon, jac_rows, jac_cols, nnzj, hess_rows, hess_cols, nnzh)
+
+  return nlcon
+end
+
+"""
+    oracles = parser_oracles(moimodel)
+
+Parse nonlinear oracles of a `MOI.ModelLike`.
+"""
+function parser_oracles(moimodel)
+    hessian_oracles_supported = true
+    oracles = Tuple{MOI.VectorOfVariables,_VectorNonlinearOracleCache}[]
+    lcon = Float64[]
+    ucon = Float64[]
+
+    # We know this pair exists from ListOfConstraintTypesPresent
+    for ci in MOI.get(
+        moimodel,
+        MOI.ListOfConstraintIndices{MOI.VectorOfVariables, MOI.VectorNonlinearOracle{Float64}}(),
+    )
+        f   = MOI.get(moimodel, MOI.ConstraintFunction(), ci)  # ::MOI.VectorOfVariables
+        set = MOI.get(moimodel, MOI.ConstraintSet(), ci)       # ::MOI.VectorNonlinearOracle{Float64}
+
+        cache = _VectorNonlinearOracleCache(set)
+        push!(oracles, (f, cache))
+
+        # Bounds: MOI.VectorNonlinearOracle stores them internally (l, u)
+        append!(lcon, set.l)
+        append!(ucon, set.u)
+
+        # Support for the Hessian
+        hessian_oracles_supported = hessian_oracles_supported && !isnothing(set.eval_hessian_lagrangian)
+    end
+
+    # Number of scalar constraints represented by all oracles
+    ncon = length(lcon)
+
+    # Number of nonzeros for the Jacobian and Hessian
+    nnzj = 0
+    nnzh = 0
+    for (_, cache) in oracles
+        nnzj += length(cache.set.jacobian_structure)
+        # there may or may not be Hessian info
+        nnzh += length(cache.set.hessian_lagrangian_structure)
+    end
 
-  return nnln, nlcon, nl_lcon, nl_ucon
+    return Oracles(oracles, ncon, lcon, ucon, nnzj, nnzh, hessian_oracles_supported)
 end
 
 """
@@ -742,6 +851,7 @@ function parser_nonlinear_expression(cmodel, nvar, F; hessian::Bool = true)
 
   # Nonlinear least squares model
   F_is_array_of_containers = F isa Array{<:AbstractArray}
+  nnlnequ = 0
   if F_is_array_of_containers
     nnlnequ = sum(sum(isa(Fi, NLE) for Fi in FF) for FF in F)
     if nnlnequ > 0
@@ -778,7 +888,7 @@ function parser_nonlinear_expression(cmodel, nvar, F; hessian::Bool = true)
   Fhess_cols = hessian && Feval ≠ nothing ? getindex.(Fhess, 2) : Int[]
   nl_Fnnzh = length(Fhess)
 
-  nlequ = NonLinearStructure(Fjac_rows, Fjac_cols, nl_Fnnzj, Fhess_rows, Fhess_cols, nl_Fnnzh)
+  nlequ = NonLinearStructure(nnlnequ, Float64[], Float64[], Fjac_rows, Fjac_cols, nl_Fnnzj, Fhess_rows, Fhess_cols, nl_Fnnzh)
 
   return Feval, nlequ, nnlnequ
 end
diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl
index b649357..e3739d0 100644
--- a/test/nlp_consistency.jl
+++ b/test/nlp_consistency.jl
@@ -1,17 +1,19 @@
-for problem in nlp_problems
+for problem in [nlp_problems; nlp_oracle_problems]
   @testset "Problem $problem" begin
-    nlp_manual = eval(Symbol(problem))()
-    problem_f = eval(Symbol(lowercase(problem)))
+    nlp_manual = eval(replace(problem, "_oracle" => "") |> uppercase |> Symbol)()
+    problem_f = eval(problem |> lowercase |> Symbol)
     nlp_moi = MathOptNLPModel(problem_f())
     nlps = [nlp_manual; nlp_moi]
-    consistent_nlps(
-      nlps,
-      linear_api = true,
-      test_slack = false,
-      test_counters = false,
-      exclude = [],
-    )
-    coord_memory_nlp(nlp_moi, linear_api = true)
+    @testset "linear_api = $linear_api" for linear_api in (false, true)
+      consistent_nlps(
+        nlps,
+        linear_api = linear_api,
+        test_slack = false,
+        test_counters = false,
+        exclude = [],
+      )
+      coord_memory_nlp(nlp_moi, linear_api = linear_api)
+    end
     view_subarray_nlp(nlp_moi)
   end
 end
diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl
index 3e657e1..200e8e4 100644
--- a/test/nlp_problems/hs10.jl
+++ b/test/nlp_problems/hs10.jl
@@ -12,3 +12,50 @@ function hs10()
 
   return nlp
 end
+
+"Problem 10 in the Hock-Schittkowski suite, but the nonlinear constraint is a VectorNonlinearOracle"
+function hs10_oracle()
+    model = Model()
+
+    @variable(model, x[1:2])
+    set_start_value(x[1], -10)
+    set_start_value(x[2], 10)
+
+    @objective(model, Min, x[1] - x[2])
+
+    # g(x) = -3 x1^2 + 2 x1 x2 - x2^2 + 1 ≥ 0
+    # Bounds: 0 ≤ g(x) ≤ +∞
+    set = MOI.VectorNonlinearOracle(;
+        dimension = 2,                # number of input variables (x1, x2)
+        l = [-1.0],                   # lower bound on g(x)
+        u = [Inf],                    # upper bound on g(x)
+        eval_f = (ret, xv) -> begin
+            # ret[1] = g(x)
+            ret[1] = -3 * xv[1]^2 + 2 * xv[1] * xv[2] - xv[2]^2
+        end,
+        # Jacobian of g(x): ∇g = [∂g/∂x1, ∂g/∂x2]
+        # ∂g/∂x1 = -6 x1 + 2 x2
+        # ∂g/∂x2 =  2 x1 - 2 x2
+        jacobian_structure = [(1, 1), (1, 2)],  # (row, col) for each entry of ret
+        eval_jacobian = (ret, xv) -> begin
+            ret[1] = -6 * xv[1] + 2 * xv[2]     # d g / d x1
+            ret[2] =  2 * xv[1] - 2 * xv[2]     # d g / d x2
+        end,
+        # Hessian of g(x) (constant):
+        # ∂²g/∂x1² = -6
+        # ∂²g/∂x1∂x2 = ∂²g/∂x2∂x1 = 2
+        # ∂²g/∂x2² = -2
+        hessian_lagrangian_structure = [(1, 1), (2, 1), (2, 2)],
+        eval_hessian_lagrangian = (ret, xv, μ) -> begin
+            # Hessian of μ[1] * g(x)
+            ret[1] = μ[1] * (-6.0)  # (1,1)
+            ret[2] = μ[1] * ( 2.0)  # (2,1)
+            ret[3] = μ[1] * (-2.0)  # (2,2)
+        end,
+    )
+
+    # Same constraint, but expressed as a VectorNonlinearOracle on [x[1], x[2]]
+    @constraint(model, c, [x[1], x[2]] in set)
+
+    return model
+end
diff --git a/test/nlp_problems/hs100.jl b/test/nlp_problems/hs100.jl
index 92156cb..75e640b 100644
--- 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
index e349d98..68f78d2 100644
--- a/test/nlp_problems/hs14.jl
+++ b/test/nlp_problems/hs14.jl
@@ -14,3 +14,56 @@ function hs14()
 
   return nlp
 end
+
+"Problem 14 in the Hock–Schittkowski suite, with inequality as a VectorNonlinearOracle"
+function hs14_oracle()
+    model = Model()
+
+    @variable(model, x[1:2])
+    set_start_value(x[1], 2)
+    set_start_value(x[2], 2)
+
+    @objective(model, Min, (x[1] - 2)^2 + (x[2] - 1)^2)
+
+    # Inequality: -x1^2/4 - x2^2 + 1 ≥ 0
+    # g(x) = -x1^2/4 - x2^2 + 1,  0 ≤ g(x) ≤ +∞
+    set = MOI.VectorNonlinearOracle(;
+        dimension = 2,                # 2 inputs: x1, x2
+        l = [-1.0],                    # lower bound on g(x)
+        u = [Inf],                    # upper bound on g(x)
+        eval_f = (ret, xv) -> begin
+            # ret[1] = g(x)
+            ret[1] = -0.25 * xv[1]^2 - xv[2]^2
+        end,
+        # Jacobian of g(x): ∇g = [∂g/∂x1, ∂g/∂x2]
+        # ∂g/∂x1 = -0.5 x1
+        # ∂g/∂x2 = -2 x2
+        jacobian_structure = [(1, 1), (1, 2)],
+        eval_jacobian = (ret, xv) -> begin
+            ret[1] = -0.5 * xv[1]     # d g / d x1
+            ret[2] = -2.0 * xv[2]     # d g / d x2
+        end,
+        # Hessian of g(x):
+        # ∂²g/∂x1² = -0.5
+        # ∂²g/∂x1∂x2 = 0
+        # ∂²g/∂x2² = -2
+        #
+        # We store only lower-triangular entries: (1,1), (2,1), (2,2)
+        hessian_lagrangian_structure = [(1, 1), (2, 1), (2, 2)],
+        eval_hessian_lagrangian = (ret, xv, μ) -> begin
+            # Hessian of μ[1] * g(x)
+            ret[1] = μ[1] * (-0.5)   # (1,1)
+            ret[2] = μ[1] * 0.0      # (2,1)
+            ret[3] = μ[1] * (-2.0)   # (2,2)
+        end,
+    )
+
+    # Inequality as oracle
+    @constraint(model, [x[1], x[2]] in set)
+
+    # Equality stays as a standard nonlinear constraint:
+    # x1 - 2 x2 + 1 == 0
+    @constraint(model, x[1] - 2 * x[2] + 1 == 0)
+
+    return model
+end
diff --git a/test/nlp_problems/hs61.jl b/test/nlp_problems/hs61.jl
index 5e422cc..384b223 100644
--- 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/nls_consistency.jl b/test/nls_consistency.jl
index 38a2159..353d856 100644
--- a/test/nls_consistency.jl
+++ b/test/nls_consistency.jl
@@ -8,7 +8,9 @@ for problem in nls_problems
     if isdefined(Main, Symbol(spc))
       push!(nlss, eval(Meta.parse(spc))())
     end
-    consistent_nlss(nlss, linear_api = true, test_slack = false, test_counters = false)
+    @testset "linear_api = $linear_api" for linear_api in (false, true)
+      consistent_nlss(nlss, linear_api = linear_api, test_slack = false, test_counters = false)
+    end
     view_subarray_nls(nls_moi)
   end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7125c6d..7588866 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,6 +8,9 @@ nls_problems = NLPModelsTest.nls_problems
 extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon", "operatorspb", "nf"]
 extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"]
 
+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"))
 end
@@ -19,6 +22,8 @@ end
 include("test_moi_nlp_model.jl")
 include("test_moi_nls_model.jl")
 
+include("test_moi_nlp_oracle.jl")
+
 include("nlp_consistency.jl")
 include("nls_consistency.jl")
 
diff --git a/test/test_moi_nlp_model.jl b/test/test_moi_nlp_model.jl
index bfafe05..536425a 100644
--- a/test/test_moi_nlp_model.jl
+++ b/test/test_moi_nlp_model.jl
@@ -11,7 +11,7 @@ println("Testing MathOptNLPModel")
   "‖c(x₀)‖"
 )
 # Test that every problem can be instantiated.
-for prob in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_problems))
+for prob in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_problems ∪ extra_nlp_oracle_problems))
   prob_fn = eval(prob)
   nlp = MathOptNLPModel(prob_fn(), hessian = (prob != :nohesspb))
   n = nlp.meta.nvar
diff --git a/test/test_moi_nlp_oracle.jl b/test/test_moi_nlp_oracle.jl
new file mode 100644
index 0000000..49853a2
--- /dev/null
+++ b/test/test_moi_nlp_oracle.jl
@@ -0,0 +1,43 @@
+@testset "Testing Oracle vs Non-Oracle NLP Models: $prob_oracle" for prob_oracle in extra_nlp_oracle_problems
+    prob_no_oracle = replace(prob_oracle, "_oracle" => "")
+    prob_fn_no_oracle = eval(prob_no_oracle |> Symbol)
+    prob_fn_oracle = eval(prob_oracle |> Symbol)
+
+    nlp_no_oracle = MathOptNLPModel(prob_fn_no_oracle(), hessian = true)
+    nlp_oracle = MathOptNLPModel(prob_fn_oracle(), hessian = true)
+
+    n = nlp_no_oracle.meta.nvar
+    m = nlp_no_oracle.meta.ncon
+    x = nlp_no_oracle.meta.x0
+
+    # Objective value
+    fx_no_oracle   = obj(nlp_no_oracle, x)
+    fx_with_oracle = obj(nlp_oracle, x)
+    @test isapprox(fx_no_oracle, fx_with_oracle; atol = 1e-8, rtol = 1e-8)
+
+    # Gradient of objective
+    ngx_no_oracle   = grad(nlp_no_oracle, x)
+    ngx_with_oracle = grad(nlp_oracle, x)
+    @test isapprox(ngx_no_oracle, ngx_with_oracle; atol = 1e-8, rtol = 1e-8)
+
+    # Constraint values (up to ordering)
+    ncx_no_oracle   = cons(nlp_no_oracle, x)
+    ncx_with_oracle = cons(nlp_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_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_oracle, x, λ)
+
+    @test isapprox(Matrix(H_no), Matrix(H_with); atol = 1e-8, rtol = 1e-8)
+end