diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 15d21a8..dab5d3c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -45,6 +45,13 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v1
+      - name: MOI
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="MathOptInterface", rev="bl/linearity"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         env:
diff --git a/ext/IpoptMathOptInterfaceExt/IpoptMathOptInterfaceExt.jl b/ext/IpoptMathOptInterfaceExt/IpoptMathOptInterfaceExt.jl
index 3184476..1b420f7 100644
--- a/ext/IpoptMathOptInterfaceExt/IpoptMathOptInterfaceExt.jl
+++ b/ext/IpoptMathOptInterfaceExt/IpoptMathOptInterfaceExt.jl
@@ -18,62 +18,62 @@ end
 
 include("MOI_wrapper.jl")
 
-PrecompileTools.@setup_workload begin
-    PrecompileTools.@compile_workload begin
-        model = MOI.Utilities.CachingOptimizer(
-            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
-            MOI.instantiate(Optimizer; with_bridge_type = Float64),
-        )
-        # We don't want to advertise this option, but it's required so that
-        # we don't print the banner during precompilation.
-        MOI.set(model, MOI.RawOptimizerAttribute("sb"), "yes")
-        MOI.set(model, MOI.Silent(), true)
-        x = MOI.add_variables(model, 3)
-        MOI.supports(model, MOI.VariableName(), typeof(x[1]))
-        MOI.set(model, MOI.VariableName(), x[1], "x1")
-        MOI.set(model, MOI.VariablePrimalStart(), x[1], 0.0)
-        for F in (MOI.VariableIndex, MOI.ScalarAffineFunction{Float64})
-            MOI.supports_constraint(model, F, MOI.GreaterThan{Float64})
-            MOI.supports_constraint(model, F, MOI.LessThan{Float64})
-            MOI.supports_constraint(model, F, MOI.EqualTo{Float64})
-            # These return false, but it doesn't matter
-            MOI.supports_constraint(model, F, MOI.ZeroOne)
-            MOI.supports_constraint(model, F, MOI.Integer)
-        end
-        MOI.add_constraint(model, x[1], MOI.GreaterThan(0.0))
-        MOI.add_constraint(model, x[2], MOI.LessThan(0.0))
-        MOI.add_constraint(model, x[3], MOI.EqualTo(0.0))
-        f = 1.0 * x[1] + x[2] + x[3]
-        c1 = MOI.add_constraint(model, f, MOI.GreaterThan(0.0))
-        MOI.set(model, MOI.ConstraintName(), c1, "c1")
-        MOI.supports(model, MOI.ConstraintName(), typeof(c1))
-        MOI.add_constraint(model, f, MOI.LessThan(0.0))
-        MOI.add_constraint(model, f, MOI.EqualTo(0.0))
-        y, _ = MOI.add_constrained_variables(model, MOI.Nonnegatives(2))
-        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-        MOI.supports(model, MOI.ObjectiveFunction{typeof(f)}())
-        MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-        MOI.set(
-            model,
-            MOI.NLPBlock(),
-            MOI.NLPBlockData([], _EmptyNLPEvaluator(), false),
-        )
-        MOI.optimize!(model)
-        MOI.get(model, MOI.TerminationStatus())
-        MOI.get(model, MOI.PrimalStatus())
-        MOI.get(model, MOI.DualStatus())
-        MOI.get(model, MOI.VariablePrimal(), x)
-        # We put these after `optimize!` so that the error is thrown on add,
-        # not on optimize!
-        try
-            MOI.add_constraint(model, x[1], MOI.ZeroOne())
-        catch
-        end
-        try
-            MOI.add_constraint(model, x[1], MOI.Integer())
-        catch
-        end
-    end
-end
+#PrecompileTools.@setup_workload begin
+#    PrecompileTools.@compile_workload begin
+#        model = MOI.Utilities.CachingOptimizer(
+#            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+#            MOI.instantiate(Optimizer; with_bridge_type = Float64),
+#        )
+#        # We don't want to advertise this option, but it's required so that
+#        # we don't print the banner during precompilation.
+#        MOI.set(model, MOI.RawOptimizerAttribute("sb"), "yes")
+#        MOI.set(model, MOI.Silent(), true)
+#        x = MOI.add_variables(model, 3)
+#        MOI.supports(model, MOI.VariableName(), typeof(x[1]))
+#        MOI.set(model, MOI.VariableName(), x[1], "x1")
+#        MOI.set(model, MOI.VariablePrimalStart(), x[1], 0.0)
+#        for F in (MOI.VariableIndex, MOI.ScalarAffineFunction{Float64})
+#            MOI.supports_constraint(model, F, MOI.GreaterThan{Float64})
+#            MOI.supports_constraint(model, F, MOI.LessThan{Float64})
+#            MOI.supports_constraint(model, F, MOI.EqualTo{Float64})
+#            # These return false, but it doesn't matter
+#            MOI.supports_constraint(model, F, MOI.ZeroOne)
+#            MOI.supports_constraint(model, F, MOI.Integer)
+#        end
+#        MOI.add_constraint(model, x[1], MOI.GreaterThan(0.0))
+#        MOI.add_constraint(model, x[2], MOI.LessThan(0.0))
+#        MOI.add_constraint(model, x[3], MOI.EqualTo(0.0))
+#        f = 1.0 * x[1] + x[2] + x[3]
+#        c1 = MOI.add_constraint(model, f, MOI.GreaterThan(0.0))
+#        MOI.set(model, MOI.ConstraintName(), c1, "c1")
+#        MOI.supports(model, MOI.ConstraintName(), typeof(c1))
+#        MOI.add_constraint(model, f, MOI.LessThan(0.0))
+#        MOI.add_constraint(model, f, MOI.EqualTo(0.0))
+#        y, _ = MOI.add_constrained_variables(model, MOI.Nonnegatives(2))
+#        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+#        MOI.supports(model, MOI.ObjectiveFunction{typeof(f)}())
+#        MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+#        MOI.set(
+#            model,
+#            MOI.NLPBlock(),
+#            MOI.NLPBlockData([], _EmptyNLPEvaluator(), false),
+#        )
+#        MOI.optimize!(model)
+#        MOI.get(model, MOI.TerminationStatus())
+#        MOI.get(model, MOI.PrimalStatus())
+#        MOI.get(model, MOI.DualStatus())
+#        MOI.get(model, MOI.VariablePrimal(), x)
+#        # We put these after `optimize!` so that the error is thrown on add,
+#        # not on optimize!
+#        try
+#            MOI.add_constraint(model, x[1], MOI.ZeroOne())
+#        catch
+#        end
+#        try
+#            MOI.add_constraint(model, x[1], MOI.Integer())
+#        catch
+#        end
+#    end
+#end
 
 end  # module IpoptMathOptInterfaceExt
diff --git a/ext/IpoptMathOptInterfaceExt/MOI_wrapper.jl b/ext/IpoptMathOptInterfaceExt/MOI_wrapper.jl
index 74fdf4b..ae87663 100644
--- a/ext/IpoptMathOptInterfaceExt/MOI_wrapper.jl
+++ b/ext/IpoptMathOptInterfaceExt/MOI_wrapper.jl
@@ -3,7 +3,16 @@
 # Use of this source code is governed by an MIT-style license that can be found
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
-include("utils.jl")
+const QPBlockData = MOI.Nonlinear.QPBlockData
+
+"""
+    _VectorNonlinearOracle
+
+An alias for [`MOI.VectorNonlinearOracle`](@ref), kept for backwards
+compatibility. The set used to be defined here before it was moved to
+MathOptInterface.
+"""
+const _VectorNonlinearOracle = MOI.VectorNonlinearOracle
 
 const _PARAMETER_OFFSET = 0x00f0000000000000
 
@@ -15,155 +24,6 @@ function _is_parameter(term::MOI.ScalarQuadraticTerm)
     return _is_parameter(term.variable_1) || _is_parameter(term.variable_2)
 end
 
-"""
-    _VectorNonlinearOracle(;
-        dimension::Int,
-        l::Vector{Float64},
-        u::Vector{Float64},
-        eval_f::Function,
-        jacobian_structure::Vector{Tuple{Int,Int}},
-        eval_jacobian::Function,
-        hessian_lagrangian_structure::Vector{Tuple{Int,Int}} = Tuple{Int,Int}[],
-        eval_hessian_lagrangian::Union{Nothing,Function} = nothing,
-    ) <: MOI.AbstractVectorSet
-
-The set:
-```math
-S = \\{x \\in \\mathbb{R}^{dimension}: l \\le f(x) \\le u\\}
-```
-where ``f`` is defined by the vectors `l` and `u`, and the callback oracles
-`eval_f`, `eval_jacobian`, and `eval_hessian_lagrangian`.
-
-!!! warning
-    This set is experimental. We will decide by September 30, 2025, whether to
-    convert this into the public `Ipopt.VectorNonlinearOracle`, move it to
-    `MOI.VectorNonlinearOracle`, or remove it completely.
-
-## f
-
-The `eval_f` function must have the signature
-```julia
-eval_f(ret::AbstractVector, x::AbstractVector)::Nothing
-```
-which fills ``f(x)`` into the dense vector `ret`.
-
-## Jacobian
-
-The `eval_jacobian` function must have the signature
-```julia
-eval_jacobian(ret::AbstractVector, x::AbstractVector)::Nothing
-```
-which fills the sparse Jacobian ``\\nabla f(x)`` into `ret`.
-
-The one-indexed sparsity structure must be provided in the `jacobian_structure`
-argument.
-
-## Hessian
-
-The `eval_hessian_lagrangian` function is optional.
-
-If `eval_hessian_lagrangian === nothing`, Ipopt will use a Hessian approximation
-instead of the exact Hessian.
-
-If `eval_hessian_lagrangian` is a function, it must have the signature
-```julia
-eval_hessian_lagrangian(
-    ret::AbstractVector,
-    x::AbstractVector,
-    μ::AbstractVector,
-)::Nothing
-```
-which fills the sparse Hessian of the Lagrangian ``\\sum \\mu_i \\nabla^2 f_i(x)``
-into `ret`.
-
-The one-indexed sparsity structure must be provided in the
-`hessian_lagrangian_structure` argument.
-
-## Example
-
-To model the set:
-```math
-\\begin{align}
-0 \\le & x^2           \\le 1
-0 \\le & y^2 + z^3 - w \\le 0
-\\end{align}
-```
-do
-```jldoctest
-julia> import Ipopt
-
-julia> set = Ipopt._VectorNonlinearOracle(;
-           dimension = 3,
-           l = [0.0, 0.0],
-           u = [1.0, 0.0],
-           eval_f = (ret, x) -> begin
-               ret[1] = x[2]^2
-               ret[2] = x[3]^2 + x[4]^3 - x[1]
-               return
-           end,
-           jacobian_structure = [(1, 2), (2, 1), (2, 3), (2, 4)],
-           eval_jacobian = (ret, x) -> begin
-               ret[1] = 2.0 * x[2]
-               ret[2] = -1.0
-               ret[3] = 2.0 * x[3]
-               ret[4] = 3.0 * x[4]^2
-               return
-           end,
-           hessian_lagrangian_structure = [(2, 2), (3, 3), (4, 4)],
-           eval_hessian_lagrangian = (ret, x, u) -> begin
-               ret[1] = 2.0 * u[1]
-               ret[2] = 2.0 * u[2]
-               ret[3] = 6.0 * x[4] * u[2]
-               return
-           end,
-       );
-```
-"""
-struct _VectorNonlinearOracle <: MOI.AbstractVectorSet
-    input_dimension::Int
-    output_dimension::Int
-    l::Vector{Float64}
-    u::Vector{Float64}
-    eval_f::Function
-    jacobian_structure::Vector{Tuple{Int,Int}}
-    eval_jacobian::Function
-    hessian_lagrangian_structure::Vector{Tuple{Int,Int}}
-    eval_hessian_lagrangian::Union{Nothing,Function}
-    # Temporary storage
-    x::Vector{Float64}
-
-    function _VectorNonlinearOracle(;
-        dimension::Int,
-        l::Vector{Float64},
-        u::Vector{Float64},
-        eval_f::Function,
-        jacobian_structure::Vector{Tuple{Int,Int}},
-        eval_jacobian::Function,
-        # The hessian_lagrangian is optional.
-        hessian_lagrangian_structure::Vector{Tuple{Int,Int}} = Tuple{Int,Int}[],
-        eval_hessian_lagrangian::Union{Nothing,Function} = nothing,
-    )
-        @assert length(l) == length(u)
-        return new(
-            dimension,
-            length(l),
-            l,
-            u,
-            eval_f,
-            jacobian_structure,
-            eval_jacobian,
-            hessian_lagrangian_structure,
-            eval_hessian_lagrangian,
-            # Temporary storage
-            zeros(dimension),
-        )
-    end
-end
-
-MOI.dimension(s::_VectorNonlinearOracle) = s.input_dimension
-
-MOI.copy(s::_VectorNonlinearOracle) = s
-
 """
     Optimizer()
 
@@ -195,7 +55,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     ad_backend::MOI.Nonlinear.AbstractAutomaticDifferentiation
 
     vector_nonlinear_oracle_constraints::Vector{
-        Tuple{MOI.VectorOfVariables,_VectorNonlinearOracle},
+        Tuple{MOI.VectorOfVariables,MOI.VectorNonlinearOracle{Float64}},
     }
 
     function Optimizer()
@@ -221,7 +81,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
             nothing,
             0,
             MOI.Nonlinear.SparseReverseMode(),
-            Tuple{MOI.VectorOfVariables,_VectorNonlinearOracle}[],
+            Tuple{MOI.VectorOfVariables,MOI.VectorNonlinearOracle{Float64}}[],
         )
     end
 end
@@ -246,6 +106,8 @@ MOI.jacobian_structure(::_EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
 MOI.hessian_lagrangian_structure(::_EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
 MOI.eval_constraint_jacobian(::_EmptyNLPEvaluator, J, x) = nothing
 MOI.eval_hessian_lagrangian(::_EmptyNLPEvaluator, H, x, σ, μ) = nothing
+MOI.Nonlinear.num_constraints(::_EmptyNLPEvaluator) = 0
+MOI.Nonlinear.constraint_bounds(::_EmptyNLPEvaluator) = MOI.NLPBoundsPair[]
 
 function MOI.empty!(model::Optimizer)
     model.inner = nothing
@@ -312,6 +174,11 @@ function MOI.add_constrained_variable(
     push!(model.list_of_variable_indices, p)
     model.parameters[p] =
         MOI.Nonlinear.add_parameter(model.nlp_model, set.value)
+    # Register the parameter in the QP block. `QPBlockData` treats a variable
+    # as a parameter if and only if its index is a key of `parameters`, so
+    # this must happen before any structure query. The value is re-synced in
+    # `copy_parameters` before every solve.
+    model.qp_data.parameters[p.value] = set.value
     ci = MOI.ConstraintIndex{MOI.VariableIndex,typeof(set)}(p.value)
     return p, ci
 end
@@ -402,7 +269,7 @@ function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintTypesPresent)
     append!(ret, MOI.get(model.qp_data, attr))
     _add_scalar_nonlinear_constraints(ret, model.nlp_model)
     if !isempty(model.vector_nonlinear_oracle_constraints)
-        push!(ret, (MOI.VectorOfVariables, _VectorNonlinearOracle))
+        push!(ret, (MOI.VectorOfVariables, MOI.VectorNonlinearOracle{Float64}))
     end
     return ret
 end
@@ -796,19 +663,19 @@ function MOI.set(
     return
 end
 
-### MOI.VectorOfVariables in _VectorNonlinearOracle
+### MOI.VectorOfVariables in MOI.VectorNonlinearOracle
 
 function MOI.supports_constraint(
     ::Optimizer,
     ::Type{MOI.VectorOfVariables},
-    ::Type{_VectorNonlinearOracle},
+    ::Type{<:MOI.VectorNonlinearOracle},
 )
     return true
 end
 
 function MOI.is_valid(
     model::Optimizer,
-    ci::MOI.ConstraintIndex{MOI.VectorOfVariables,_VectorNonlinearOracle},
+    ci::MOI.ConstraintIndex{MOI.VectorOfVariables,<:MOI.VectorNonlinearOracle},
 )
     return 1 <= ci.value <= length(model.vector_nonlinear_oracle_constraints)
 end
@@ -816,7 +683,7 @@ end
 function MOI.get(
     model::Optimizer,
     attr::MOI.ListOfConstraintIndices{F,S},
-) where {F<:MOI.VectorOfVariables,S<:_VectorNonlinearOracle}
+) where {F<:MOI.VectorOfVariables,S<:MOI.VectorNonlinearOracle}
     n = length(model.vector_nonlinear_oracle_constraints)
     return MOI.ConstraintIndex{F,S}.(1:n)
 end
@@ -824,7 +691,7 @@ end
 function MOI.get(
     model::Optimizer,
     attr::MOI.NumberOfConstraints{F,S},
-) where {F<:MOI.VectorOfVariables,S<:_VectorNonlinearOracle}
+) where {F<:MOI.VectorOfVariables,S<:MOI.VectorNonlinearOracle}
     return length(model.vector_nonlinear_oracle_constraints)
 end
 
@@ -832,7 +699,7 @@ function MOI.add_constraint(
     model::Optimizer,
     f::F,
     s::S,
-) where {F<:MOI.VectorOfVariables,S<:_VectorNonlinearOracle}
+) where {F<:MOI.VectorOfVariables,S<:MOI.VectorNonlinearOracle{Float64}}
     model.inner = nothing
     push!(model.vector_nonlinear_oracle_constraints, (f, s))
     n = length(model.vector_nonlinear_oracle_constraints)
@@ -842,7 +709,7 @@ end
 function row(
     model::Optimizer,
     ci::MOI.ConstraintIndex{F,S},
-) where {F<:MOI.VectorOfVariables,S<:_VectorNonlinearOracle}
+) where {F<:MOI.VectorOfVariables,S<:MOI.VectorNonlinearOracle}
     offset = length(model.qp_data)
     for i in 1:(ci.value-1)
         _, s = model.vector_nonlinear_oracle_constraints[i]
@@ -856,7 +723,7 @@ function MOI.get(
     model::Optimizer,
     attr::MOI.ConstraintPrimal,
     ci::MOI.ConstraintIndex{F,S},
-) where {F<:MOI.VectorOfVariables,S<:_VectorNonlinearOracle}
+) where {F<:MOI.VectorOfVariables,S<:MOI.VectorNonlinearOracle}
     MOI.check_result_index_bounds(model, attr)
     MOI.throw_if_not_valid(model, ci)
     f, _ = model.vector_nonlinear_oracle_constraints[ci.value]
@@ -867,25 +734,35 @@ function MOI.get(
     model::Optimizer,
     attr::MOI.ConstraintDual,
     ci::MOI.ConstraintIndex{F,S},
-) where {F<:MOI.VectorOfVariables,S<:_VectorNonlinearOracle}
+) where {F<:MOI.VectorOfVariables,S<:MOI.VectorNonlinearOracle}
     MOI.check_result_index_bounds(model, attr)
     MOI.throw_if_not_valid(model, ci)
     sign = -_dual_multiplier(model)
     f, s = model.vector_nonlinear_oracle_constraints[ci.value]
     λ = model.inner.mult_g[row(model, ci)]
-    J = Tuple{Int,Int}[]
-    _jacobian_structure(J, 0, f, s)
-    J_val = zeros(length(J))
-    _eval_constraint_jacobian(J_val, 0, model.inner.x, f, s)
+    x = [model.inner.x[v.value] for v in f.variables]
+    J_val = zeros(length(s.jacobian_structure))
+    s.eval_jacobian(J_val, x)
     dual = zeros(MOI.dimension(s))
-    # dual = λ' * J(x)
-    col_to_index = Dict(x.value => j for (j, x) in enumerate(f.variables))
-    for ((row, col), J_rc) in zip(J, J_val)
-        dual[col_to_index[col]] += sign * J_rc * λ[row]
+    # dual = λ' * J(x). The columns of `s.jacobian_structure` are indices
+    # into `f.variables`.
+    for ((r, c), J_rc) in zip(s.jacobian_structure, J_val)
+        dual[c] += sign * J_rc * λ[r]
     end
     return dual
 end
 
+function MOI.get(
+    model::Optimizer,
+    attr::MOI.LagrangeMultiplier,
+    ci::MOI.ConstraintIndex{F,S},
+) where {F<:MOI.VectorOfVariables,S<:MOI.VectorNonlinearOracle}
+    MOI.check_result_index_bounds(model, attr)
+    MOI.throw_if_not_valid(model, ci)
+    s = -_dual_multiplier(model)
+    return s .* model.inner.mult_g[row(model, ci)]
+end
+
 ### UserDefinedFunction
 
 MOI.supports(model::Optimizer, ::MOI.UserDefinedFunction) = true
@@ -1138,172 +1015,6 @@ function MOI.set(
     return
 end
 
-### Eval_F_CB
-
-function MOI.eval_objective(model::Optimizer, x)
-    # TODO(odow): FEASIBILITY_SENSE could produce confusing solver output if
-    # a nonzero objective is set.
-    if model.sense == MOI.FEASIBILITY_SENSE
-        return 0.0
-    elseif model.nlp_data.has_objective
-        return MOI.eval_objective(model.nlp_data.evaluator, x)::Float64
-    end
-    return MOI.eval_objective(model.qp_data, x)
-end
-
-### Eval_Grad_F_CB
-
-function MOI.eval_objective_gradient(model::Optimizer, grad, x)
-    if model.sense == MOI.FEASIBILITY_SENSE
-        grad .= zero(eltype(grad))
-    elseif model.nlp_data.has_objective
-        MOI.eval_objective_gradient(model.nlp_data.evaluator, grad, x)
-    else
-        MOI.eval_objective_gradient(model.qp_data, grad, x)
-    end
-    return
-end
-
-### Eval_G_CB
-
-function _eval_constraint(
-    g::AbstractVector,
-    offset::Int,
-    x::AbstractVector,
-    f::MOI.VectorOfVariables,
-    s::_VectorNonlinearOracle,
-)
-    for i in 1:s.input_dimension
-        s.x[i] = x[f.variables[i].value]
-    end
-    ret = view(g, offset .+ (1:s.output_dimension))
-    s.eval_f(ret, s.x)
-    return offset + s.output_dimension
-end
-
-function MOI.eval_constraint(model::Optimizer, g, x)
-    MOI.eval_constraint(model.qp_data, g, x)
-    offset = length(model.qp_data)
-    for (f, s) in model.vector_nonlinear_oracle_constraints
-        offset = _eval_constraint(g, offset, x, f, s)
-    end
-    g_nlp = view(g, (offset+1):length(g))
-    MOI.eval_constraint(model.nlp_data.evaluator, g_nlp, x)
-    return
-end
-
-### Eval_Jac_G_CB
-
-function _jacobian_structure(
-    ret::AbstractVector,
-    row_offset::Int,
-    f::MOI.VectorOfVariables,
-    s::_VectorNonlinearOracle,
-)
-    for (i, j) in s.jacobian_structure
-        push!(ret, (row_offset + i, f.variables[j].value))
-    end
-    return row_offset + s.output_dimension
-end
-
-function MOI.jacobian_structure(model::Optimizer)
-    J = MOI.jacobian_structure(model.qp_data)
-    offset = length(model.qp_data)
-    for (f, s) in model.vector_nonlinear_oracle_constraints
-        offset = _jacobian_structure(J, offset, f, s)
-    end
-    if length(model.nlp_data.constraint_bounds) > 0
-        J_nlp = MOI.jacobian_structure(
-            model.nlp_data.evaluator,
-        )::Vector{Tuple{Int64,Int64}}
-        for (row, col) in J_nlp
-            push!(J, (row + offset, col))
-        end
-    end
-    return J
-end
-
-function _eval_constraint_jacobian(
-    values::AbstractVector,
-    offset::Int,
-    x::AbstractVector,
-    f::MOI.VectorOfVariables,
-    s::_VectorNonlinearOracle,
-)
-    for i in 1:s.input_dimension
-        s.x[i] = x[f.variables[i].value]
-    end
-    nnz = length(s.jacobian_structure)
-    s.eval_jacobian(view(values, offset .+ (1:nnz)), s.x)
-    return offset + nnz
-end
-
-function MOI.eval_constraint_jacobian(model::Optimizer, values, x)
-    offset = MOI.eval_constraint_jacobian(model.qp_data, values, x)
-    offset -= 1  # .qp_data returns one-indexed offset
-    for (f, s) in model.vector_nonlinear_oracle_constraints
-        offset = _eval_constraint_jacobian(values, offset, x, f, s)
-    end
-    nlp_values = view(values, (offset+1):length(values))
-    MOI.eval_constraint_jacobian(model.nlp_data.evaluator, nlp_values, x)
-    return
-end
-
-### Eval_H_CB
-
-function _hessian_lagrangian_structure(
-    ret::AbstractVector,
-    f::MOI.VectorOfVariables,
-    s::_VectorNonlinearOracle,
-)
-    for (i, j) in s.hessian_lagrangian_structure
-        push!(ret, (f.variables[i].value, f.variables[j].value))
-    end
-    return
-end
-
-function MOI.hessian_lagrangian_structure(model::Optimizer)
-    H = MOI.hessian_lagrangian_structure(model.qp_data)
-    for (f, s) in model.vector_nonlinear_oracle_constraints
-        _hessian_lagrangian_structure(H, f, s)
-    end
-    append!(H, MOI.hessian_lagrangian_structure(model.nlp_data.evaluator))
-    return H
-end
-
-function _eval_hessian_lagrangian(
-    H::AbstractVector,
-    H_offset::Int,
-    x::AbstractVector,
-    μ::AbstractVector,
-    μ_offset::Int,
-    f::MOI.VectorOfVariables,
-    s::_VectorNonlinearOracle,
-)
-    for i in 1:s.input_dimension
-        s.x[i] = x[f.variables[i].value]
-    end
-    H_nnz = length(s.hessian_lagrangian_structure)
-    H_view = view(H, H_offset .+ (1:H_nnz))
-    μ_view = view(μ, μ_offset .+ (1:s.output_dimension))
-    s.eval_hessian_lagrangian(H_view, s.x, μ_view)
-    return H_offset + H_nnz, μ_offset + s.output_dimension
-end
-
-function MOI.eval_hessian_lagrangian(model::Optimizer, H, x, σ, μ)
-    offset = MOI.eval_hessian_lagrangian(model.qp_data, H, x, σ, μ)
-    offset -= 1  # .qp_data returns one-indexed offset
-    μ_offset = length(model.qp_data)
-    for (f, s) in model.vector_nonlinear_oracle_constraints
-        offset, μ_offset =
-            _eval_hessian_lagrangian(H, offset, x, μ, μ_offset, f, s)
-    end
-    H_nlp = view(H, (offset+1):length(H))
-    μ_nlp = view(μ, (μ_offset+1):length(μ))
-    MOI.eval_hessian_lagrangian(model.nlp_data.evaluator, H_nlp, x, σ, μ_nlp)
-    return
-end
-
 ### MOI.AutomaticDifferentiationBackend
 
 MOI.supports(::Optimizer, ::MOI.AutomaticDifferentiationBackend) = true
@@ -1327,54 +1038,27 @@ end
 
 ### MOI.optimize!
 
-function _setup_model(model::Optimizer)
-    vars = MOI.get(model.variables, MOI.ListOfVariableIndices())
-    if isempty(vars)
-        # Don't attempt to create a problem because Ipopt will error.
-        model.invalid_model = true
-        return
-    end
-    if model.nlp_model !== nothing
-        model.nlp_data = MOI.NLPBlockData(
-            MOI.Nonlinear.Evaluator(model.nlp_model, model.ad_backend, vars),
-        )
-    end
-    has_quadratic_constraints =
-        any(isequal(_kFunctionTypeScalarQuadratic), model.qp_data.function_type)
-    has_nlp_constraints =
-        !isempty(model.nlp_data.constraint_bounds) ||
-        !isempty(model.vector_nonlinear_oracle_constraints)
-    has_hessian = :Hess in MOI.features_available(model.nlp_data.evaluator)
-    for (_, s) in model.vector_nonlinear_oracle_constraints
-        if s.eval_hessian_lagrangian === nothing
-            has_hessian = false
-            break
-        end
-    end
-    init_feat = [:Grad]
-    if has_hessian
-        push!(init_feat, :Hess)
-    end
-    if has_nlp_constraints
-        push!(init_feat, :Jac)
-    end
-    MOI.initialize(model.nlp_data.evaluator, init_feat)
-    jacobian_sparsity = MOI.jacobian_structure(model)
+# Function barrier as the type of `evaluator` depends
+# on the sparsity structure of the hessian through the
+# size of its seed matrix.
+function _setup_callbacks(model::Optimizer, evaluator, vars, has_hessian)
+    let evaluator = evaluator
+    jacobian_sparsity = MOI.jacobian_structure(evaluator)
     hessian_sparsity = if has_hessian
-        MOI.hessian_lagrangian_structure(model)
+        MOI.hessian_lagrangian_structure(evaluator)
     else
         Tuple{Int,Int}[]
     end
-    eval_f_cb(x) = MOI.eval_objective(model, x)
-    eval_grad_f_cb(x, grad_f) = MOI.eval_objective_gradient(model, grad_f, x)
-    eval_g_cb(x, g) = MOI.eval_constraint(model, g, x)
+    eval_f_cb(x) = MOI.eval_objective(evaluator, x)
+    eval_grad_f_cb(x, grad_f) = MOI.eval_objective_gradient(evaluator, grad_f, x)
+    eval_g_cb(x, g) = MOI.eval_constraint(evaluator, g, x)
     function eval_jac_g_cb(x, rows, cols, values)
         if values === nothing
             for i in 1:length(jacobian_sparsity)
                 rows[i], cols[i] = jacobian_sparsity[i]
             end
         else
-            MOI.eval_constraint_jacobian(model, values, x)
+            MOI.eval_constraint_jacobian(evaluator, values, x)
         end
         return
     end
@@ -1384,19 +1068,11 @@ function _setup_model(model::Optimizer)
                 rows[i], cols[i] = hessian_sparsity[i]
             end
         else
-            MOI.eval_hessian_lagrangian(model, values, x, obj_factor, lambda)
+            MOI.eval_hessian_lagrangian(evaluator, values, x, obj_factor, lambda)
         end
         return
     end
-    g_L, g_U = copy(model.qp_data.g_L), copy(model.qp_data.g_U)
-    for (_, s) in model.vector_nonlinear_oracle_constraints
-        append!(g_L, s.l)
-        append!(g_U, s.u)
-    end
-    for bound in model.nlp_data.constraint_bounds
-        push!(g_L, bound.lower)
-        push!(g_U, bound.upper)
-    end
+    g_L, g_U = _constraint_bounds(model)
     model.inner = Ipopt.CreateIpoptProblem(
         length(vars),
         model.variables.lower,
@@ -1412,6 +1088,91 @@ function _setup_model(model::Optimizer)
         eval_jac_g_cb,
         has_hessian ? eval_h_cb : nothing,
     )
+end
+end
+
+# The number of rows before the rows of `model.nlp_data`, that is, the rows
+# of the QP block and of the vector-nonlinear-oracle constraints.
+function _nlp_block_offset(model::Optimizer)
+    offset = length(model.qp_data)
+    for (_, s) in model.vector_nonlinear_oracle_constraints
+        offset += s.output_dimension
+    end
+    return offset
+end
+
+function _constraint_bounds(model::Optimizer)
+    g_L, g_U = copy(model.qp_data.g_L), copy(model.qp_data.g_U)
+    for (_, s) in model.vector_nonlinear_oracle_constraints
+        append!(g_L, s.l)
+        append!(g_U, s.u)
+    end
+    for bound in model.nlp_data.constraint_bounds
+        push!(g_L, bound.lower)
+        push!(g_U, bound.upper)
+    end
+    return g_L, g_U
+end
+
+function _setup_model(model::Optimizer)
+    vars = MOI.get(model.variables, MOI.ListOfVariableIndices())
+    if isempty(vars)
+        # Don't attempt to create a problem because Ipopt will error.
+        model.invalid_model = true
+        return
+    end
+    inner_model = MOI.Nonlinear.Model()
+    if model.nlp_model !== nothing
+        inner_model = model.nlp_model
+        model.nlp_data = MOI.NLPBlockData(
+            MOI.Nonlinear.Evaluator(model.nlp_model, model.ad_backend, vars),
+        )
+    end
+    # Assemble the MOI.Nonlinear layer stack over the existing storage. The
+    # row order is [qp, oracle, nlp], matching the layers' own-rows-first
+    # convention with the quad layer outermost.
+    oracle_model = MOI.Nonlinear.ModelWithOracles{Float64}(
+        model.vector_nonlinear_oracle_constraints,
+        inner_model,
+    )
+    objective_sink = if model.sense == MOI.FEASIBILITY_SENSE
+        # TODO(odow): FEASIBILITY_SENSE could produce confusing solver output
+        # if a nonzero objective is set.
+        :none
+    elseif model.nlp_data.has_objective
+        :inner
+    else
+        :quad
+    end
+    quad_model = MOI.Nonlinear.ModelWithQuad{Float64}(
+        model.qp_data,
+        oracle_model;
+        objective_sink = objective_sink,
+    )
+    evaluator = MOI.Nonlinear.EvaluatorWithQuad(
+        quad_model,
+        MOI.Nonlinear.EvaluatorWithOracles(
+            oracle_model,
+            model.nlp_data.evaluator,
+            vars,
+        ),
+        vars,
+    )
+    has_nlp_constraints =
+        !isempty(model.nlp_data.constraint_bounds) ||
+        !isempty(model.vector_nonlinear_oracle_constraints)
+    # The oracle layer removes :Hess from `MOI.features_available` if an
+    # oracle does not implement `eval_hessian_lagrangian`.
+    has_hessian = :Hess in MOI.features_available(evaluator)
+    init_feat = [:Grad]
+    if has_hessian
+        push!(init_feat, :Hess)
+    end
+    if has_nlp_constraints
+        push!(init_feat, :Jac)
+    end
+    MOI.initialize(evaluator, init_feat)
+    _setup_callbacks(model, evaluator, vars, has_hessian)
     if model.sense == MOI.MIN_SENSE
         Ipopt.AddIpoptNumOption(model.inner, "obj_scaling_factor", 1.0)
     elseif model.sense == MOI.MAX_SENSE
@@ -1430,14 +1191,14 @@ function _setup_model(model::Optimizer)
             "limited-memory",
         )
     end
-    if !has_nlp_constraints && !has_quadratic_constraints
+    linearity = MOI.Nonlinear.constraint_linearity(evaluator)
+    if linearity !== nothing && all(==(MOI.Nonlinear.LINEAR), linearity)
         Ipopt.AddIpoptStrOption(model.inner, "jac_c_constant", "yes")
         Ipopt.AddIpoptStrOption(model.inner, "jac_d_constant", "yes")
-        if !model.nlp_data.has_objective
-            # We turn on this option if all constraints are linear and the
-            # objective is linear or quadratic. From the documentation, it's
-            # unclear if it may also apply if the constraints are at most
-            # quadratic.
+        if MOI.Nonlinear.objective_linearity(evaluator) <=
+           MOI.Nonlinear.QUADRATIC
+            # All constraints are linear and the objective has a constant
+            # Hessian.
             Ipopt.AddIpoptStrOption(model.inner, "hessian_constant", "yes")
         end
     end
@@ -1495,9 +1256,9 @@ function MOI.optimize!(model::Optimizer)
     for (i, start) in enumerate(model.qp_data.mult_g)
         inner.mult_g[i] = _dual_start(model, start, -1)
     end
-    offset = length(model.qp_data.mult_g)
+    inner.mult_g[(length(model.qp_data.mult_g)+1):end] .= 0.0
+    offset = _nlp_block_offset(model)
     if model.nlp_dual_start === nothing
-        inner.mult_g[(offset+1):end] .= 0.0
         for (key, val) in model.mult_g_nlp
             inner.mult_g[offset+key.value] = val
         end
@@ -1588,12 +1349,8 @@ function _manually_evaluated_primal_status(model::Optimizer)
     x, g = model.inner.x, model.inner.g
     m, n = length(g), length(x)
     x_L, x_U = model.variables.lower, model.variables.upper
-    g_L, g_U = copy(model.qp_data.g_L), copy(model.qp_data.g_U)
-    # Assuming constraints are guaranteed to be in the order [qp_cons, nlp_cons]
-    for bound in model.nlp_data.constraint_bounds
-        push!(g_L, bound.lower)
-        push!(g_U, bound.upper)
-    end
+    # Constraints are in the order [qp_cons, oracle_cons, nlp_cons]
+    g_L, g_U = _constraint_bounds(model)
     # 1e-8 is the default tolerance
     tol = get(model.options, "tol", 1e-8)
     if all(x_L[i] - tol <= x[i] <= x_U[i] + tol for i in 1:n) &&
@@ -1784,7 +1541,7 @@ end
 function MOI.get(model::Optimizer, attr::MOI.NLPBlockDual)
     MOI.check_result_index_bounds(model, attr)
     s = -_dual_multiplier(model)
-    return s .* model.inner.mult_g[(length(model.qp_data)+1):end]
+    return s .* model.inner.mult_g[(_nlp_block_offset(model)+1):end]
 end
 
 ### Ipopt.CallbackFunction
diff --git a/ext/IpoptMathOptInterfaceExt/utils.jl b/ext/IpoptMathOptInterfaceExt/utils.jl
deleted file mode 100644
index cedaed5..0000000
--- a/ext/IpoptMathOptInterfaceExt/utils.jl
+++ /dev/null
@@ -1,541 +0,0 @@
-# Copyright (c) 2013: Iain Dunning, Miles Lubin, and contributors
-#
-# Use of this source code is governed by an MIT-style license that can be found
-# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
-
-# !!! warning
-#
-#     The contents of this file are experimental.
-#
-#     Until this message is removed, breaking changes to the functions and
-#     types, including their deletion, may be introduced in any minor or patch
-#     release of Ipopt.
-
-@enum(
-    _FunctionType,
-    _kFunctionTypeVariableIndex,
-    _kFunctionTypeScalarAffine,
-    _kFunctionTypeScalarQuadratic,
-)
-
-function _function_type_to_func(::Type{T}, k::_FunctionType) where {T}
-    if k == _kFunctionTypeVariableIndex
-        return MOI.VariableIndex
-    elseif k == _kFunctionTypeScalarAffine
-        return MOI.ScalarAffineFunction{T}
-    else
-        @assert k == _kFunctionTypeScalarQuadratic
-        return MOI.ScalarQuadraticFunction{T}
-    end
-end
-
-_function_info(::MOI.VariableIndex) = _kFunctionTypeVariableIndex
-_function_info(::MOI.ScalarAffineFunction) = _kFunctionTypeScalarAffine
-_function_info(::MOI.ScalarQuadraticFunction) = _kFunctionTypeScalarQuadratic
-
-@enum(
-    _BoundType,
-    _kBoundTypeLessThan,
-    _kBoundTypeGreaterThan,
-    _kBoundTypeEqualTo,
-    _kBoundTypeInterval,
-)
-
-_set_info(s::MOI.LessThan) = _kBoundTypeLessThan, -Inf, s.upper
-_set_info(s::MOI.GreaterThan) = _kBoundTypeGreaterThan, s.lower, Inf
-_set_info(s::MOI.EqualTo) = _kBoundTypeEqualTo, s.value, s.value
-_set_info(s::MOI.Interval) = _kBoundTypeInterval, s.lower, s.upper
-
-function _bound_type_to_set(::Type{T}, k::_BoundType) where {T}
-    if k == _kBoundTypeEqualTo
-        return MOI.EqualTo{T}
-    elseif k == _kBoundTypeLessThan
-        return MOI.LessThan{T}
-    elseif k == _kBoundTypeGreaterThan
-        return MOI.GreaterThan{T}
-    else
-        @assert k == _kBoundTypeInterval
-        return MOI.Interval{T}
-    end
-end
-
-mutable struct QPBlockData{T}
-    objective::Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}}
-    objective_function_type::_FunctionType
-    constraints::Vector{
-        Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}},
-    }
-    g_L::Vector{T}
-    g_U::Vector{T}
-    mult_g::Vector{Union{Nothing,T}}
-    function_type::Vector{_FunctionType}
-    bound_type::Vector{_BoundType}
-    parameters::Dict{Int64,T}
-
-    function QPBlockData{T}() where {T}
-        return new(
-            zero(MOI.ScalarQuadraticFunction{T}),
-            _kFunctionTypeScalarAffine,
-            Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}}[],
-            T[],
-            T[],
-            Union{Nothing,T}[],
-            _FunctionType[],
-            _BoundType[],
-            Dict{Int64,T}(),
-        )
-    end
-end
-
-function _value(variable::MOI.VariableIndex, x::Vector, p::Dict)
-    if _is_parameter(variable)
-        return p[variable.value]
-    else
-        return x[variable.value]
-    end
-end
-
-function eval_function(
-    f::MOI.ScalarQuadraticFunction{T},
-    x::Vector{T},
-    p::Dict{Int64,T},
-)::T where {T}
-    y = f.constant
-    for term in f.affine_terms
-        y += term.coefficient * _value(term.variable, x, p)
-    end
-    for term in f.quadratic_terms
-        v1 = _value(term.variable_1, x, p)
-        v2 = _value(term.variable_2, x, p)
-        if term.variable_1 == term.variable_2
-            y += term.coefficient * v1 * v2 / 2
-        else
-            y += term.coefficient * v1 * v2
-        end
-    end
-    return y
-end
-
-function eval_function(
-    f::MOI.ScalarAffineFunction{T},
-    x::Vector{T},
-    p::Dict{Int64,T},
-)::T where {T}
-    y = f.constant
-    for term in f.terms
-        y += term.coefficient * _value(term.variable, x, p)
-    end
-    return y
-end
-
-function eval_dense_gradient(
-    ∇f::Vector{T},
-    f::MOI.ScalarQuadraticFunction{T},
-    x::Vector{T},
-    p::Dict{Int64,T},
-)::Nothing where {T}
-    for term in f.affine_terms
-        if !_is_parameter(term.variable)
-            ∇f[term.variable.value] += term.coefficient
-        end
-    end
-    for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1)
-            v = _value(term.variable_2, x, p)
-            ∇f[term.variable_1.value] += term.coefficient * v
-        end
-        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
-            v = _value(term.variable_1, x, p)
-            ∇f[term.variable_2.value] += term.coefficient * v
-        end
-    end
-    return
-end
-
-function eval_dense_gradient(
-    ∇f::Vector{T},
-    f::MOI.ScalarAffineFunction{T},
-    x::Vector{T},
-    p::Dict{Int64,T},
-)::Nothing where {T}
-    for term in f.terms
-        if !_is_parameter(term.variable)
-            ∇f[term.variable.value] += term.coefficient
-        end
-    end
-    return
-end
-
-function append_sparse_gradient_structure!(
-    f::MOI.ScalarQuadraticFunction,
-    J,
-    row,
-)
-    for term in f.affine_terms
-        if !_is_parameter(term.variable)
-            push!(J, (row, term.variable.value))
-        end
-    end
-    for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1)
-            push!(J, (row, term.variable_1.value))
-        end
-        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
-            push!(J, (row, term.variable_2.value))
-        end
-    end
-    return
-end
-
-function append_sparse_gradient_structure!(f::MOI.ScalarAffineFunction, J, row)
-    for term in f.terms
-        if !_is_parameter(term.variable)
-            push!(J, (row, term.variable.value))
-        end
-    end
-    return
-end
-
-function eval_sparse_gradient(
-    ∇f::AbstractVector{T},
-    f::MOI.ScalarQuadraticFunction{T},
-    x::Vector{T},
-    p::Dict{Int64,T},
-)::Int where {T}
-    i = 0
-    for term in f.affine_terms
-        if !_is_parameter(term.variable)
-            i += 1
-            ∇f[i] = term.coefficient
-        end
-    end
-    for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1)
-            v = _value(term.variable_2, x, p)
-            i += 1
-            ∇f[i] = term.coefficient * v
-        end
-        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
-            v = _value(term.variable_1, x, p)
-            i += 1
-            ∇f[i] = term.coefficient * v
-        end
-    end
-    return i
-end
-
-function eval_sparse_gradient(
-    ∇f::AbstractVector{T},
-    f::MOI.ScalarAffineFunction{T},
-    x::Vector{T},
-    p::Dict{Int64,T},
-)::Int where {T}
-    i = 0
-    for term in f.terms
-        if !_is_parameter(term.variable)
-            i += 1
-            ∇f[i] = term.coefficient
-        end
-    end
-    return i
-end
-
-function append_sparse_hessian_structure!(f::MOI.ScalarQuadraticFunction, H)
-    for term in f.quadratic_terms
-        if _is_parameter(term.variable_1) || _is_parameter(term.variable_2)
-            continue
-        end
-        push!(H, (term.variable_1.value, term.variable_2.value))
-    end
-    return
-end
-
-append_sparse_hessian_structure!(::MOI.ScalarAffineFunction, H) = nothing
-
-function eval_sparse_hessian(
-    ∇²f::AbstractVector{T},
-    f::MOI.ScalarQuadraticFunction{T},
-    σ::T,
-)::Int where {T}
-    i = 0
-    for term in f.quadratic_terms
-        if _is_parameter(term.variable_1) || _is_parameter(term.variable_2)
-            continue
-        end
-        i += 1
-        ∇²f[i] = term.coefficient * σ
-    end
-    return i
-end
-
-function eval_sparse_hessian(
-    ∇²f::AbstractVector{T},
-    f::MOI.ScalarAffineFunction{T},
-    σ::T,
-)::Int where {T}
-    return 0
-end
-
-Base.length(block::QPBlockData) = length(block.bound_type)
-
-function MOI.set(
-    block::QPBlockData{T},
-    ::MOI.ObjectiveFunction{F},
-    f::F,
-) where {T,F<:Union{MOI.VariableIndex,MOI.ScalarAffineFunction{T}}}
-    block.objective = convert(MOI.ScalarAffineFunction{T}, f)
-    block.objective_function_type = _function_info(f)
-    return
-end
-
-function MOI.set(
-    block::QPBlockData{T},
-    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}},
-    f::MOI.ScalarQuadraticFunction{T},
-) where {T}
-    block.objective = f
-    block.objective_function_type = _function_info(f)
-    return
-end
-
-function MOI.get(block::QPBlockData{T}, ::MOI.ObjectiveFunctionType) where {T}
-    return _function_type_to_func(T, block.objective_function_type)
-end
-
-function MOI.get(block::QPBlockData{T}, ::MOI.ObjectiveFunction{F}) where {T,F}
-    return convert(F, block.objective)
-end
-
-function MOI.get(
-    block::QPBlockData{T},
-    ::MOI.ListOfConstraintTypesPresent,
-) where {T}
-    constraints = Set{Tuple{Type,Type}}()
-    for i in 1:length(block)
-        F = _function_type_to_func(T, block.function_type[i])
-        S = _bound_type_to_set(T, block.bound_type[i])
-        push!(constraints, (F, S))
-    end
-    return collect(constraints)
-end
-
-function MOI.is_valid(
-    block::QPBlockData{T},
-    ci::MOI.ConstraintIndex{F,S},
-) where {
-    T,
-    F<:Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}},
-    S<:Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T},MOI.Interval{T}},
-}
-    return 1 <= ci.value <= length(block)
-end
-
-function MOI.get(
-    block::QPBlockData{T},
-    ::MOI.ListOfConstraintIndices{F,S},
-) where {
-    T,
-    F<:Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}},
-    S<:Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T},MOI.Interval{T}},
-}
-    ret = MOI.ConstraintIndex{F,S}[]
-    for i in 1:length(block)
-        if _bound_type_to_set(T, block.bound_type[i]) != S
-            continue
-        elseif _function_type_to_func(T, block.function_type[i]) != F
-            continue
-        end
-        push!(ret, MOI.ConstraintIndex{F,S}(i))
-    end
-    return ret
-end
-
-function MOI.get(
-    block::QPBlockData{T},
-    ::MOI.NumberOfConstraints{F,S},
-) where {
-    T,
-    F<:Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}},
-    S<:Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T},MOI.Interval{T}},
-}
-    return length(MOI.get(block, MOI.ListOfConstraintIndices{F,S}()))
-end
-
-function MOI.add_constraint(
-    block::QPBlockData{T},
-    f::Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}},
-    s::Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T},MOI.Interval{T}},
-) where {T}
-    push!(block.constraints, f)
-    bound_type, l, u = _set_info(s)
-    push!(block.g_L, l)
-    push!(block.g_U, u)
-    push!(block.mult_g, nothing)
-    push!(block.bound_type, bound_type)
-    push!(block.function_type, _function_info(f))
-    return MOI.ConstraintIndex{typeof(f),typeof(s)}(length(block.bound_type))
-end
-
-function MOI.get(
-    block::QPBlockData{T},
-    ::MOI.ConstraintFunction,
-    c::MOI.ConstraintIndex{F,S},
-) where {T,F,S}
-    return convert(F, block.constraints[c.value])
-end
-
-function MOI.get(
-    block::QPBlockData{T},
-    ::MOI.ConstraintSet,
-    c::MOI.ConstraintIndex{F,S},
-) where {T,F,S}
-    row = c.value
-    if block.bound_type[row] == _kBoundTypeEqualTo
-        return MOI.EqualTo(block.g_L[row])
-    elseif block.bound_type[row] == _kBoundTypeLessThan
-        return MOI.LessThan(block.g_U[row])
-    elseif block.bound_type[row] == _kBoundTypeGreaterThan
-        return MOI.GreaterThan(block.g_L[row])
-    else
-        @assert block.bound_type[row] == _kBoundTypeInterval
-        return MOI.Interval(block.g_L[row], block.g_U[row])
-    end
-end
-
-function MOI.set(
-    block::QPBlockData{T},
-    ::MOI.ConstraintSet,
-    c::MOI.ConstraintIndex{F,MOI.LessThan{T}},
-    set::MOI.LessThan{T},
-) where {T,F}
-    row = c.value
-    block.g_U[row] = set.upper
-    return
-end
-
-function MOI.set(
-    block::QPBlockData{T},
-    ::MOI.ConstraintSet,
-    c::MOI.ConstraintIndex{F,MOI.GreaterThan{T}},
-    set::MOI.GreaterThan{T},
-) where {T,F}
-    row = c.value
-    block.g_L[row] = set.lower
-    return
-end
-
-function MOI.set(
-    block::QPBlockData{T},
-    ::MOI.ConstraintSet,
-    c::MOI.ConstraintIndex{F,MOI.EqualTo{T}},
-    set::MOI.EqualTo{T},
-) where {T,F}
-    row = c.value
-    block.g_L[row] = set.value
-    block.g_U[row] = set.value
-    return
-end
-
-function MOI.set(
-    block::QPBlockData{T},
-    ::MOI.ConstraintSet,
-    c::MOI.ConstraintIndex{F,MOI.Interval{T}},
-    set::MOI.Interval{T},
-) where {T,F}
-    row = c.value
-    block.g_L[row] = set.lower
-    block.g_U[row] = set.upper
-    return
-end
-
-function MOI.get(
-    block::QPBlockData{T},
-    ::MOI.ConstraintDualStart,
-    c::MOI.ConstraintIndex{F,S},
-) where {T,F,S}
-    return block.mult_g[c.value]
-end
-
-function MOI.set(
-    block::QPBlockData{T},
-    ::MOI.ConstraintDualStart,
-    c::MOI.ConstraintIndex{F,S},
-    value,
-) where {T,F,S}
-    block.mult_g[c.value] = value
-    return
-end
-
-function MOI.eval_objective(
-    block::QPBlockData{T},
-    x::AbstractVector{T},
-) where {T}
-    return eval_function(block.objective, x, block.parameters)
-end
-
-function MOI.eval_objective_gradient(
-    block::QPBlockData{T},
-    ∇f::AbstractVector{T},
-    x::AbstractVector{T},
-) where {T}
-    ∇f .= zero(T)
-    eval_dense_gradient(∇f, block.objective, x, block.parameters)
-    return
-end
-
-function MOI.eval_constraint(
-    block::QPBlockData{T},
-    g::AbstractVector{T},
-    x::AbstractVector{T},
-) where {T}
-    for (i, constraint) in enumerate(block.constraints)
-        g[i] = eval_function(constraint, x, block.parameters)
-    end
-    return
-end
-
-function MOI.jacobian_structure(block::QPBlockData)
-    J = Tuple{Int,Int}[]
-    for (row, constraint) in enumerate(block.constraints)
-        append_sparse_gradient_structure!(constraint, J, row)
-    end
-    return J
-end
-
-function MOI.eval_constraint_jacobian(
-    block::QPBlockData{T},
-    J::AbstractVector{T},
-    x::AbstractVector{T},
-) where {T}
-    i = 1
-    for constraint in block.constraints
-        ∇f = view(J, i:length(J))
-        i += eval_sparse_gradient(∇f, constraint, x, block.parameters)
-    end
-    return i
-end
-
-function MOI.hessian_lagrangian_structure(block::QPBlockData)
-    H = Tuple{Int,Int}[]
-    append_sparse_hessian_structure!(block.objective, H)
-    for constraint in block.constraints
-        append_sparse_hessian_structure!(constraint, H)
-    end
-    return H
-end
-
-function MOI.eval_hessian_lagrangian(
-    block::QPBlockData{T},
-    H::AbstractVector{T},
-    x::AbstractVector{T},
-    σ::T,
-    μ::AbstractVector{T},
-) where {T}
-    i = 1
-    i += eval_sparse_hessian(H, block.objective, σ)
-    for (row, constraint) in enumerate(block.constraints)
-        ∇²f = view(H, i:length(H))
-        i += eval_sparse_hessian(∇²f, constraint, μ[row])
-    end
-    return i
-end
diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl
index ccf73eb..a8044db 100644
--- a/test/MOI_wrapper.jl
+++ b/test/MOI_wrapper.jl
@@ -853,7 +853,7 @@ function test_scalar_nonlinear_function_attributes()
 end
 
 function test_vector_nonlinear_oracle()
-    set = Ipopt._VectorNonlinearOracle(;
+    set = MOI.VectorNonlinearOracle(;
         dimension = 5,
         l = [0.0, 0.0],
         u = [0.0, 0.0],
@@ -887,7 +887,7 @@ function test_vector_nonlinear_oracle()
         end,
     )
     @test MOI.dimension(set) == 5
-    @test MOI.copy(set) === set
+    @test MOI.copy(set) == set
     model = Ipopt.Optimizer()
     MOI.set(model, MOI.Silent(), true)
     x = MOI.add_variables(model, 3)
@@ -895,7 +895,7 @@ function test_vector_nonlinear_oracle()
     y = MOI.add_variables(model, 2)
     MOI.optimize!(model)
     f = MOI.VectorOfVariables([x; y])
-    F, S = MOI.VectorOfVariables, Ipopt._VectorNonlinearOracle
+    F, S = MOI.VectorOfVariables, MOI.VectorNonlinearOracle{Float64}
     @test MOI.supports_constraint(model, F, S)
     @test !((F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()))
     @test isempty(MOI.get(model, MOI.ListOfConstraintIndices{F,S}()))
@@ -919,7 +919,7 @@ function test_vector_nonlinear_oracle()
 end
 
 function test_vector_nonlinear_oracle_two()
-    set = Ipopt._VectorNonlinearOracle(;
+    set = MOI.VectorNonlinearOracle(;
         dimension = 5,
         l = [0.0, 0.0],
         u = [0.0, 0.0],
@@ -958,7 +958,7 @@ function test_vector_nonlinear_oracle_two()
     c_z = MOI.add_constraint(model, f_z, set)
     @test MOI.is_valid(model, c_y)
     @test MOI.is_valid(model, c_z)
-    F, S = MOI.VectorOfVariables, Ipopt._VectorNonlinearOracle
+    F, S = MOI.VectorOfVariables, MOI.VectorNonlinearOracle{Float64}
     @test (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
     @test MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) == [c_y, c_z]
     @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 2
@@ -976,7 +976,7 @@ function test_vector_nonlinear_oracle_two()
 end
 
 function test_vector_nonlinear_oracle_optimization()
-    set = Ipopt._VectorNonlinearOracle(;
+    set = MOI.VectorNonlinearOracle(;
         dimension = 4,
         l = [0.0, 0.0],
         u = [0.0, 0.0],
@@ -1037,7 +1037,7 @@ function test_vector_nonlinear_oracle_optimization()
 end
 
 function test_vector_nonlinear_oracle_optimization_min_sense()
-    set = Ipopt._VectorNonlinearOracle(;
+    set = MOI.VectorNonlinearOracle(;
         dimension = 4,
         l = [0.0, 0.0],
         u = [0.0, 0.0],
@@ -1135,7 +1135,7 @@ function test_vector_nonlinear_oracle_scalar_nonlinear_equivalent()
 end
 
 function test_vector_nonlinear_oracle_no_hessian()
-    set = Ipopt._VectorNonlinearOracle(;
+    set = MOI.VectorNonlinearOracle(;
         dimension = 5,
         l = [0.0, 0.0],
         u = [0.0, 0.0],