Initialize a struct without a type-parameterized field and add the field later


#1

I am trying to do something that may, or may not be possible. I want to use ForwardDiff, and create a Jacobian function. This function will be part of a struct. A parameter container is also part of this struct. Now, I want the Jacobian function to close over the entire struct, not just the parameter container. A short description of that is present here: https://github.com/JuliaDiff/ForwardDiff.jl/issues/296

Unfortuantely, I do not know how to do it. At the moment here is my try:

# Here f must be of the form: f(x) -> SVector (ONE ARGUMENT!)
function generate_jacobian_oop(f::F, x::X) where {F, X}
    # Test f structure:
    L = length(x)
    y = f(x)
    @assert typeof(y) <: SVector
    @assert length(y) == L
    # Setup config
    cfg = ForwardDiff.JacobianConfig(f, x)
    FDjac(x, p) = ForwardDiff.jacobian(f, x, cfg)
    return FDjac
end

mutable struct DiscreteDS{D, T<:Number, F, P, J}
    state::SVector{D,T}
    eom::F
    p::P
    jacob::J
    function DiscreteDS{D, T, F, P}(state::SVector{D,T}, eom::F, p::P) where {D, T, F, P}
       f = new(state, eom, p)
       reducedeom = (x) -> eom(x, f.p)
       f.jacob = generate_jacobian_oop(reducedeom, state)
       return f
   end
end

function DiscreteDS(state::SVector{D,T}, eom::F, p::P) where {D, T, F, P}
    return DiscreteDS{D, T, F, P}(state, eom, p)
end

@inline henon_eom(x, p) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
p = [1.4, 0.3]

ds = DiscreteDS(SVector(0.0, 0.0), henon_eom, p)

This doesn’t run at the moment, because it tells me that there are too few parameters specified in new.

But more importantly, I also do not know how to enforce the type parameter in the line f.jacob = generate_jacobian_oop(reducedeom, state)

Any help would be greatly appreciated!


#2

Just create jacob before new?


#3

How can I then close it over f.p, where f comes from new?


#4

UPDATE: I have managed to do this using 2 different structs. I would like to avoid this if possible, but it is a last resort:

mutable struct DiscreteLaw{D, T, F, P}
    state::SVector{D,T}
    eom::F
    p::P
end

mutable struct DiscreteDS{D, T, F, P, J}
    system::DiscreteLaw{D, T, F, P}
    jacobian::J
end

function DiscreteDS(state::SVector{D,T}, eom::F, p::P) where {D, T, F, P}
    system = DiscreteLaw(state, eom, p)
    reducedeom = (x) -> eom(x, system.p)
    jacob = generate_jacobian_oop(reducedeom, state)
    return DiscreteDS(system, jacob)
end


@inline henon_eom(x, p) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
p = [1.4, 0.3]

ds = DiscreteDS(SVector(0.0, 0.0), henon_eom, p)