text:
# TODO: have just removed some converts() because the promote_type obviously doesn't work when input empty

# Note: Might need to relax the type on mutation arguments for the mm_ functions, not the exported functions though.
# Might also need to relax the typing on arguments for the local functions in mm_ functions.

# Need to ensure convert() is used consistently

"""
    mpccmodel_setup_forwarddiff_dense(config::MPCCModelConfig)

Accepts model config and returns a struct with all the nice Julia functions to
calculate Jacobians, Hessians, and parametric derivatives using
ForwardDiff.jl.

Both standard and mutating functions are constructed. Most of the functions
accept an index as the last argument to specify a subset of the constraints.
For F constraints, this should be a vector of either tuple pairs or
CartesianIndexes.
"""
function mpccmodel_setup_forwarddiff_dense(config::MPCCModelConfig)

    @unpack dimspec, fns = config

    # Full f functions
    function local_f(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return fns.f(x, pr, ps)[1]
    end

    function local_f!(out_f::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        # fns.f!(@view(out_f[1]), x, pr, ps)
        fns.f!(out_f, x, pr, ps)
        return nothing
    end

    # Full ce functions
    function local_ce(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        # return convert(Vector{promote_type(S, T)}, fns.ce(x, pr, ps))
        return fns.ce(x, pr, ps)
    end

    function local_ce!(out_ce::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        # fns.ce!(@view(out_ce[:]), x, pr, ps)
        fns.ce!(out_ce, x, pr, ps)
        return nothing
    end

    # Indexed ce functions
    function local_ce(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_ce_i(dimspec, fns.ce_i, x, pr, ps, idxs_ce)
    end

    function local_ce!(out_ce_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_ce_i!(out_ce_i, dimspec, fns.ce_i!, x, pr, ps, idxs_ce)
        return nothing
    end

    # Full ci functions
    function local_ci(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        # return convert(Vector{promote_type(S, T)}, fns.ci(x, pr, ps))
        return fns.ci(x, pr, ps)
    end

    function local_ci!(out_ci::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        # fns.ci!(@view(out_ci[:]), x, pr, ps)
        fns.ci!(out_ci, x, pr, ps)
        return nothing
    end

    # Indexed ci functions
    function local_ci(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_ci_i(dimspec, fns.ci_i, x, pr, ps, idxs_ci)
    end

    function local_ci!(out_ci_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_ci!(out_ci_i, dimspec, fns.ci_i!, x, pr, ps, idxs_ci)
        return nothing
    end

    # Full F functions
    function local_F(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}  # ::Matrix{S}
        # return convert(Matrix{promote_type(S, T)}, fns.F(x, pr, ps))
        return fns.F(x, pr, ps)
    end

    function local_F!(out_F::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        # fns.F!(@view(out_F[:,:]), x, pr, ps)
        fns.F!(out_F, x, pr, ps)
        return nothing
    end

    # Indexed F functions; NOTE, these yeild a vector in the order that the (l,q) index tuples were in!
    function local_F(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {S <: Real, T <: Real}
        return mm_fd_dn_F_i(dimspec, fns.F_i, x, pr, ps, idxs_F)
    end

    function local_F(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}}) where {S <: Real, T <: Real}
        cart_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        return mm_fd_dn_F_i(dimspec, fns.F_i, x, pr, ps, cart_idxs_F)
    end

    function local_F!(out_F_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_F_i!(out_F_i, dimspec, fns.F_i!, x, pr, ps, idxs_F)
        return nothing
    end

    function local_F!(out_F_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}})::Nothing where {S <: Real, T <: Real}
        cart_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        mm_fd_dn_F_i!(out_F_i, dimspec, fns.F_i!, x, pr, ps, cart_idxs_F)
        return nothing
    end

    # Full Fq functions
    function local_Fq(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        # return convert(Vector{promote_type(S, T)}, fns.Fq(x, pr, ps))
        return fns.Fq(x, pr, ps)
    end

    function local_Fq!(out_Fq::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        # fns.Fq!(@view(out_Fq[:]), x, pr, ps)
        fns.Fq!(out_Fq, x, pr, ps)
        return nothing
    end

    # Indexed Fq functions
    function local_Fq(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_Fq_i(dimspec, fns.Fq_i, x, pr, ps, idxs_Fq)
    end

    function local_Fq!(out_Fq_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_Fq!(out_Fq_i, dimspec, fns.Fq_i!, x, pr, ps, idxs_Fq)
        return nothing
    end

    # gradf functions
    function local_gradf(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradf(dimspec, fns.f, x, pr, ps)
    end

    function local_gradf!(out_gradf::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradf!(out_gradf, dimspec, fns.f, x, pr, ps)
        return nothing
    end

    # Full gradce functions
    function local_gradce(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradce(dimspec, fns.ce, x, pr, ps)
    end

    function local_gradce!(out_gradce::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradce!(out_gradce, dimspec, fns.ce!, x, pr, ps)
        return nothing
    end

    # Indexed gradce functions
    function local_gradce(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradce_i(dimspec, fns.ce_i, x, pr, ps, idxs_ce)
    end

    function local_gradce!(out_gradce_i::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradce_i!(out_gradce_i, dimspec, fns.ce_i!, x, pr, ps, idxs_ce)
        return nothing
    end

    # Full gradci functions
    function local_gradci(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradci(dimspec, fns.ci, x, pr, ps)
    end

    function local_gradci!(out_gradci::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradci!(out_gradci, dimspec, fns.ci!, x, pr, ps)
        return nothing
    end

    # Indexed gradci functions
    function local_gradci(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradci_i(dimspec, fns.ci_i, x, pr, ps, idxs_ci)
    end

    function local_gradci!(out_gradci_i::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradci_i!(out_gradci_i, dimspec, fns.ci_i!, x, pr, ps, idxs_ci)
        return nothing
    end

    # Full gradF functions
    function local_gradF(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradF(dimspec, fns.F, x, pr, ps)
    end

    function local_gradF!(out_gradF::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradF!(out_gradF, dimspec, fns.F!, x, pr, ps)
        return nothing
    end

    # Indexed gradF functions
    # NOTE, these yeild a vector in the order that the (l,q) index tuples were in!
    function local_gradF(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradF_i(dimspec, fns.F_i, x, pr, ps, idxs_F)
    end

    function local_gradF(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}}) where {S <: Real, T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        return mm_fd_dn_gradF_i(dimspec, fns.F_i, x, pr, ps, new_idxs_F)
    end

    function local_gradF!(out_gradF_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradF_i!(out_gradF_i, dimspec, fns.F_i!, x, pr, ps, idxs_F)
        return nothing
    end

    function local_gradF!(out_gradF_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}})::Nothing where {S <: Real, T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        mm_fd_dn_gradF_i!(out_gradF_i, dimspec, fns.F_i!, x, pr, ps, new_idxs_F)
        return nothing
    end

    # Full gradFq functions
    function local_gradFq(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradFq(dimspec, fns.Fq, x, pr, ps)
    end

    function local_gradFq!(out_gradF::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradFq!(out_gradF, dimspec, fns.Fq!, x, pr, ps)
        return nothing
    end

    # Indexed gradF functions
    function local_gradFq(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradFq_i(dimspec, fns.Fq_i, x, pr, ps, idxs_Fq)
    end

    function local_gradFq!(out_gradFq_i::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradFq_i!(out_gradFq_i, dimspec, fns.Fq_i!, x, pr, ps, idxs_Fq)
        return nothing
    end

    # Hessian f functions (we don't do indexed version for this)
    function local_hessf(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessf(dimspec, fns.f, x, pr, ps)
    end

    function local_hessf!(out_hessf::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_hessf!(out_hessf, dimspec, fns.f, x, pr, ps)
        return nothing
    end

    # Full hessce functions
        # Vector of length me of 2d Hessian matrices
    function local_hessce(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessce(dimspec, fns.ce_i, x, pr, ps)
    end

    # function local_hessce!(out_hessce::AbstractVector{AbstractMatrix{R}}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {R <: Real, S <: Real, T <: Real}
    function local_hessce!(out_hessce::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_hessce!(out_hessce, dimspec, fns.ce_i, x, pr, ps)
        return nothing
    end

    # Indexed hessce functions
    function local_hessce(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessce_i(dimspec, fns.ce_i, x, pr, ps, idxs_ce)
    end

    function local_hessce!(out_hessce_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_hessce_i!(out_hessce_i, dimspec, fns.ce_i, x, pr, ps, idxs_ce)
        return nothing
    end

    # Full hessci functions
    # Vector of length me of 2d Hessian matricis
    function local_hessci(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessci(dimspec, fns.ci_i, x, pr, ps)
    end

    function local_hessci!(out_hessci::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_hessci!(out_hessci, dimspec, fns.ci_i, x, pr, ps)
        return nothing
    end

    # Indexed hessci functions
    function local_hessci(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessci_i(dimspec, fns.ci_i, x, pr, ps, idxs_ci)
    end

    function local_hessci!(out_hessci_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_hessci_i!(out_hessci_i, dimspec, fns.ci_i, x, pr, ps, idxs_ci)
        return nothing
    end

    # Full hessF functions
    function local_hessF(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}  # ::Matrix{Matrix{S}}
        return mm_fd_dn_hessF(dimspec, fns.F_i, x, pr, ps)
    end

    function local_hessF!(out_hessF::AbstractMatrix, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_hessF!(out_hessF, dimspec, fns.F_i, x, pr, ps)
        return nothing
    end

    # Indexed hessF functions
    # NOTE, these yeild a vector in the order that the (l,q) index tuples were in!
    function local_hessF(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessF_i(dimspec, fns.F_i, x, pr, ps, idxs_F)
    end

    function local_hessF(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}}) where {S <: Real, T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        return mm_fd_dn_hessF_i(dimspec, fns.F_i, x, pr, ps, new_idxs_F)
    end

    function local_hessF!(out_hessF_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_hessF_i!(out_hessF_i, dimspec, fns.F_i, x, pr, ps, idxs_F)
        return nothing
    end

    function local_hessF!(out_hessF_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}})::Nothing where {S <: Real, T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        mm_fd_dn_hessF_i!(out_hessF_i, dimspec, fns.F_i, x, pr, ps, new_idxs_F)
        return nothing
    end

    # Full hessFq functions

    function local_hessFq(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessFq(dimspec, fns.Fq_i, x, pr, ps)
    end

    function local_hessFq!(out_hessFq::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        # mm_fd_dn_hessFq!(out_hessFq, dimspec, fns.Fq_i, x, pr, ps)
        error("not implemented")
        return nothing
    end

    # Indexed hessFq functions
    function local_hessFq(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_hessFq_i(dimspec, fns.Fq_i, x, pr, ps, idxs_Fq)
    end

    function local_hessFq!(out_hessFq_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        # mm_fd_dn_hessFq_i!(out_hessFq_i, dimspec, fns.Fq_i, x, pr, ps, idxs_Fq)
        error("not implemented")
        return nothing
    end

    # fdp functions
    function local_fdp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}   # ::AbstractVector{T}
        return mm_fd_dn_fdp(dimspec, fns.f, x, pr, ps)
    end

    function local_fdp!(out_fdp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_fdp!(out_fdp, dimspec, fns.f, x, pr, ps)
        return nothing
    end

    # Full cedp functions
    function local_cedp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}  # ::Vector{Vector{S}}
        return mm_fd_dn_cedp(dimspec, fns.ce, x, pr, ps)
    end

    function local_cedp!(out_cedp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_cedp!(out_cedp, dimspec, fns.ce!, x, pr, ps)
        return nothing
    end

    # Indexed cedp functions
    function local_cedp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {S <: Real, T <: Real}  # Vector{Vector} ret vector of len r, of vectors of len idx
        return mm_fd_dn_cedp_i(dimspec, fns.ce_i, x, pr, ps, idxs_ce)
    end

    function local_cedp!(out_cedp_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_cedp_i!(out_cedp_i, dimspec, fns.ce_i!, x, pr, ps, idxs_ce)
        return nothing
    end

    # Full cidp functions
    function local_cidp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}  # ::Vector{Vector{S}}
        return mm_fd_dn_cidp(dimspec, fns.ci, x, pr, ps)
    end

    function local_cidp!(out_cidp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_cidp!(out_cidp, dimspec, fns.ci!, x, pr, ps)
        return nothing
    end

    # Indexed cidp functions
    function local_cidp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {S <: Real, T <: Real}  # Vector{Vector} ret vector of len r, of vectors of len idx
        return mm_fd_dn_cidp_i(dimspec, fns.ci_i, x, pr, ps, idxs_ci)
    end

    function local_cidp!(out_cidp_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_cidp_i!(out_cidp_i, dimspec, fns.ci_i!, x, pr, ps, idxs_ci)
        return nothing
    end

    # Full Fdp functions
    function local_Fdp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}   # ::Vector{Matrix{T}}
        return mm_fd_dn_Fdp(dimspec, fns.F, x, pr, ps)
    end

    function local_Fdp!(out_Fdp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_Fdp!(out_Fdp, dimspec, fns.F!, x, pr, ps)
        return nothing
    end

    # Indexed Fdp functions
    # NOTE, these yeild a vector in the order that the (l,q) index tuples were in!
    function local_Fdp(x::AbstractVector{T}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Vector{Vector{T}} where {T <: Real}
        return mm_fd_dn_Fdp_i(dimspec, fns.F_i, x, pr, ps, idxs_F)
    end

    function local_Fdp(x::AbstractVector{T}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}})::Vector{Vector{T}} where {T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        return mm_fd_dn_Fdp_i(dimspec, fns.F_i, x, pr, ps, new_idxs_F)
    end

    function local_Fdp!(out_Fdp_i::AbstractVector, x::AbstractVector{T}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {T <: Real}
        mm_fd_dn_Fdp_i!(out_Fdp_i, dimspec, fns.F_i!, x, pr, ps, idxs_F)
        return nothing
    end

    function local_Fdp!(out_Fdp_i::AbstractVector, x::AbstractVector{T}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}})::Nothing where {T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        mm_fd_dn_Fdp_i!(out_Fdp_i, dimspec, fns.F_i!, x, pr, ps, new_idxs_F)
        return nothing
    end

    # Full gradfdp functions
    function local_gradfdp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}   # ::Vector{Vector{T}}
        return mm_fd_dn_gradfdp(dimspec, fns.f, x, pr, ps)
    end

    function local_gradfdp!(out_Fdp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradfdp!(out_Fdp, dimspec, fns.f, x, pr, ps)
        return nothing
    end

    # Full gradcedp functions
    function local_gradcedp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}  # ::Vector{Matrix{S}}
        return mm_fd_dn_gradcedp(dimspec, fns.ce, x, pr, ps)
    end

    function local_gradcedp!(out_gradcedp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradcedp!(out_gradcedp, dimspec, fns.ce!, x, pr, ps)
        return nothing
    end

    # Indexed gradcedp functions
    function local_gradcedp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Vector{Matrix{T}} where {S <: Real, T <: Real}
        return mm_fd_dn_gradcedp_i(dimspec, fns.ce_i, x, pr, ps, idxs_ce)
    end

    function local_gradcedp!(out_gradcedp_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradcedp_i!(out_gradcedp_i, dimspec, fns.ce_i!, x, pr, ps, idxs_ce)
        return nothing
    end

    # Full gradcidp functions
    function local_gradcidp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradcidp(dimspec, fns.ci, x, pr, ps)
    end

    function local_gradcidp!(out_gradcidp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradcidp!(out_gradcidp, dimspec, fns.ci!, x, pr, ps)
        return nothing
    end

    # Indexed gradcidp functions
    function local_gradcidp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Vector{Matrix{S}} where {S <: Real, T <: Real}
        return mm_fd_dn_gradcidp_i(dimspec, fns.ci_i, x, pr, ps, idxs_ci)
    end

    function local_gradcidp!(out_gradcidp_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradcidp_i!(out_gradcidp_i, dimspec, fns.ci_i!, x, pr, ps, idxs_ci)
        return nothing
    end

    # Full gradFdp functions
    function local_gradFdp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradFdp(dimspec, fns.F_i, x, pr, ps)
    end

    function local_gradFdp!(out_gradFdp::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradFdp!(out_gradFdp, dimspec, fns.F_i, x, pr, ps)
        return nothing
    end

    # Indexed gradFdp functions
    function local_gradFdp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {S <: Real, T <: Real}
        return mm_fd_dn_gradFdp_i(dimspec, fns.F_i, x, pr, ps, idxs_F)
    end

    function local_gradFdp(x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}}) where {S <: Real, T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        return mm_fd_dn_gradFdp_i(dimspec, fns.F_i, x, pr, ps, new_idxs_F)
    end

    function local_gradFdp!(out_gradFdp_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {S <: Real, T <: Real}
        mm_fd_dn_gradFdp_i!(out_gradFdp_i, dimspec, fns.F_i, x, pr, ps, idxs_F)
        return nothing
    end

    function local_gradFdp!(out_gradFdp_i::AbstractVector, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}})::Nothing where {S <: Real, T <: Real}
        new_idxs_F = map(x -> CartesianIndex{2}(x), idxs_F)
        mm_fd_dn_gradFdp_i!(out_gradFdp_i, dimspec, fns.F_i, x, pr, ps, new_idxs_F)
        return nothing
    end

    return MPCCModel(   config,
                        local_f, local_f!,
                        local_ce, local_ce!,
                        local_ci, local_ci!,
                        local_F, local_F!,
                        local_Fq, local_Fq!,
                        local_gradf, local_gradf!,
                        local_gradce, local_gradce!,
                        local_gradci, local_gradci!,
                        local_gradF, local_gradF!,
                        local_gradFq, local_gradFq!,
                        local_hessf, local_hessf!,
                        local_hessce, local_hessce!,
                        local_hessci, local_hessci!,
                        local_hessF, local_hessF!,
                        local_hessFq, local_hessFq!,
                        local_fdp, local_fdp!,
                        local_cedp, local_cedp!,
                        local_cidp, local_cidp!,
                        local_Fdp, local_Fdp!,
                        local_gradfdp, local_gradfdp!,
                        local_gradcedp, local_gradcedp!,
                        local_gradcidp, local_gradcidp!,
                        local_gradFdp, local_gradFdp!
                        # local_do_self_test,
                        # nz_mask
                    )
end

# Indexed ce functions
function mm_fd_dn_ce_i(dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    res_ce_i = [ ce_i[idxs_ce[lp_ce]](x, pr, ps)[1] for lp_ce in eachindex(idxs_ce) ]
    return res_ce_i
end

function mm_fd_dn_ce_i!(out_ce_i::AbstractArray, dimspec::MPCCDimSpec, ce_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    for lp_ce in eachindex(idxs_ce)
        ce_i![lp_ce](@view(out_ce_i[lp_ce:lp_ce]), x, pr, ps)        # As first argument we pass a 1-element vector, so that it can be mutated
    end
    return nothing
end

# Indexed ci functions
function mm_fd_dn_ci_i(dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    res_ci_i = [ ci_i[idxs_ci[lp_ci]](x, pr, ps)[1] for lp_ci in eachindex(idxs_ci) ]
    return res_ci_i
end

function mm_fd_dn_ci_i!(out_ci_i::AbstractArray, dimspec::MPCCDimSpec, ci_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    for lp_ci in eachindex(idxs_ci)
        ci_i![lp_ci](@view(out_ci_i[lp_ci:lp_ci]), x, pr, ps)        # As first argument we pass a 1-element vector, so that it can be mutated
    end
    return nothing
end

# Indexed F functions; NOTE, these yeild a vector in the order that the (l,q) index tuples were in!
function mm_fd_dn_F_i(dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {F <: Function, S <: Real, T <: Real}
    res_F_i = [ F_i[idxs_F[lp_F]](x, pr, ps)[1] for lp_F in eachindex(idxs_F) ]
    return res_F_i
end

function mm_fd_dn_F_i!(out_F_i::AbstractArray, dimspec::MPCCDimSpec, F_i!::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {F <: Function, S <: Real, T <: Real}
    for lp_F in eachindex(idxs_F)
        # This is a bawhair away from being a contestant in an "obfuscated Julia" competition...
        F_i![idxs_F[lp_F]](@view(out_F_i[lp_F:lp_F]), x, pr, ps)
    end
    return nothing
end

# Indexed Fq functions
function mm_fd_dn_Fq_i(dimspec::MPCCDimSpec, Fq_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    res_Fq_i = [ Fq_i[idxs_Fq[lp_Fq]](x, pr, ps)[1] for lp_Fq in eachindex(idxs_Fq) ]
    return res_Fq_i
end

function mm_fd_dn_Fq!(out_Fq_i::AbstractArray, dimspec::MPCCDimSpec, Fq_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    for lp_Fq in eachindex(idxs_Fq)
        Fq_i![lp_Fq](@view(out_Fq_i[lp_Fq:lp_Fq]), x, pr, ps)        # As first argument we pass a 1-element vector, so that it can be mutated
    end
    return nothing
end

# gradf functions

function mm_fd_dn_gradf(dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    local_f(z) = f(z, pr, ps)[1]
    return ForwardDiff.gradient(local_f, x)
end

function mm_fd_dn_gradf!(out_gradf::AbstractArray, dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    local_f(z) = f(z, pr, ps)[1]
    ForwardDiff.gradient!(out_gradf, local_f, x)
    return nothing
end

# Full gradce functions
function mm_fd_dn_gradce(dimspec::MPCCDimSpec, ce::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    local_ce = (z::AbstractArray) -> ce(z, pr, ps)
    # return convert(Matrix{promote_type(S, T)}, ForwardDiff.jacobian(local_ce, x))
    return ForwardDiff.jacobian(local_ce, x)
end

function mm_fd_dn_gradce!(out_gradce::AbstractArray, dimspec::MPCCDimSpec, ce!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, me = dimspec
    local_ce! = (y::AbstractArray, z::AbstractArray) -> ce!(y, z, pr, ps)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), me)
    ForwardDiff.jacobian!(out_gradce, local_ce!, y, x)
    return nothing
end

# Indexed gradce functions
function mm_fd_dn_gradce_i(dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    local_ce = (z::AbstractArray) -> mm_fd_dn_ce_i(dimspec, ce_i, z, pr, ps, idxs_ce)
    # return convert(Matrix{promote_type(S, T)}, ForwardDiff.jacobian(local_ce, x))
    return ForwardDiff.jacobian(local_ce, x)
end

function mm_fd_dn_gradce_i!(out_gradce_i::AbstractArray, dimspec::MPCCDimSpec, ce_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n = dimspec
    len_idxs_ce_i = length(idxs_ce)
    local_ce! = (y::AbstractArray, z::AbstractArray) -> mm_fd_dn_ce_i!(y, dimspec, ce_i!, z, pr, ps, idxs_ce)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), len_idxs_ce_i)
    ForwardDiff.jacobian!(out_gradce_i, local_ce!, y, x)
    return nothing
end

# Full gradci functions
function mm_fd_dn_gradci(dimspec::MPCCDimSpec, ci::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    local_ci = (z::AbstractArray) -> ci(z, pr, ps)
    # return convert(Matrix{promote_type(S, T)}, ForwardDiff.jacobian(local_ci, x))
    return ForwardDiff.jacobian(local_ci, x)
end

function mm_fd_dn_gradci!(out_gradci::AbstractArray, dimspec::MPCCDimSpec, ci!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, mi = dimspec
    local_ci! = (y::AbstractArray, z::AbstractArray) -> ci!(y, z, pr, ps)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), mi)
    ForwardDiff.jacobian!(out_gradci, local_ci!, y, x)
    return nothing
end

# Indexed gradci functions
function mm_fd_dn_gradci_i(dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    local_ci = (z::AbstractArray) -> mm_fd_dn_ci_i(dimspec, ci_i, z, pr, ps, idxs_ci)
    # return convert(Matrix{promote_type(S, T)}, ForwardDiff.jacobian(local_ci, x))
    return ForwardDiff.jacobian(local_ci, x)
end

function mm_fd_dn_gradci_i!(out_gradci_i::AbstractArray, dimspec::MPCCDimSpec, ci_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n = dimspec
    len_idxs_ci_i = length(idxs_ci)
    local_ci! = (y::AbstractArray, z::AbstractArray) -> mm_fd_dn_ci_i!(y, dimspec, ci_i!, z, pr, ps, idxs_ci)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), len_idxs_ci_i)
    ForwardDiff.jacobian!(out_gradci_i, local_ci!, y, x)
    return nothing
end

# Full gradF functions
# # Matrix of dim l x q of n length vectors (the gradients wrt x)

function mm_fd_dn_gradF(dimspec::MPCCDimSpec, F::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}   # ::Matrix{Vector{S}}
    @unpack n, l, q = dimspec
    # ForwardDiff.jacobian can work with multidimensional arrays, so we do that
    local_F = (z::AbstractArray) -> F(z, pr, ps)
    gradF_flat = ForwardDiff.jacobian(local_F, x)
    gradF = [ gradF_flat[l*(lp_q-1)+lp_l, :] for lp_l in 1:l, lp_q in 1:q ]
    return gradF
end

function mm_fd_dn_gradF!(out_gradF::AbstractArray, dimspec::MPCCDimSpec, F!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, l, q = dimspec
    # ForwardDiff.jacobian can work with multidimensional arrays, so we do that
    local_F! = (y::AbstractArray, z::AbstractArray) -> F!(y, z, pr, ps)

    # Alloc these locally for now, should really be done by caller to realise efficiency gains
    y = zeros(promote_type(S, T), l, q)
    gradF_flat = zeros(promote_type(S, T), l*q, n)
    ForwardDiff.jacobian!(gradF_flat, local_F!, y, x)
    for lp_q=1:q
        for lp_l=1:l
            out_gradF[lp_l, lp_q] = gradF_flat[l*(lp_q-1)+lp_l, :]      # Surpringly, I think this is the simplest way to do this
        end
    end
    return nothing
end

# Indexed gradF functions
function mm_fd_dn_gradF_i(dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q = dimspec
    len_idxs_F_i = length(idxs_F)
    local_F = (z::AbstractArray) -> mm_fd_dn_F_i(dimspec, F_i, z, pr, ps, idxs_F)
    gradF_flat = ForwardDiff.jacobian(local_F, x)
    gradF = [ gradF_flat[lp_idx, :] for lp_idx=1:len_idxs_F_i ]
    return gradF
end

function mm_fd_dn_gradF_i!(out_gradF_i::AbstractArray, dimspec::MPCCDimSpec, F_i!::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q = dimspec
    len_idxs_F_i = length(idxs_F)
    local_F! = (y::AbstractArray, z::AbstractArray) -> mm_fd_dn_F_i!(y, dimspec, F_i!, z, pr, ps, idxs_F)
    # Alloc these locally for now, should really be done by caller to realise efficiency gains
    y = zeros(promote_type(S, T), len_idxs_F_i)
    gradF_flat = zeros(promote_type(S, T), len_idxs_F_i, n)
    ForwardDiff.jacobian!(gradF_flat, local_F!, y, x)
    for lp_idx=1:len_idxs_F_i
        out_gradF_i[lp_idx] = gradF_flat[lp_idx, :]
    end
    return nothing
end

# Full gradFq functions
function mm_fd_dn_gradFq(dimspec::MPCCDimSpec, Fq::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    local_Fq = (z::AbstractArray) -> Fq(z, pr, ps)
    # return convert(Matrix{promote_type(S, T)}, ForwardDiff.jacobian(local_Fq, x))
    return ForwardDiff.jacobian(local_Fq, x)
end

function mm_fd_dn_gradFq!(out_gradFq::AbstractArray, dimspec::MPCCDimSpec, Fq!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, q = dimspec
    local_Fq! = (y::AbstractArray, z::AbstractArray) -> Fq!(y, z, pr, ps)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), q)
    ForwardDiff.jacobian!(out_gradFq, local_Fq!, y, x)
    return nothing
end

# Indexed gradFq functions
function mm_fd_dn_gradFq_i(dimspec::MPCCDimSpec, Fq_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    local_Fq = (z::AbstractArray) -> mm_fd_dn_Fq_i(dimspec, Fq_i, z, pr, ps, idxs_Fq)
    return ForwardDiff.jacobian(local_Fq, x)
end

function mm_fd_dn_gradFq_i!(out_gradFq_i::AbstractArray, dimspec::MPCCDimSpec, Fq_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n = dimspec
    len_idxs_Fq_i = length(idxs_Fq)
    local_Fq! = (y::AbstractArray, z::AbstractArray) -> mm_fd_dn_Fq!(y, dimspec, Fq_i!, z, pr, ps, idxs_Fq)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), len_idxs_Fq_i)

    # Unfortuantely, forewarddiff.jacobian returns the transpose of what we want.
    # So, currently, we're stuck allocating again: TODO sort this
    ForwardDiff.jacobian!(out_gradFq_i, local_Fq!, y, x)
    return nothing
end

# Hessian of f

function mm_fd_dn_hessf(dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    local_f(z) = f(z, pr, ps)[1]
    return ForwardDiff.hessian(local_f, x)
end

function mm_fd_dn_hessf!(out_hessf::AbstractArray, dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    local_f(z) = f(z, pr, ps)[1]
    ForwardDiff.hessian!(out_hessf, local_f, x)
    return nothing
end

# Full hessce functions

function mm_fd_dn_hessce(dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    # Call the indexed version...
    idxs_ce = collect(1:dimspec.me)
    return mm_fd_dn_hessce_i(dimspec, ce_i, x, pr, ps, idxs_ce::AbstractVector{Int64})
end

function mm_fd_dn_hessce!(out_hessce::AbstractArray, dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    # Call the indexed version...
    idxs_ce = collect(1:dimspec.me)
    mm_fd_dn_hessce_i!(out_hessce, dimspec, ce_i, x, pr, ps, idxs_ce)
    return nothing
end

# Indexed hessce functions

function mm_fd_dn_hessce_i(dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}

    @unpack n = dimspec

    len_idxs_ce_i = length(idxs_ce)
    local_ce_i = Vector{Function}(undef, len_idxs_ce_i)
    hessce_i = Vector{Matrix{T}}(undef, len_idxs_ce_i)

    for lp_ce in eachindex(idxs_ce)
        local_ce_i[lp_ce] = (z::AbstractArray) -> ce_i[idxs_ce[lp_ce]](z, pr, ps)[1]
        hessce_i[lp_ce] = ForwardDiff.hessian(local_ce_i[lp_ce], x)
    end
    return hessce_i
end

function mm_fd_dn_hessce_i!(out_hessce_i::AbstractArray, dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    len_idxs_ce_i = length(idxs_ce)
    local_ce_i = Vector{Function}(undef, len_idxs_ce_i)

    for lp_ce in eachindex(idxs_ce)
        local_ce_i[lp_ce] = (z::AbstractArray) -> ce_i[idxs_ce[lp_ce]](z, pr, ps)[1]
        ForwardDiff.hessian!(@view(out_hessce_i[lp_ce][:]), local_ce_i[lp_ce], x)
    end
    return nothing
end

# Full hessci functions

function mm_fd_dn_hessci(dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    # Call the indexed version...
    idxs_ci = collect(1:dimspec.mi)
    return mm_fd_dn_hessci_i(dimspec, ci_i, x, pr, ps, idxs_ci::AbstractVector{Int64})
end

function mm_fd_dn_hessci!(out_hessci::AbstractArray, dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    # Call the indexed version...
    idxs_ci = collect(1:dimspec.mi)
    mm_fd_dn_hessci_i!(out_hessci, dimspec, ci_i, x, pr, ps, idxs_ci)
    return nothing
end

# Indexed hessci functions

function mm_fd_dn_hessci_i(dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}

    @unpack n = dimspec

    len_idxs_ci_i = length(idxs_ci)
    local_ci_i = Vector{Function}(undef, len_idxs_ci_i)
    hessci_i = Vector{Matrix{T}}(undef, len_idxs_ci_i)

    for lp_ci in eachindex(idxs_ci)
        local_ci_i[lp_ci] = (z::AbstractArray) -> ci_i[idxs_ci[lp_ci]](z, pr, ps)[1]
        hessci_i[lp_ci] = ForwardDiff.hessian(local_ci_i[lp_ci], x)
    end
    return hessci_i
end

function mm_fd_dn_hessci_i!(out_hessci_i::AbstractArray, dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    len_idxs_ci_i = length(idxs_ci)
    local_ci_i = Vector{Function}(undef, len_idxs_ci_i)

    for lp_ci in eachindex(idxs_ci)
        local_ci_i[lp_ci] = (z::AbstractArray) -> ci_i[idxs_ci[lp_ci]](z, pr, ps)[1]
        ForwardDiff.hessian!(@view(out_hessci_i[lp_ci][:]), local_ci_i[lp_ci], x)
    end
    return nothing
end

# Full hessF functions
# # Matrix of dim l x q of n length vectors (the hessians wrt x)

function mm_fd_dn_hessF(dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}   # ::Matrix{Matrix{S}}

    @unpack n, l, q = dimspec

    local_F = Matrix{Function}(undef, l, q)       # TODO do we need to actually store all of them in an array?
    hessF = Matrix{Matrix{T}}(undef, l, q)

    for lp_q=1:q
        for lp_l=1:l
            local_F[lp_l, lp_q] = (z::AbstractArray) -> F_i[lp_l, lp_q](z, pr, ps)[1]
            hessF[lp_l, lp_q] = ForwardDiff.hessian(local_F[lp_l, lp_q], x)
        end
    end
    return hessF
end

function mm_fd_dn_hessF!(out_hessF::AbstractArray, dimspec::MPCCDimSpec,  F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q = dimspec

    local_F = Matrix{Function}(undef, l, q)

    for lp_q=1:q
        for lp_l=1:l
            local_F[lp_l, lp_q] = (z::AbstractArray) -> F_i[lp_l, lp_q](z, pr, ps)[1]
            ForwardDiff.hessian!(@view(out_hessF[lp_l, lp_q][:,:]), local_F[lp_l, lp_q], x)
        end
    end
    return nothing
end

# Indexed hessF functions
function mm_fd_dn_hessF_i(dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q = dimspec
    len_idxs_F_i = length(idxs_F)
    local_F_i = Vector{Function}(undef, len_idxs_F_i)
    hessF_i = Vector{Matrix{T}}(undef, len_idxs_F_i)
    for lp_F in eachindex(idxs_F)
        local_F_i[lp_F] = (z::AbstractArray) -> F_i[idxs_F[lp_F]](z, pr, ps)[1]
        hessF_i[lp_F] = ForwardDiff.hessian(local_F_i[lp_F], x)
    end
    return hessF_i
end

function mm_fd_dn_hessF_i!(out_hessF_i::AbstractArray, dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q = dimspec
    len_idxs_F_i = length(idxs_F)
    local_F_i = Vector{Function}(undef, len_idxs_F_i)       # TODO do we need to actually store all of them in an array?
    for lp_F in eachindex(idxs_F)
        local_F_i[lp_F] = (z::AbstractArray) -> F_i[idxs_F[lp_F]](z, pr, ps)[1]
        ForwardDiff.hessian!(@view(out_hessF_i[lp_F][:,:]), local_F_i[lp_F], x)
    end
    return nothing
end

# Full hessFq functions

function mm_fd_dn_hessFq(dimspec::MPCCDimSpec, Fq_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    # Call the indexed version...
    idxs_Fq = collect(1:dimspec.q)
    return mm_fd_dn_hessFq_i(dimspec, Fq_i, x, pr, ps, idxs_Fq::AbstractVector{Int64})
end

# Indexed hessFq functions

function mm_fd_dn_hessFq_i(dimspec::MPCCDimSpec, Fq_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_Fq::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    @unpack n = dimspec

    len_idxs_Fq_i = length(idxs_Fq)
    local_Fq_i = Vector{Function}(undef, len_idxs_Fq_i)
    hessFq_i = Vector{Matrix{T}}(undef, len_idxs_Fq_i)

    for lp_Fq in eachindex(idxs_Fq)
        local_Fq_i[lp_Fq] = (z::AbstractArray) -> Fq_i[idxs_Fq[lp_Fq]](z, pr, ps)[1]
        hessFq_i[lp_Fq] = ForwardDiff.hessian(local_Fq_i[lp_Fq], x)
    end
    return hessFq_i
end

# fdp functions

function mm_fd_dn_fdp(dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}   # ::AbstractVector{T}
    local_f(qr) = f(x, qr, ps)[1]
    return ForwardDiff.gradient(local_f, pr)
end

function mm_fd_dn_fdp!(out_fdp::AbstractArray, dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    local_f(z) = f(z, pr, ps)[1]
    ForwardDiff.gradient!(out_fdp, local_f, x)
    return nothing
end

# Full cedp functions
# # Vector of length r of vectors of length me
function mm_fd_dn_cedp(dimspec::MPCCDimSpec, ce::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    @unpack me, r = dimspec
    local_ce = (qr::AbstractArray) -> ce(x, qr, ps)
    temp_cedp = ForwardDiff.jacobian(local_ce, pr)
    cedp = [ [ temp_cedp[lp_me, lp_r] for lp_me in 1:me ] for lp_r in 1:r ]
    return cedp
end

function mm_fd_dn_cedp!(out_cedp::AbstractArray, dimspec::MPCCDimSpec, ce!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack me, r = dimspec
    local_ce! = (y::AbstractArray, qr::AbstractArray) -> ce!(y, x, qr, ps)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), me)
    temp_cedp = zeros(promote_type(S, T), me, r)
    ForwardDiff.jacobian!(temp_cedp, local_ce!, y, pr)
    for lp_r=1:r
        for lp_me=1:me
            out_cedp[lp_r][lp_me] = temp_cedp[lp_me, lp_r]
        end
    end
    return nothing
end

# Indexed cedp functions
function mm_fd_dn_cedp_i(dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    @unpack r = dimspec
    len_idxs_ce_i = length(idxs_ce)
    local_ce = (qr::AbstractArray) -> mm_fd_dn_ce_i(dimspec, ce_i, x, qr, ps, idxs_ce)
    temp_cedp = ForwardDiff.jacobian(local_ce, pr)
    cedp = [ [ temp_cedp[lp_me, lp_r] for lp_me in 1:len_idxs_ce_i ] for lp_r in 1:r ]
    return cedp
end

function mm_fd_dn_cedp_i!(out_cedp::AbstractArray, dimspec::MPCCDimSpec, ce_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack me, r = dimspec
    len_idxs_ce_i = length(idxs_ce)
    local_ce! = (y::AbstractArray, qr::AbstractArray) -> mm_fd_dn_ce_i!(y, dimspec, ce_i!, x, qr, ps, idxs_ce)
    # mm_fd_dn_ce_i!(out_ce_i::AbstractVector{R}, dimspec::MPCCDimSpec, ce_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {F <: Function, R <: Real, S <: Real, T <: Real}
    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage memory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), len_idxs_ce_i)
    temp_cedp = zeros(promote_type(S, T), len_idxs_ce_i, r)
    ForwardDiff.jacobian!(temp_cedp, local_ce!, y, pr)
    for lp_r=1:r
        for lp_me=1:len_idxs_ce_i
            out_cedp[lp_r][lp_me] = temp_cedp[lp_me, lp_r]
        end
    end
    return nothing
end

# Full cidp functions
# # Vector of length r of vectors of length mi
function mm_fd_dn_cidp(dimspec::MPCCDimSpec, ci::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    @unpack mi, r = dimspec
    local_ci = (qr::AbstractArray) -> ci(x, qr, ps)
    temp_cidp = ForwardDiff.jacobian(local_ci, pr)
    cidp = [ [ temp_cidp[lp_mi, lp_r] for lp_mi in 1:mi ] for lp_r in 1:r ]
    return cidp
end

function mm_fd_dn_cidp!(out_cidp::AbstractArray, dimspec::MPCCDimSpec, ci!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack mi, r = dimspec
    local_ci! = (y::AbstractArray, qr::AbstractArray) -> ci!(y, x, qr, ps)

    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage mimory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), mi)
    temp_cidp = zeros(promote_type(S, T), mi, r)
    ForwardDiff.jacobian!(temp_cidp, local_ci!, y, pr)
    for lp_r=1:r
        for lp_mi=1:mi
            out_cidp[lp_r][lp_mi] = temp_cidp[lp_mi, lp_r]
        end
    end
    return nothing
end

# Indexed cidp functions
function mm_fd_dn_cidp_i(dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    @unpack r = dimspec
    len_idxs_ci_i = length(idxs_ci)
    local_ci = (qr::AbstractArray) -> mm_fd_dn_ci_i(dimspec, ci_i, x, qr, ps, idxs_ci)
    temp_cidp = ForwardDiff.jacobian(local_ci, pr)
    cidp = [ [ temp_cidp[lp_mi, lp_r] for lp_mi in 1:len_idxs_ci_i ] for lp_r in 1:r ]
    return cidp
end

function mm_fd_dn_cidp_i!(out_cidp::AbstractArray, dimspec::MPCCDimSpec, ci_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack mi, r = dimspec
    len_idxs_ci_i = length(idxs_ci)
    local_ci! = (y::AbstractArray, qr::AbstractArray) -> mm_fd_dn_ci_i!(y, dimspec, ci_i!, x, qr, ps, idxs_ci)
    # mm_fd_dn_ci_i!(out_ci_i::AbstractVector{R}, dimspec::MPCCDimSpec, ci_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {F <: Function, R <: Real, S <: Real, T <: Real}
    # ForwardDiff asks for a temporary storage area, presumably so that caller can manage mimory allocs
    # We do this here for now, perhaps move it upstream later.
    y = zeros(promote_type(S, T), len_idxs_ci_i)
    temp_cidp = zeros(promote_type(S, T), len_idxs_ci_i, r)
    ForwardDiff.jacobian!(temp_cidp, local_ci!, y, pr)
    for lp_r=1:r
        for lp_mi=1:len_idxs_ci_i
            out_cidp[lp_r][lp_mi] = temp_cidp[lp_mi, lp_r]
        end
    end
    return nothing
end

# Full Fdp functions

function mm_fd_dn_Fdp(dimspec::MPCCDimSpec, F::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    @unpack n, l, q, r = dimspec

    # ForwardDiff.jacobian can work with multidimensional arrays, so we do that
    local_F = (qr::AbstractArray) -> F(x, qr, ps)
    Fdp = Vector{Matrix{T}}(undef, r)
    Fdp_flat = ForwardDiff.jacobian(local_F, pr)
    for lp_pr=1:r
        Fdp[lp_pr] = reshape(Fdp_flat[:, lp_pr], (l, q))
    end
    return Fdp
end

function mm_fd_dn_Fdp!(out_Fdp::AbstractArray, dimspec::MPCCDimSpec, F!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, l, q, r = dimspec

    # ForwardDiff.jacobian can work with multidimensional arrays, so we do that
    local_F! = (z::AbstractArray, qr::AbstractArray) -> F!(z, x, qr, ps)

    # Alloc these locally for now, should really be done by caller to realise efficiency gains
    y = zeros(promote_type(S, T), l, q)
    Fdp_flat = zeros(promote_type(S, T), l*q, r)
    ForwardDiff.jacobian!(Fdp_flat, local_F!, y, pr)
    for lp_pr=1:r
        out_Fdp[lp_pr] = reshape(Fdp_flat[:, lp_pr], (l, q))
    end
    return nothing
end

# Indexed Fdp
function mm_fd_dn_Fdp_i(dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q, r = dimspec

    len_idxs_F_i = length(idxs_F)
    local_F = (qr::AbstractArray) -> mm_fd_dn_F_i(dimspec, F_i, x, qr, ps, idxs_F)
    Fdp = Vector{Vector{T}}(undef, r)
    Fdp_flat = ForwardDiff.jacobian(local_F, pr)
    for lp_pr=1:r
        Fdp[lp_pr] = Fdp_flat[(1:len_idxs_F_i)*lp_pr]
    end
    return Fdp
end

function mm_fd_dn_Fdp_i!(out_Fdp::AbstractArray, dimspec::MPCCDimSpec, F_i!::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q, r = dimspec

    len_idxs_F_i = length(idxs_F)
    local_F! = (z::AbstractArray, qr::AbstractArray) -> mm_fd_dn_F_i!(z, dimspec, F_i!, x, qr, ps, idxs_F)

    # Alloc these locally for now, should really be done by caller to realise efficiency gains
    y = zeros(promote_type(S, T), len_idxs_F_i)
    Fdp_flat = zeros(promote_type(S, T), len_idxs_F_i * r)
    ForwardDiff.jacobian!(Fdp_flat, local_F!, y, pr)
    for lp_pr=1:r
        for lp_F=1:len_idxs_F_i
            out_Fdp[lp_pr][lp_F] = Fdp_flat[lp_F+((lp_pr-1)*len_idxs_F_i)]
        end
    end
    return nothing
end

# Full gradfdp functions

function mm_fd_dn_gradfdp(dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    @unpack n, r = dimspec
    local_gradf = (qr::AbstractArray) -> mm_fd_dn_gradf(dimspec, f, x, qr, ps)
    gradfdp_flat = ForwardDiff.jacobian(local_gradf, pr)
    gradfdp = [ gradfdp_flat[:, lp_r] for lp_r in 1:r ]
    return gradfdp
end

function mm_fd_dn_gradfdp!(out_gradfdp::AbstractArray, dimspec::MPCCDimSpec, f::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, r = dimspec
    local_gradf! = (z::AbstractArray, qr::AbstractArray) -> mm_fd_dn_gradf!(z, dimspec, f, x, qr, ps)
    y = zeros(promote_type(S, T), n)
    gradfdp_flat = zeros(promote_type(S, T), n, r)
    ForwardDiff.jacobian!(gradfdp_flat, local_gradf!, y, pr)
    for lp_pr=1:r
        out_gradfdp[lp_pr] = gradfdp_flat[:, lp_pr]
    end
    return nothing
end

# Full gradcedp functions
function mm_fd_dn_gradcedp(dimspec::MPCCDimSpec, ce::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    @unpack n, me, r = dimspec
    local_gradce = (qr::AbstractArray) -> mm_fd_dn_gradce(dimspec, ce, x, qr, ps)
    gradcedp_flat = ForwardDiff.jacobian(local_gradce, pr)
    gradcedp = [ reshape(gradcedp_flat[:, lp_pr], (me, n)) for lp_pr in 1:r ]
    return gradcedp
end

function mm_fd_dn_gradcedp!(out_gradcedp::AbstractArray, dimspec::MPCCDimSpec, ce!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, me, r = dimspec
    local_gradce! = (z::AbstractArray, qr::AbstractArray) -> mm_fd_dn_gradce!(z, dimspec, ce!, x, qr, ps)

    # Local alloc, should perhaps do in caller
    y = zeros(promote_type(S, T), me, n)
    gradcedp_flat = zeros(promote_type(S, T), n*me, r)
    ForwardDiff.jacobian!(gradcedp_flat, local_gradce!, y, pr)
    for lp_pr=1:r
        out_gradcedp[lp_pr] = reshape(gradcedp_flat[:, lp_pr], (me, n))
    end
    return nothing
end

# Indexed gradcedp
function mm_fd_dn_gradcedp_i(dimspec::MPCCDimSpec, ce_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    @unpack n, me, r = dimspec
    len_idxs_ce_i = length(idxs_ce)
    local_gradce = (qr::AbstractArray) -> mm_fd_dn_gradce_i(dimspec, ce_i, x, qr, ps, idxs_ce)
    gradcedp_flat = ForwardDiff.jacobian(local_gradce, pr)
    gradcedp = [ reshape(gradcedp_flat[:, lp_pr], (len_idxs_ce_i, n)) for lp_pr in 1:r ]
    return gradcedp
end

function mm_fd_dn_gradcedp_i!(out_gradcedp::AbstractArray, dimspec::MPCCDimSpec, ce_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ce::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n, me, r = dimspec
    len_idxs_ce_i = length(idxs_ce)
    local_gradce! = (z::AbstractArray, qr::AbstractArray) -> mm_fd_dn_gradce_i!(z, dimspec, ce_i!, x, qr, ps, idxs_ce)

    # Local alloc, should perhaps do in caller
    y = zeros(promote_type(S, T), len_idxs_ce_i, n)
    gradcedp_flat = zeros(promote_type(S, T), n*me, r)
    ForwardDiff.jacobian!(gradcedp_flat, local_gradce!, y, pr)
    for lp_pr=1:r
        out_gradcedp[lp_pr] = reshape(gradcedp_flat[:, lp_pr], (len_idxs_ce_i, n))
    end
    return nothing
end

# Full gradcidp functions
function mm_fd_dn_gradcidp(dimspec::MPCCDimSpec, ci::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {S <: Real, T <: Real}
    @unpack n, mi, r = dimspec
    local_gradci = (qr::AbstractArray) -> mm_fd_dn_gradci(dimspec, ci, x, qr, ps)
    gradcidp_flat = ForwardDiff.jacobian(local_gradci, pr)
    gradcidp = [ reshape(gradcidp_flat[:, lp_pr], (mi, n)) for lp_pr in 1:r ]
    return gradcidp
end

function mm_fd_dn_gradcidp!(out_gradcidp::AbstractArray, dimspec::MPCCDimSpec, ci!::Function, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64})::Nothing where {S <: Real, T <: Real}
    @unpack n, mi, r = dimspec
    local_gradci! = (z::AbstractArray, qr::AbstractArray) -> mm_fd_dn_gradci!(z, dimspec, ci!, x, qr, ps)

    # Local alloc, should perhaps do in caller
    y = zeros(promote_type(S, T), mi, n)
    gradcidp_flat = zeros(promote_type(S, T), n*mi, r)
    ForwardDiff.jacobian!(gradcidp_flat, local_gradci!, y, pr)
    for lp_pr=1:r
        out_gradcidp[lp_pr] = reshape(gradcidp_flat[:, lp_pr], (mi, n))
    end
    return nothing
end

# Indexed gradcidp
function mm_fd_dn_gradcidp_i(dimspec::MPCCDimSpec, ci_i::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    @unpack n, mi, r = dimspec
    len_idxs_ci_i = length(idxs_ci)
    local_gradci = (qr::AbstractArray) -> mm_fd_dn_gradci_i(dimspec, ci_i, x, qr, ps, idxs_ci)
    gradcidp_flat = ForwardDiff.jacobian(local_gradci, pr)
    gradcidp = [ reshape(gradcidp_flat[:, lp_pr], (len_idxs_ci_i, n)) for lp_pr in 1:r ]
    return gradcidp
end

function mm_fd_dn_gradcidp_i!(out_gradcidp::AbstractArray, dimspec::MPCCDimSpec, ci_i!::AbstractVector{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_ci::AbstractVector{Int64})::Nothing where {F <: Function, S <: Real, T <: Real}
    @unpack n, mi, r = dimspec
    len_idxs_ci_i = length(idxs_ci)
    local_gradci! = (z::AbstractArray, qr::AbstractArray) -> mm_fd_dn_gradci_i!(z, dimspec, ci_i!, x, qr, ps, idxs_ci)

    # Local alloc, should perhaps do in caller
    y = zeros(promote_type(S, T), len_idxs_ci_i, n)
    gradcidp_flat = zeros(promote_type(S, T), n*mi, r)
    ForwardDiff.jacobian!(gradcidp_flat, local_gradci!, y, pr)
    for lp_pr=1:r
        out_gradcidp[lp_pr] = reshape(gradcidp_flat[:, lp_pr], (len_idxs_ci_i, n))
    end
    return nothing
end

# Full gradFdp functions
function mm_fd_dn_gradFdp(dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}   # ::Vector{Matrix{Vector{S}}}
    @unpack n, l, q, r = dimspec
    local_gradF_i = Matrix{Function}(undef, l, q)
    gradFdp = [ [ zeros(promote_type(S, T), n) for lp_l in 1:l, lp_q in 1:q ] for lp_r in 1:r ]
    for lp_q=1:q
        for lp_l=1:l
            local_gradF_i[lp_l, lp_q] = (qr::AbstractArray) -> mm_fd_dn_gradF_i(dimspec, F_i, x, qr, ps, [CartesianIndex(lp_l, lp_q)])[1]
            temp_graddp_flat_i = ForwardDiff.jacobian(local_gradF_i[lp_l, lp_q], pr)
            for lp_r=1:r
                gradFdp[lp_r][lp_l, lp_q][:] = temp_graddp_flat_i[:, lp_r]
            end
        end
    end
    return gradFdp
end

function mm_fd_dn_gradFdp!(out_gradFdp::AbstractArray, dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}) where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q, r = dimspec
    local_gradF_i = Matrix{Function}(undef, l, q)
    # Local alloc, should perhaps do in caller
    y = zeros(promote_type(S, T), n)
    gradFdp_flat_i = zeros(promote_type(S, T), n, r)
    for lp_q=1:q
        for lp_l=1:l
            # Can't get it working fully mutating because it wants to mutate into a vector of vectors,
            # and forwarddiff needs just a plain vector.  So doesn't work.

            # local_gradF_i![lp_l, lp_q] = function (z::AbstractArray, qr::AbstractArray)
            #     vec_of_z = [ z ]
            #     mm_fd_dn_gradF_i!(vec_of_z, dimspec, F_i!, x, qr, ps, [CartesianIndex(lp_l, lp_q)])
            #     println("vec of z: ", vec_of_z)
            # end
            # ForwardDiff.jacobian!(gradFdp_flat_i, local_gradF_i![lp_l, lp_q], y, pr)
            local_gradF_i[lp_l, lp_q] = (qr::AbstractArray) -> mm_fd_dn_gradF_i(dimspec, F_i, x, qr, ps, [CartesianIndex(lp_l, lp_q)])[1]
            ForwardDiff.jacobian!(gradFdp_flat_i, local_gradF_i[lp_l, lp_q], pr)
            for lp_r=1:r
                out_gradFdp[lp_r][lp_l, lp_q][:] = gradFdp_flat_i[:, lp_r]
            end
        end
    end
    return nothing
end

# Indexed gradFdp functions
function mm_fd_dn_gradFdp_i(dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {F <: Function, S <: Real, T <: Real}   # ::Vector{Matrix{Vector{S}}}
    @unpack n, l, q, r = dimspec
    len_idxs_F_i = length(idxs_F)
    local_gradF_i = Vector{Function}(undef, len_idxs_F_i)

    # This time we need to do some construction work...
    gradFdp = [ [ zeros(promote_type(S, T), n) for lp_idx in 1:len_idxs_F_i ] for lp_r in 1:r ]

    # Can probably do this in one go rather than loop TODO
    for lp_idx=1:len_idxs_F_i
        local_gradF_i[lp_idx] = (qr::AbstractArray) -> mm_fd_dn_gradF_i(dimspec, F_i, x, qr, ps, [ idxs_F[lp_idx] ])[1]
        temp_graddp_flat_i = ForwardDiff.jacobian(local_gradF_i[lp_idx], pr)
        for lp_r=1:r
            gradFdp[lp_r][lp_idx][:] = temp_graddp_flat_i[:, lp_r]
        end
    end
    return gradFdp
end

function mm_fd_dn_gradFdp_i!(out_gradFdp::AbstractArray, dimspec::MPCCDimSpec, F_i::AbstractMatrix{F}, x::AbstractVector{S}, pr::AbstractVector{T}, ps::AbstractVector{Int64}, idxs_F::Vector{CartesianIndex{2}}) where {F <: Function, S <: Real, T <: Real}
    @unpack n, l, q, r = dimspec
    len_idxs_F_i = length(idxs_F)
    local_gradF_i = Vector{Function}(undef, len_idxs_F_i)
    gradFdp_flat_i = zeros(promote_type(S, T), n, r)

    # Can probably do this in one go rather than loop TODO
    for lp_idx=1:len_idxs_F_i
        local_gradF_i[lp_idx] = (qr::AbstractArray) -> mm_fd_dn_gradF_i(dimspec, F_i, x, qr, ps, [ idxs_F[lp_idx] ])[1]
        ForwardDiff.jacobian!(gradFdp_flat_i, local_gradF_i[lp_idx], pr)
        for lp_r=1:r
            out_gradFdp[lp_r][lp_idx][:] = gradFdp_flat_i[:, lp_r]
        end

    end
    return nothing
end


text:
import FactorGraph
using FactorGraph.ContractionOrder
using FactorGraph.ContractionOrder: analyze_contraction, contract_pair!, evaluate_costs
using OMEinsum

using Test
@testset "analyze contraction" begin
    incidence_list = IncidenceList(Dict('A' => ['a', 'b', 'k', 'o', 'f'], 'B'=>['a', 'c', 'd', 'm', 'f'], 'C'=>['b', 'c', 'e', 'f'], 'D'=>['e'], 'E'=>['d', 'f']), openedges=['c', 'f', 'o'])
    info = analyze_contraction(incidence_list, 'A', 'B')
    @test Set(info.l1) == Set(['k'])
    @test Set(info.l2) == Set(['m'])
    @test Set(info.l12) == Set(['a'])
    @test Set(info.l01) == Set(['b','o'])
    @test Set(info.l02) == Set(['c', 'd'])
    @test Set(info.l012) == Set(['f'])
end

@testset "tree greedy" begin
    incidence_list = IncidenceList(Dict('A' => ['a', 'b'], 'B'=>['a', 'c', 'd'], 'C'=>['b', 'c', 'e', 'f'], 'D'=>['e'], 'E'=>['d', 'f']))
    log2_edge_sizes = Dict([c=>i for (i,c) in enumerate(['a', 'b', 'c', 'd', 'e', 'f'])]...)
    il = copy(incidence_list)
    contract_pair!(il, 'A', 'B', log2_edge_sizes)
    target = IncidenceList(Dict('A' => ['b', 'c', 'd'], 'C'=>['b', 'c', 'e', 'f'], 'D'=>['e'], 'E'=>['d', 'f']))
    @test il.v2e == target.v2e
    @test length(target.e2v) == length(il.e2v)
    for (k,v) in il.e2v
        @test sort(target.e2v[k]) == sort(v)
    end
    costs = evaluate_costs(MinSpaceOut(), incidence_list, log2_edge_sizes)
    @test costs == Dict(('A', 'B')=>9-0.01, ('A', 'C')=>15-0.02, ('B','C')=>18-0.03, ('B','E')=>10-0.04, ('C','D')=>11-0.05, ('C', 'E')=>14-0.06)
    tree, log2_tcs, log2_scs = tree_greedy(incidence_list, log2_edge_sizes)
    @test log2_tcs == [10.0, 16, 15, 9]
    @test tree == ContractionTree(ContractionTree('A', 'B'), ContractionTree(ContractionTree('C', 'D'), 'E'))
    @test timespace_complexity(incidence_list, tree, log2_edge_sizes) == (log2(exp2(10)+exp2(16)+exp2(15)+exp2(9)), 11)
    vertices = ['A', 'B', 'C', 'D', 'E']
    optcode1 = FactorGraph.parse_eincode(incidence_list, tree, vertices=vertices)
    @test optcode1 isa OMEinsum.NestedEinsum
    tree2 = FactorGraph.parse_tree(optcode1, vertices)
    @test tree2 == tree

    eincode = ein"ab,acd,bcef,e,df->"
    size_dict = Dict([c=>(1<<i) for (i,c) in enumerate(['a', 'b', 'c', 'd', 'e', 'f'])]...)
    optcode2 = FactorGraph.optimizecode_greedy(eincode, size_dict)
    tc, sc = timespace_complexity(optcode2, log2_edge_sizes)
    @test tc ≈ log2(exp2(10)+exp2(16)+exp2(15)+exp2(9))
    @test sc == 11
    @test optcode1 == optcode2
end

ext: jl


text:
using Documenter, GPLikelihoods

makedocs(;
    sitename="GPLikelihoods.jl",
    format=Documenter.HTML(),
    modules=[GPLikelihoods],
    pages=["Home" => "index.md", "API" => "api.md"],
    repo="https://github.com/JuliaGaussianProcesses/GPLikelihoods.jl/blob/{commit}{path}#L{line}",
    authors="JuliaGaussianProcesses organization",
    strict=true,
    checkdocs=:exports,
    #doctestfilters=JuliaGPsDocs.DOCTEST_FILTERS,
)

deploydocs(; repo="github.com/JuliaGaussianProcesses/GPLikelihoods.jl", push_preview=true)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CGAL2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CGAL2")
JLLWrappers.@generate_main_file("CGAL2", UUID("a133c068-ba04-5466-9207-ec1c2ac43820"))
end  # module CGAL2_jll


text:
module Visualize_Beam

import Modia3D

mat1  = Modia3D.VisualMaterial(color="LightBlue", transparency=0.5)

world = Modia3D.Object3D(
    feature = Modia3D.Visual(shape=Modia3D.CoordinateSystem(length=1.0)))
beam  = Modia3D.Object3D(parent=world, translation=[1.0, 1.0, 1.0],
    feature = Modia3D.Visual(shape=Modia3D.Beam(axis=1, length=0.4, width=0.5, thickness=0.3), visualMaterial=mat1))

Modia3D.visualizeWorld!(world,
    sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.7, enableContactDetection=false))

println("... test/old/Visualize_Beam.jl completed.")

end


text:
using Test

# import CLIMAParameters
import CLIMAParameters
const CP = CLIMAParameters

# read parameters needed for tests
full_parameter_set = CP.create_toml_dict(Float64; dict_type = "alias")

const CPP = CP.Planet
struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const param_set_cpp = EarthParameterSet()

universal_constant_aliases = [
    "gas_constant",
    "light_speed",
    "h_Planck",
    "k_Boltzmann",
    "Stefan",
    "astro_unit",
    "avogad",
]

# CP modules list:
module_names = [
    CP,
    CP.Planet,
    CP.SubgridScale,
    # CP.Atmos.EDMF,
    # CP.Atmos.SubgridScale,
    CP.Atmos.Microphysics_0M,
    CP.Atmos.Microphysics,
    CP.SurfaceFluxes.UniversalFunctions,
]

CP_parameters = Dict(mod => String.(names(mod)) for mod in module_names)
logfilepath1 = joinpath(@__DIR__, "toml", "log_file_test_1.toml")

@testset "parameter file interface tests" begin

    @testset "load with name or alias" begin
        @test_throws AssertionError CP.create_toml_dict(
            Float64;
            dict_type = "not name or alias",
        )

    end

    @testset "TOML - CliMAParameters.jl consistency" begin
        # tests to check parameter consistency of new toml files with existing
        # CP defaults.

        k_found = [0]
        for (k, v) in full_parameter_set #iterates over data (by alias)

            for mod in module_names
                k_pair = CP.get_parameter_values(full_parameter_set, k)
                k_value = last(k_pair)
                if k in CP_parameters[mod]
                    k_found[1] = 1
                    cp_k = getfield(mod, Symbol(k))
                    if !(k in universal_constant_aliases)
                        @test (k_value ≈ cp_k(param_set_cpp))
                    else #universal parameters have no argument
                        @test (k_value ≈ cp_k())
                    end
                    #for the logfile test later:
                    CP.get_parameter_values!(
                        full_parameter_set,
                        k,
                        string(nameof(mod)),
                    )
                end
            end
            if k_found[1] == 0
                println("on trying alias: ", k)
                @warn("did not find in any modules")
            end

            k_found[1] = 0
        end

        #create a dummy log file listing where CLIMAParameter lives
        CP.write_log_file(full_parameter_set, logfilepath1)
    end

    @testset "Parameter logging" begin

        #read in log file as new parameter file and rerun test.
        full_parameter_set_from_log = CP.create_toml_dict(
            Float64;
            override_file = logfilepath1,
            dict_type = "alias",
        )
        k_found = [0]
        for (k, v) in full_parameter_set_from_log #iterates over data (by alias)
            for mod in module_names
                k_pair = CP.get_parameter_values(full_parameter_set_from_log, k)
                k_value = last(k_pair)
                if k in CP_parameters[mod]
                    k_found[1] = 1
                    cp_k = getfield(mod, Symbol(k))
                    if !(k in universal_constant_aliases)
                        @test (k_value ≈ cp_k(param_set_cpp))
                    else #universal parameters have no argument
                        @test (k_value ≈ cp_k())
                    end
                end
            end
            if k_found[1] == 0
                println("on trying alias: ", k)
                @warn("did not find in any modules")
            end
            k_found[1] = 0
        end
    end

    @testset "parameter arrays" begin
        # Tests to check if file parsing, extracting and logging of parameter
        # values also works with array-valued parameters

        # Create parameter dict consisting of the parameters contained in the
        # default parameter file ("parameters.toml") and additional (array valued)
        # parameters ("array_parameters.toml").
        path_to_array_params =
            joinpath(@__DIR__, "toml", "array_parameters.toml")
        # parameter struct of type Float64 (default)
        toml_dict = CP.create_toml_dict(
            Float64;
            override_file = path_to_array_params,
            dict_type = "name",
        )
        # parameter struct of type Float32
        toml_dict_f32 = CP.create_toml_dict(
            Float32;
            override_file = path_to_array_params,
            dict_type = "name",
        )

        # true parameter values (used to check if values are correctly read from
        # the toml file)
        true_param_1 = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
        true_param_2 = [0.0, 1.0, 1.0, 2.0, 3.0, 5.0, 8.0, 13.0, 21.0, 34.0]
        true_param_3 = [9.81, 10.0]
        true_param_4 = 299792458
        true_params = [true_param_1, true_param_2, true_param_3, true_param_4]
        param_names = [
            "array_parameter_1",
            "array_parameter_2",
            "gravitational_acceleration",
            "light_speed",
        ]

        # Let's assume that parameter_vector_1 and parameter_vector_2 are used
        # in a module called "Test"
        mod = "Test"

        # Get parameter values and add information on the module where the
        # parameters are used.
        for i in range(1, stop = length(true_params))

            param_pair =
                CP.get_parameter_values!(toml_dict, param_names[i], mod)
            param = last(param_pair)
            @test param == true_params[i]
            # Check if the parameter is of the correct type. It should have
            # the same type as the ParamDict, which is specified by the
            # `float_type` argument to `create_toml_dict`.
            @test eltype(param) == Float64

            param_f32_pair =
                CP.get_parameter_values!(toml_dict_f32, param_names[i], mod)
            param_f32 = last(param_f32_pair)
            @test eltype(param_f32) == Float32

        end

        # Get several parameter values (scalar and arrays) at once
        params = CP.get_parameter_values(toml_dict, param_names)
        for j in 1:length(param_names)
            param_val = last(params[j])
            @test param_val == true_params[j]
        end

        # Write parameters to log file
        mktempdir(@__DIR__) do path
            logfilepath2 = joinpath(path, "log_file_test_2.toml")
            CP.write_log_file(toml_dict, logfilepath2)
        end

        # `toml_dict` and `full_param_set` contain different values for the
        # `gravitational_acceleration` parameter. The merged parameter set should
        # contain the value from `toml_dict`.
        full_param_set = CP.create_toml_dict(Float64; dict_type = "name")
        merged_param_set =
            CP.merge_override_default_values(toml_dict, full_param_set)
        grav_pair = CP.get_parameter_values(
            merged_param_set,
            "gravitational_acceleration",
        )
        grav = last(grav_pair)
        @test grav == true_param_3
    end

    @testset "checks for overrides" begin
        full_param_set = CP.create_toml_dict(
            Float64;
            override_file = joinpath(@__DIR__, "toml", "override_typos.toml"),
            dict_type = "name",
        )
        mod = "test_module_name"
        CP.get_parameter_values!(full_param_set, "light_speed", mod)

        mktempdir(@__DIR__) do path
            logfilepath3 = joinpath(path, "log_file_test_3.toml")
            @test_logs (:warn,) CP.log_parameter_information(
                full_param_set,
                logfilepath3,
            )
            @test_throws ErrorException CP.log_parameter_information(
                full_param_set,
                logfilepath3,
                strict = true,
            )
        end
    end

end

rm(logfilepath1; force = true)

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule H2AgModel_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("H2AgModel")
JLLWrappers.@generate_main_file("H2AgModel", UUID("a822e7dc-5958-54de-9eb2-fa386864eb49"))
end  # module H2AgModel_jll


text:
using AnimatedPlots

sin_graph = StaticGraph(sin, color=SFML.red)
sin_graph.accuracy = 5
cos_graph = StaticGraph(cos, color=SFML.blue)
sin_graph.accuracy = 5
tan_graph = StaticGraph(x -> tan(0.1x), color=SFML.green)
sin_graph.accuracy = 5

plot(sin_graph)
plot(cos_graph)
plot(tan_graph)

waitfor(current_window())


text:
pages = [
    "Home" => "index.md",
    "Tutorials" => Any[
        "tutorials/numerical_integrals.md",
        "tutorials/differentiating_integrals.md"
    ],
    "Basics" => Any[
        "basics/IntegralProblem.md",
        "basics/solve.md",
        "basics/FAQ.md"
    ],
    "Solvers" => Any[
        "solvers/IntegralSolvers.md"
    ]
]


text:
module AeroMDAO

## Libraries
#==========================================================================================#

using LinearAlgebra
using StaticArrays
using CoordinateTransformations
using Rotations
using ForwardDiff: jacobian!
using DiffResults: JacobianResult, jacobian, value
using PrettyTables
using StructArrays

using Statistics: mean

using SplitApplyCombine: combinedimsview, combinedims
export combinedimsview, combinedims

using ComponentArrays

export ComponentVector, ComponentArray, valkeys

using Setfield
using LabelledArrays

## Methods to be extended in submodules
#==========================================================================================#

function collocation_point end

function velocity end

function solve_linear end

function solve_nonlinear end

function solve_linear! end

function solve_nonlinear! end

function solve_system end

function properties end

function surface_velocities end

function surface_coefficients end

export solve_linear, solve_nonlinear, solve_nonlinear!, collocation_point

## Math tools
#==========================================================================================#

include("Tools/MathTools.jl")
import .MathTools: forward_sum, forward_difference, forward_division, weighted_vector, vectarray, slope, splitat, adj3, columns, extend_yz, reflect_mapper, cosine_interp, structtolist, inverse_rotation, rotation, affine_2D, Point2D, Point3D, x, y, z, reshape_array, midpair_map, partition, uniform_spacing, linear_spacing, cosine_spacing, sine_spacing

# export forward_sum, forward_difference, forward_division, weighted_vector, vectarray, slope, splitat, adj3, columns, extend_yz, reflect_mapper, cosine_interp, structtolist, inverse_rotation, rotation, affine_2D, Point2D, Point3D, reshape_array, midpair_map, partition, uniform_spacing, linear_spacing, cosine_spacing, sine_spacing

## Non-dimensionalization
#==========================================================================================#

include("Tools/NonDimensional.jl")
import .NonDimensional: dynamic_pressure, force_coefficient, moment_coefficient, rate_coefficient, pressure_coefficient, aerodynamic_coefficients, reynolds_number, force, moment

export dynamic_pressure, force_coefficient, moment_coefficient, rate_coefficient, pressure_coefficient, reynolds_number, force, moment

## Panels
#==========================================================================================#

include("Geometry/PanelGeometry/PanelGeometry.jl")
import .PanelGeometry: AbstractPanel, AbstractPanel2D, Panel2D, WakePanel2D, AbstractPanel3D, Panel3D, p1, p2, p3, p4, zero, panel_length, transform_panel, transform_panel_points, panel_angle, panel_tangent, panel_normal, panel_location, panel_area, panel_coordinates, transform, midpoint, panel_points, wake_panel, wake_panels, reverse_panel, panel_velocity, panel_scalar, trailing_edge_panel, get_surface_values, panel_vector, distance, average_chord, average_width, wetted_area, make_panels, local_coordinate_system

export AbstractPanel, AbstractPanel2D, Panel2D, WakePanel2D, AbstractPanel3D, Panel3D, p1, p2, p3, p4, transform, panel_normal, midpoint, panel_location, panel_tangent, panel_points, distance, wake_panel, wake_panels, panel_area, reverse_panel, panel_length, transform_panel, panel_angle, panel_vector, panel_velocity, panel_scalar, trailing_edge_panel, get_surface_values, average_chord, average_width, wetted_area, make_panels, local_coordinate_system

## Aircraft geometry
#==========================================================================================#

include("Geometry/AircraftGeometry/AircraftGeometry.jl")
import .AircraftGeometry: AbstractAircraft, AbstractWing, AbstractFoil, Foil, arc_length, kulfan_CST, naca4, camber_CST, make_panels, read_foil, leading_edge_index, upper_surface, lower_surface, split_surface, coordinates_to_camber_thickness, camber_thickness_to_coordinates, camber_thickness, camber_thickness_to_coordinates, cosine_spacing, camber_thickness_to_CST, coordinates_to_CST, maximum_thickness_to_chord, translate, rotate, affine, scale, camber, camber_line, control_surface, # 2D setups
Fuselage, projected_area, length, cosine_spacing, # Fuselage
HalfWing, HalfWingSection, Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, aspect_ratio, projected_area, taper_ratio, leading_edge, trailing_edge, chop_leading_edge, chop_trailing_edge, chop_wing, chop_sections, chop_coordinates, chop_spanwise_sections, chop_chords, chop_spans, wing_bounds, make_panels, mesh_chords, mesh_wing, mesh_cambers, max_thickness_to_chord_ratio_sweeps, mean_aerodynamic_center, panel_wing, number_of_spanwise_panels, symmetric_spacing, coordinates, chord_coordinates, camber_coordinates, surface_coordinates, foils, chords, twists, spans, dihedrals, sweeps, position, orientation, WingMesh, chord_panels, camber_panels, normal_vectors, surface_panels, AbstractSpacing, Sine, Cosine, Uniform, properties

export AbstractAircraft, AbstractWing, AbstractFoil, Foil, arc_length, kulfan_CST, naca4, camber_CST, make_panels, read_foil, leading_edge_index, upper_surface, lower_surface, split_surface, coordinates_to_camber_thickness, camber_thickness_to_coordinates, camber_thickness, camber_thickness_to_coordinates, cosine_spacing, camber_thickness_to_CST, coordinates_to_CST, maximum_thickness_to_chord, translate, rotate, affine, scale, camber, camber_line, control_surface, # 2D setups
Fuselage, projected_area, length, cosine_spacing, # Fuselage
HalfWing, HalfWingSection, Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, aspect_ratio, projected_area, taper_ratio, leading_edge, trailing_edge, chop_leading_edge, chop_trailing_edge, chop_wing, chop_sections, chop_coordinates, chop_spanwise_sections, chop_chords, chop_spans, wing_bounds, make_panels, mesh_chords, mesh_wing, mesh_cambers, max_thickness_to_chord_ratio_sweeps, mean_aerodynamic_center, panel_wing, number_of_spanwise_panels, symmetric_spacing, coordinates, chord_coordinates, camber_coordinates, surface_coordinates, foils, chords, twists, spans, dihedrals, sweeps, position, orientation, WingMesh, chord_panels, camber_panels, normal_vectors, surface_panels, AbstractSpacing, Sine, Cosine, Uniform, properties

make_horseshoes(wing :: WingMesh) = StructArray(Horseshoe.(chord_panels(wing), normal_vectors(wing)))

make_vortex_rings(wing :: WingMesh) = StructArray(VortexRing.(camber_panels(wing)))

export make_horseshoes, make_vortex_rings

## Laplace
#==========================================================================================#

include("Tools/Laplace.jl")
import .Laplace: Uniform2D, potential, stream,  cartesian_to_freestream, freestream_to_cartesian

export Uniform2D, stream, vortex_stream_1, vortex_stream_2, source_stream, cartesian_to_freestream, freestream_to_cartesian

## Aerodynamic analyses
#==========================================================================================#

## Doublet-source panel method

include("Aerodynamics/DoubletSource/DoubletSource.jl")
import .DoubletSource: doublet_matrix, source_matrix, boundary_vector, wake_panels, source_strengths, surface_velocities, lift_coefficient

export doublet_matrix, source_matrix, boundary_vector, wake_panels, source_strengths, surface_velocities, lift_coefficient

## Linear-strength source and vorticity panel method

include("Aerodynamics/LinearVortexSource/LinearVortexSource.jl")
import .LinearVortexSource: total_velocity, source_velocity, vortex_velocity, vortex_influence_matrix, source_influence_matrix, neumann_boundary_condition, kutta_condition, two_point_matrix, linear_source_matrix, linear_vortex_matrix, constant_source_matrix, constant_source_boundary_condition

export total_velocity, source_velocity, vortex_velocity, vortex_influence_matrix, source_influence_matrix, neumann_boundary_condition, kutta_condition, two_point_matrix, linear_source_matrix, linear_vortex_matrix, constant_source_matrix, constant_source_boundary_condition

## Vortex lattice

include("Aerodynamics/VortexLattice/VortexLattice.jl")
import .VortexLattice: Horseshoe, VortexLatticeSystem, References, AbstractAxisSystem, Stability, Wind, Body, Geometry, streamlines, influence_coefficient, influence_matrix, boundary_condition, solve_system, transform, bound_leg, bound_leg_center, bound_leg_vector, r1, r2, Horseshoe, surface_velocity, surface_forces, surface_moments, nearfield_drag, geometry_to_wind_axes, geometry_to_stability_axes, stability_to_geometry_axes, wind_to_geometry_axes,  rate_coefficient, nearfield, farfield, farfield_forces, surface_velocities, surface_forces, surface_dynamics, surface_coefficients, nearfield_coefficients, farfield_coefficients, VortexRing, Freestream, velocity, body_frame_velocity, kinematic_viscosity, mach_number

export Horseshoe, VortexLatticeSystem, References, AbstractAxisSystem, Stability, Wind, Body, Geometry, streamlines, influence_coefficient, influence_matrix, boundary_condition, solve_system, transform, bound_leg, bound_leg_center, bound_leg_vector, r1, r2, Horseshoe, surface_velocity, surface_forces, surface_moments, nearfield_drag, geometry_to_wind_axes, geometry_to_stability_axes, stability_to_geometry_axes, wind_to_geometry_axes,  rate_coefficient, nearfield, farfield, farfield_forces, surface_velocities, surface_forces, surface_dynamics, surface_coefficients, nearfield_coefficients, farfield_coefficients, VortexRing, Freestream, velocity, body_frame_velocity, kinematic_viscosity, mach_number

## Profile drag estimation

include("Aerodynamics/profile_drag.jl")

export wetted_area_drag, profile_drag_coefficient, local_dissipation_drag

## Viscous airfoil analysis

include("Aerodynamics/ViscFoil/ViscFoil.jl")

import .ViscFoil: BoundaryLayer2D, solve_inviscid_doublets, solve_inviscid_vortices, defect_block, edge_velocities, solve_viscous_case

export BoundaryLayer2D, solve_inviscid_doublets, solve_inviscid_vortices, defect_block, edge_velocities, solve_viscous_case

## Cases

include("Aerodynamics/Cases/printing.jl")

export print_case, print_info, print_coefficients, print_derivatives

include("Aerodynamics/Cases/cases.jl")

export solve_case, spanwise_loading, triangle_connectivities, extrapolate_point_mesh

include("Aerodynamics/Cases/stability_cases.jl")

export solve_case_derivatives, longitudinal_stability_derivatives, longitudinal_stability_matrix, lateral_stability_derivatives, lateral_stability_matrix

include("Aerodynamics/Cases/foil_cases.jl")

## Structural analyses
#==========================================================================================#

include("Structures/Beams.jl")
import .Beams: Material, Tube, Beam, radii, area, moment_of_inertia, polar_moment_of_inertia, J_coeffs, Iyy_coeffs, Izz_coeffs, tube_stiffness_matrix, bending_stiffness_matrix, axial_stiffness_matrix, build_stiffness_matrix, solve_cantilever_beam, elastic_modulus, shear_modulus, yield_stress, density, principal_stress, torsional_stress, von_mises_stress

export Material, Tube, Beam, radii, area, moment_of_inertia, polar_moment_of_inertia, J_coeffs, Iyy_coeffs, Izz_coeffs, tube_stiffness_matrix, bending_stiffness_matrix, axial_stiffness_matrix, build_stiffness_matrix, solve_cantilever_beam, elastic_modulus, shear_modulus, yield_stress, density, principal_stress, torsional_stress, von_mises_stress

## Aerostructural analyses
#==========================================================================================#

include("Aerostructural/Aerostructural.jl")
import .Aerostructural: AerostructWing, make_beam_mesh, transform_stiffy, permute_stiffy, build_big_stiffy, adjacent_adder, section_moments, compute_loads, fem_load_vector, rotation_matrix, transfer_displacements, mesh_translation, mesh_rotation, new_horseshoes, solve_coupled_residual!, aerostruct_gauss_seidel

export AerostructWing, make_beam_mesh, transform_stiffy, permute_stiffy, build_big_stiffy, adjacent_adder, section_moments, compute_loads, fem_load_vector, rotation_matrix, transfer_displacements, mesh_translation, mesh_rotation, new_horseshoes, solve_coupled_residual!, aerostruct_gauss_seidel

## Post-processing
#==========================================================================================#

include("Tools/plot_tools.jl")

export plot_panels, plot_streamlines, plot_planform, plot_surface

end


text:
## Marker style

#md # ![FILE_NAME.png](images/FILE_NAME.png)

## ref https://jp.mathworks.com/matlabcentral/fileexchange/35231-matlab-plot-gallery-standard-plot-markers?focused=6792907&tab=example

using SpecialFunctions

x = 0:0.2:10;
y0 = besselj.(0,x);
y1 = besselj.(1,x);
y2 = besselj.(2,x);
y3 = besselj.(3,x);
y4 = besselj.(4,x);

plot(x, [y0 y1 y2 y3 y4],
    marker=[:circle :rect :star5 :diamond :hexagon],
    line=false)

#-

Plots.supported_markers()


text:
const dictRest          = Dict( "basic"     => (1,1),
                                "basic2"    => (2,2),
                                "basic3"    => (3,2),
                                "nCol"      => (4,2),
                                "ret1"      => (5,2),
                                "ret2"      => (6,3),
                                "ret3"      => (7,3),
                                "noCol2"    => (8,3)  )


text:
"""
    randomBigInt()

Compute a random `Int64` and return it coverted it to a `BigInt`.
"""
function randomBigInt()
    big(rand(Int64))
end

"""
    randomBigFloat()

Compute a random `Float64` and return it coverted it to a `BigFloat`.
"""
function randomBigFloat()
    big(rand(Float64))
end

"""
    randomBigRational

Compute two random `Int64`, convert them to `BigInt`, and return a `Rational` made with them.
"""
function randomBigRational()
    randomBigInt() // randomBigInt()
end


text:
"""
    bind!(dest::Signal, src::Signal)

For every update to `src` also update `dest` with the same value.
"""
function bind!(dest::Signal, src::Signal)
    push!(dest.binders, Signal(x -> dest(x), src))
end
export bind!

"""
    unbind!(dst::Signal,src::Signal)

Remove bindings from `src` to `dst` that were previously created using `bind`.

    unbind!(dst::Signal)

Remove all bindings that were previously created using `bind` that will cause `dst` to update.
"""
function unbind!(dst::Signal, src::Signal)
    b_idx = findfirst(b -> b.action.args[1] == src, dst.binders)
    b_idx === 0 && return
    binder = dst.binders[b_idx]
    detach(binder)
    filter!(dst.binders) do b
        b != binder
    end
end
export unbind!

function unbind!(dst::Signal)
    foreach(detach, dst.binders)
    empty!(dst.binders)
end

import Base.detach
"""
    detach(s::Signal)

Detach a signal `s` from the signal graph making it unreachable and ready for GC.
"""
function detach(s::Signal)
    parents = Iterators.filter(x -> isa(x, Signal), s.action.args)
    foreach(parents) do p
        Iterators.filter!(x -> x!=s, p.children)
    end
    unbind!(s)
    nothing
end