Hey everyone. ComponentArrays.jl is my first package and I’m looking for some feedback and to gauge interest before I think about going down the path of registering it. I wrote this package mostly for myself for my work in control systems engineering where it’s important to be able to swap out components of your model (controllers, actuators, subsystems, etc.) quickly or pull them out and run them as standalone systems for debugging. This is easy to do in modeling languages like Simulink, but I wanted to do this with DifferentialEquations.jl.
MultiScaleArrays.jl serves the same general purpose and has some cool features that ComponentArrays.jl does not, such as the ability to add components on the fly (which, for my models at least, would be a neat way to handle rocket staging and spacecraft separation). The advantage of ComponentArrays.jl is simpler syntax and easy incorporation of existing models that weren’t written with ComponentArrays in mind. Composing models is just regular old function composition.
CArrays are the main exported type of ComponentArrays.jl. In their one-dimensional form, they act like arbitrarily nested NamedTuples that can be passed through a solver. I say “solver” and not “differential equation solver” because they work nicely with other types of solvers as well, such as those from Optim.jl. Here is what they look and act like:
julia> c = (a=2, b=[1, 2]);
julia> x = CArray(a=1, b=[2, 1, 4], c=c)
CArray{Float64}(a = 1.0, b = [2.0, 1.0, 4.0], c = (a = 2.0, b = [1.0, 2.0]))
julia> x.c.a = 400; x
CArray{Float64}(a = 1.0, b = [2.0, 1.0, 4.0], c = (a = 400.0, b = [1.0, 2.0]))
julia> x[5]
400.0
julia> collect(x)
7-element Array{Float64,1}:
1.0
2.0
1.0
4.0
400.0
1.0
2.0
julia> typeof(similar(x, Int32)) === typeof(CArray{Int32}(a=1, b=[2, 1, 4], c=c))
true
At the moment, all indexing, whether symbolic or numeric, uses views rather than slices because that’s what makes most sense for composable modeling. This might seem surprising if you try to slice into them like x[1:4]. Maybe I’ll play around with rules for that later, but it seems like the right way to handle it right now. I’ll also probably make a slice function and @slice macro at some point.
Higher dimensional CArrays can be created too, but it’s a little messy right now. It’s probably best to avoid creating them directly. The nice thing for modeling, though, is that dimension expansion through broadcasted operations can create higher-dimensional CArrays automatically, so Jacobian cache arrays that are created internally with false .* x .* x' will be CArrays with proper axes.
julia> x2 = x .* x'
7×7 CArray{Tuple{Axis{(a = 1, b = 2:4, c = (5:7, (a = 1, b = 2:3)))},Axis{(a = 1, b = 2:4, c = (5:7, (a = 1, b = 2:3)))}},Float64,2,Array{Float64,2}}:
1.0 2.0 1.0 4.0 2.0 1.0 2.0
2.0 4.0 2.0 8.0 4.0 2.0 4.0
1.0 2.0 1.0 4.0 2.0 1.0 2.0
4.0 8.0 4.0 16.0 8.0 4.0 8.0
2.0 4.0 2.0 8.0 4.0 2.0 4.0
1.0 2.0 1.0 4.0 2.0 1.0 2.0
2.0 4.0 2.0 8.0 4.0 2.0 4.0
julia> x2[:c,:c]
3×3 CArray{Tuple{Axis{(a = 1, b = 2:3)},Axis{(a = 1, b = 2:3)}},Float64,2,SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}}:
4.0 2.0 4.0
2.0 1.0 2.0
4.0 2.0 4.0
julia> x2[:a,:a]
1.0
julia> x2[:a,:c]
CArray{Float64}(a = 2.0, b = [1.0, 2.0])
julia> x2[:b,:c]
3×3 CArray{Tuple{Axis{NamedTuple()},Axis{(a = 1, b = 2:3)}},Float64,2,SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}}:
4.0 2.0 4.0
2.0 1.0 2.0
8.0 4.0 8.0
As you might have been able to see above, idea behind CArrays are to store the data as a dense array and abuse the value-as-type ability of julia’s type system and some @generated functions for indexing. The nice thing about doing it this way is large homogeneous arrays of CArrays (such as those that are built in differential equations solves) don’t have to all have to separately store the information needed to index into them. Or at least I think they don’t, right? It also makes broadcasting, similars and copy/copytos much faster than all of the other ways I tried (I went through about 6 iterations of this, including keeping the NamedTuple data structure and replacing its elements with views into an array).
Here’s an completely arbitrary and non-physical example of the use case. There are two subsystems governed by the Lorenz and Lotka-Volterra equations and a top-level system that handles coupling between the subsystems. Note that the composed system could have its own state variables or could be itself a component of a higher-level system, but it doesn’t and it isn’t because this is just a simple example.
using ComponentArrays
using DifferentialEquations
using Parameters: @unpack
tspan = (0.0, 20.0)
## Lorenz system
function lorenz!(D, u, (p, f), t)
@unpack σ, ρ, β = p
@unpack x, y, z = u
D.x = σ*(y - x)
D.y = x*(ρ - z) - y - f
D.z = x*y - β*z
return nothing
end
lorenz_p = (σ=10.0, ρ=28.0, β=8/3)
lorenz_ic = CArray(x=0.0, y=0.0, z=0.0)
lorenz_prob = ODEProblem(lorenz!, lorenz_ic, tspan, (lorenz_p, 0.0))
## Lotka-Volterra system
function lotka!(D, u, (p, f), t)
@unpack α, β, γ, δ = p
@unpack x, y = u
D.x = α*x - β*x*y + f
D.y = -γ*y + δ*x*y
return nothing
end
lotka_p = (α=2/3, β=4/3, γ=1.0, δ=1.0)
lotka_ic = CArray(x=1.0, y=1.0)
lotka_prob = ODEProblem(lotka!, lotka_ic, tspan, (lotka_p, 0.0))
## Composed Lorenz and Lotka-Volterra system
function composed!(D, u, p, t)
c = p.c #coupling parameter
@unpack lorenz, lotka = u
lorenz!(D.lorenz, lorenz, (p.lorenz, c*lotka.x), t)
lotka!(D.lotka, lotka, (p.lotka, c*lorenz.x), t)
return nothing
end
comp_p = (lorenz=lorenz_p, lotka=lotka_p, c=0.01)
comp_ic = CArray(lorenz=lorenz_ic, lotka=lotka_ic)
comp_prob = ODEProblem(composed!, comp_ic, tspan, comp_p)
## Solve problem
# We can solve the composed system...
comp_sol = solve(comp_prob)
# ...or we can unit test one of the component systems
lotka_sol = solve(lotka_prob)
There is an example in the examples folder that shows this same problem, but with composed Jacobians as well.
As far as speed goes, there is generally very little overhead penalty to using CArrays over regular Arrays. I’ve been benchmarking as I go along and don’t have anything formal yet, but that will come if there is interest in the package. The one issue is indexing like x[:a] rather than x.a is a little spotty speed-wise even though both are handled the same way. Sometimes the constant propagation works and sometimes it doesn’t and I don’t have a good feel for when or why. Since multidimensional CArrays don’t have the x.a option, you can use Vals as indexes like x[Val(:a), Val(:b)]. If anyone has some ideas on this, let me know. I might also just be benchmarking them incorrectly.