removed: old dead code

Author: franckgaga

src/estimator/mhe/construct.jl

@@ -1280,20 +1280,17 @@ function init_optimization!(
         end
     end
     # constraints with vector nonlinear oracle, objective function with splatting:
-    g_oracle, J_func, ∇J_func! = get_nonlinops(estim, optim)
-    @operator(optim, J, nZ̃, J_func, ∇J_func!)
-    @objective(optim, Min, J(Z̃var...))
+    g_oracle, J_op = get_nonlinops(estim, optim)
+    optim[:J_op] = J_op
+    @objective(optim, Min, J_op(Z̃var...))
     set_nonlincon!(estim, model, optim, g_oracle)
     return nothing
 end
 
-
 """
-    get_optim_functions(
-        estim::MovingHorizonEstimator, optim::JuMP.GenericModel,
-    ) -> Jfunc, ∇Jfunc!, gfuncs, ∇gfuncs!
+    get_nonlinops(estim::MovingHorizonEstimator, optim) -> g_oracle, J_op
 
-Return the functions for the nonlinear optimization of [`MovingHorizonEstimator`](@ref).
+Return the operators for the nonlinear optimization of [`MovingHorizonEstimator`](@ref).
 
 Return `g_oracle`, the [`VectorNonlinearOracle`](@extref MOI MathOptInterface.VectorNonlinearOracle)
 for the ineqaulity constraints. Note that `g_oracle` only includes the non-`Inf`
@@ -1308,117 +1305,8 @@ on the splatting syntax. This method is really intricate and that's because of 2
   and memoization to avoid redundant computations. This is already complex, but it's even
   worse knowing that the automatic differentiation tools do not support splatting.
 """
-function get_optim_functions(
-    estim::MovingHorizonEstimator, ::JuMP.GenericModel{JNT},
-) where {JNT <: Real}
-    # ----------- common cache for Jfunc and gfuncs  --------------------------------------
-    model, con = estim.model, estim.con
-    grad, jac = estim.gradient, estim.jacobian
-    nx̂, nym, nŷ, nu, nε, nk = estim.nx̂, estim.nym, model.ny, model.nu, estim.nε, model.nk
-    He = estim.He
-    nV̂, nX̂, ng, nZ̃ = He*nym, He*nx̂, length(con.i_g), length(estim.Z̃)
-    strict = Val(true)
-    myNaN = convert(JNT, NaN)
-    J::Vector{JNT}                   = zeros(JNT, 1)
-    V̂::Vector{JNT},  X̂0::Vector{JNT} = zeros(JNT, nV̂), zeros(JNT, nX̂)
-    k0::Vector{JNT}                  = zeros(JNT, nk)
-    û0::Vector{JNT}, ŷ0::Vector{JNT} = zeros(JNT, nu), zeros(JNT, nŷ)
-    g::Vector{JNT}                   = zeros(JNT, ng)
-    x̄::Vector{JNT}                   = zeros(JNT, nx̂)
-    # --------------------- objective functions -------------------------------------------
-    function Jfunc!(Z̃, V̂, X̂0, û0, k0, ŷ0, g, x̄)
-        update_prediction!(V̂, X̂0, û0, k0, ŷ0, g, estim, Z̃)
-        return obj_nonlinprog!(x̄, estim, model, V̂, Z̃)
-    end
-    Z̃_∇J = fill(myNaN, nZ̃)      # NaN to force update_predictions! at first call
-    ∇J_context = (
-        Cache(V̂),  Cache(X̂0), Cache(û0), Cache(k0), Cache(ŷ0),
-        Cache(g),
-        Cache(x̄),
-    )
-    # temporarily "fill" the estimation window for the preparation of the gradient: 
-    estim.Nk[] = He
-    ∇J_prep = prepare_gradient(Jfunc!, grad, Z̃_∇J, ∇J_context...; strict)
-    estim.Nk[] = 0
-    ∇J = Vector{JNT}(undef, nZ̃)
-    function update_objective!(J, ∇J, Z̃_∇J, Z̃arg)
-        if isdifferent(Z̃arg, Z̃_∇J)
-            Z̃_∇J .= Z̃arg
-            J[], _ = value_and_gradient!(Jfunc!, ∇J, ∇J_prep, grad, Z̃_∇J, ∇J_context...)
-        end
-    end
-    function J_func(Z̃arg::Vararg{T, N}) where {N, T<:Real}
-        update_objective!(J, ∇J, Z̃_∇J, Z̃arg)
-        return J[]::T
-    end
-    ∇J_func! = function (∇Jarg::AbstractVector{T}, Z̃arg::Vararg{T, N}) where {N, T<:Real}
-        # only the multivariate syntax of JuMP.@operator, univariate is impossible for MHE
-        # since Z̃ comprises the arrival state estimate AND the estimated process noise
-        update_objective!(J, ∇J, Z̃_∇J, Z̃arg)
-        return ∇Jarg .= ∇J
-    end
-    # --------------------- inequality constraint functions -------------------------------
-    function gfunc!(g, Z̃, V̂, X̂0, û0, k0, ŷ0)
-        return update_prediction!(V̂, X̂0, û0, k0, ŷ0, g, estim, Z̃)
-    end
-    Z̃_∇g = fill(myNaN, nZ̃)      # NaN to force update_predictions! at first call
-    ∇g_context = (
-        Cache(V̂), Cache(X̂0), Cache(û0), Cache(k0), Cache(ŷ0),
-    )
-    # temporarily enable all the inequality constraints for sparsity detection:
-    estim.con.i_g .= true  
-    estim.Nk[] = He
-    ∇g_prep  = prepare_jacobian(gfunc!, g, jac, Z̃_∇g, ∇g_context...; strict)
-    estim.con.i_g .= false
-    estim.Nk[] = 0
-    ∇g = init_diffmat(JNT, jac, ∇g_prep, nZ̃, ng)
-    function update_con!(g, ∇g, Z̃_∇g, Z̃arg)
-        if isdifferent(Z̃arg, Z̃_∇g)
-            Z̃_∇g .= Z̃arg
-            value_and_jacobian!(gfunc!, g, ∇g, ∇g_prep, jac, Z̃_∇g, ∇g_context...)
-        end
-    end
-    g_funcs = Vector{Function}(undef, ng)
-    for i in eachindex(g_funcs)
-        gfunc_i = function (Z̃arg::Vararg{T, N}) where {N, T<:Real}
-            update_con!(g, ∇g, Z̃_∇g, Z̃arg)
-            return g[i]::T
-        end
-        g_funcs[i] = gfunc_i
-    end
-    ∇g_funcs! = Vector{Function}(undef, ng)
-    for i in eachindex(∇g_funcs!)
-        ∇g_funcs![i] = function (∇g_i, Z̃arg::Vararg{T, N}) where {N, T<:Real}
-            # only the multivariate syntax of JuMP.@operator, see above for the explanation
-            update_con!(g, ∇g, Z̃_∇g, Z̃arg)
-            return ∇g_i .= @views ∇g[i, :]
-        end
-    end
-    return J_func, ∇J_func!, g_funcs, ∇g_funcs!
-end
-
-"""
-    get_nonlinops(estim::MovingHorizonEstimator, optim) -> g_oracle, J_op
-
-Return the operators for the nonlinear optimization of [`MovingHorizonEstimator`](@ref).
-
-Return `g_oracle`, the inequality [`VectorNonlinearOracle`](@extref MathOptInterface MathOptInterface.VectorNonlinearOracle),
-. Note that `g_oracle` only includes the non-`Inf` inequality constraints, thus
-it must be re-constructed if they change. Also return `J_op`, the [`NonlinearOperator`](@extref JuMP NonlinearOperator)
-for the objective function, based on the splatting syntax. This method is really intricate
-and that's because of 3 elements:
-
-- These functions are used inside the nonlinear optimization, so they must be type-stable
-  and as efficient as possible. All the function outputs and derivatives are cached and
-  updated in-place if required to use the efficient [`value_and_jacobian!`](@extref DifferentiationInterface DifferentiationInterface.value_and_jacobian!).
-- The splatting syntax for objective functions implies the use of `Vararg{T,N}` (see the [performance tip](@extref Julia Be-aware-of-when-Julia-avoids-specializing))
-  and memoization to avoid redundant computations. This is already complex, but it's even
-  worse knowing that the automatic differentiation tools do not support splatting.
-- The signature of gradient and hessian functions is not the same for univariate (`nZ̃ == 1`)
-  and multivariate (`nZ̃ > 1`) operators in `JuMP`. Both must be defined.
-"""
 function get_nonlinops(
-    estim::MovingHorizonEstimator, ::JuMP.GenericModel{JNT}
+    estim::MovingHorizonEstimator, optim::JuMP.GenericModel{JNT}
 ) where JNT<:Real
     # ----------- common cache for Jfunc and gfuncs  --------------------------------------
     model, con = estim.model, estim.con
@@ -1503,18 +1391,14 @@ function get_nonlinops(
         update_objective!(J, ∇J, Z̃_∇J, Z̃_arg)
         return J[]::T
     end
-    ∇J_func! = if nZ̃ == 1        # univariate syntax (see JuMP.@operator doc):
-        function (Z̃_arg)
-            update_objective!(J, ∇J, Z̃_∇J, Z̃_arg)
-            return ∇J[]
-        end
-    else                        # multivariate syntax (see JuMP.@operator doc):
-        function (∇J_arg::AbstractVector{T}, Z̃_arg::Vararg{T, N}) where {N, T<:Real}
-            update_objective!(J, ∇J, Z̃_∇J, Z̃_arg)
-            return ∇J_arg .= ∇J
-        end
+    function ∇J_func!(∇J_arg::AbstractVector{T}, Z̃_arg::Vararg{T, N}) where {N, T<:Real}
+        # only the multivariate syntax of JuMP.@operator, univariate is impossible for MHE
+        # since Z̃ comprises the arrival state estimate AND the estimated process noise
+        update_objective!(J, ∇J, Z̃_∇J, Z̃_arg)
+        return ∇J_arg .= ∇J
     end
-    g_oracle, J_func, ∇J_func!
+    J_op = JuMP.add_nonlinear_operator(optim, nZ̃, J_func, ∇J_func!, name=:J_op)
+    g_oracle, J_op
 end
 
 "By default, there is no nonlinear constraint, thus do nothing."