Best approach when the parameters object includes several types

[rvignolo]

Hi @ChrisRackauckas, hope you are having a great day.

What I am going to ask is not an issue, but I believe it can be useful for future reference.

I just wanted to know if you could point to me which of the following practices are the best when using DifferentialEquations.

For an ODE, the f(u, p, t) may depend on several parameters, such as scalars, vectors, matrices, functions. In this context, p cannot be a simple array with values of the same type.

Let's write a concrete example, where f(u, t) = a + h(t) + v[1] * u + v[2] * u * u. In this case, the object p can be a named tuple, so:

function f(du, u, p, t)

    a = p.a
    v = p.v
    h = p.h

    du = a + h(t) + v[1] * u + v[2] * u * u
end

On the other hand, pcan be an instance of a struct such as:

struct Container

    """The container variables."""
    variables  :: Dict{Symbol, Variable_t}

    """The container vectors."""
    vectors    :: Dict{Symbol, Vector_t}

end

yielding to the following function f:

function f(du, u, p, t)

    a = p.variables[:a].value
    v = p.vectors[:v].value

    du = a + h(t) + v[1] * u + v[2] * u * u
end

In this case, each time the solver calls f would involve searching elements on dictionaries (instead of tuples as in the previous example).

So, which one is the best approach or both are okay?

Another important difference between both examples is that in one case, a pointer to the function h is added to the tuple while in the other case, the Container struct does not contain such pointers. This is because I can define h and make it visible to f so I don't need to pass the pointer. Is that okay?

Thanks!

[ChrisRackauckas]

I would stay away from dictionaries if you can. Can't you just make a struct or a named tuple with the right variables with usable names?

[rvignolo]

As you suggest, dictionaries are much slower, so I am using named tuples instead.

Can't you just make a struct or a named tuple with the right variables with usable names?

Yes and no. Let me show you why.

So, I finally decided that in order to use your package (DifferentialEquations.jl) and take the most of it, I have to code my infrastructure in Julia. In that context, I am writing a couple of macros to describe a system of SDEs as simple as I can. I took some inspiration from JuMP.jl, ModelingToolkit.jl and Pumas.jl packages (among others). Imagine you just want to describe the following SDE system:

dr(t) = κ(t) * (θ(t) - r(t)) * dt + σ(t) * sqrt(α(t) + β(t) * r(t))) * dWᵣ(t),
dX(t) = r(t) * X(t) * dt + σₓ * X(t) * dWₓ(t),
with r(0) = r₀ and X(0) = X₀.

I have coded the following DSL to accomplish that task:

# This is a container of processes or SDEs.
dynamics = ModelDynamics()

# Just write each SDE as if you were writing on a piece of paper
@process(dynamics, r, r₀, κ(t) * (θ(t) - r), σ(t) * sqrt(α(t) + β(t) * r))
@process(dynamics, X, X₀, r * X, σₓ * X)

Now, specify each parameter in the previous equations and store them in a container of mathematical objects:

# This is a container of mathematical objects.
params = Container()

# define variables
@variable(params, r₀, 2.e-02)
@variable(params, X₀, 1.e+00)
# notice they can depend on other defined objects:
@variable(params, a, 5.e-03)
@variable(params, b, 5.e-03)
@variable(params, σₓ, a + b)

# define algebraic or interpolated functions
@functixn(params, κ(t), a + t)
@functixn(params, σ(t), b * t)
@functixn(params, θ(t), sin(t))
# define the remaining functions as well (α(t), β(t), ...)

Then, in order to simulate the system, just get the model dynamics and the parameters:

# algorithmic data should also be specified, such as number 
# of trajectories, time span, algorithm, seed, etc
@simulate(dynamics, params)

These variables, vectors, matrices, etc, are stored in the Container object, which has the following structure:

struct Container

    """The container variables."""
    variables  :: Dict{Symbol, Variable_t}

    """The container vectors."""
    vectors    :: Dict{Symbol, Vector_t}

    # ...
end

Since we have to append a new object to the corresponding dictionary each time a new mathematical object is defined, we cannot just define named tuples. However, before calling the SDE solver, we can easily translate the Container to a named tuple (this is what I am doing at the moment).

But there is a catch to it. Defined functions may depend on any other mathematical objects. So, when they are built, they have a pointer or reference to the container as argument, for example, function κ(t) above would be:

function κ(t, p=Ref(params))

    a = p[].variables[:a].value
    b = p[].variables[:b].value
    X₀ = p[].variables[:X₀].value
    rₓ = p[].variables[:rₓ].value
    σₓ = p[].variables[:σₓ].value

    return a + t
end

which can be called using only one argument, i.e., κ(t) for any t.The drift in process r has been given in that way by the user.

Also, more importantly, the definition of the function κ(t) would change if I change the value of the variable a at any point in the future. This behavior is wanted for optimization problems for example, when variables are changing and all mathematical objects that depend on those variables must also change. I cannot use a pointer to a named tuple as second argument in the function definition.

I think I have a workaround, by maybe you can provide some insight! Thanks a lot!

[ChrisRackauckas]

But this isn't anywhere near Pumas/JuMP/ModelingToolkit. Those all generate the equations at compile time and run on a flatted version. So no dictionaries of variables at runtime as parameters or any of that, and it can un-dynamically access NamedTuples because the compile-time gives a barrier at which the representation is then set.

[rvignolo]

Okay, I don't think I have explain myself. But never mind!

Just another question: Is it possible to simulate a process that has a diffusion that at time t depends not just in u(t) but also in u(ti) with ti a set of times such that ti < t for any i?

I wanted to simulate two processes, with the second process diffusion depending on the realized volatility of the first one.

[ChrisRackauckas]

Alright, you do you 🤷‍♂ .

Is it possible to simulate a process that has a diffusion that at time t depends not just in u(t) but also in u(ti) with ti a set of times such that ti < t for any i?

Sounds like you're looking for a stochastic delay differential equation? It's not done yet, but https://github.com/JuliaDiffEq/StochasticDelayDiffEq.jl is in progress. We'll add SDDEs to the docs when it's ready.

[rvignolo]

Amazing! thank you so much! I will be waiting for that package!

I will try to show you an example once I have a better development (by including your remarks) of my code. Thank you Chris!