Sparse variable arrays and creation of ODESystem's

Specific Domains Modelling & Simulations
,
1

Sketch of the problem: a dumb example

I have a question about the ability of ‘enforcing’ sparsity of a specific variable array that is created with the @variables macro.

When running experiments, I have noticed that creating (in fact, completing) an ODESystem is slow when a specific variable is very large. The reason for this, I believe, is that when calling complete(sys), there is a for-loop over all individual elements of the variables in the variable array.

To illustrate, consider the following stupid examples

function create_odesystem(S::Int)
    @variables (x(t))[1:S]
    @parameters p[1:S]

    eqns = [D(x[i]) ~ p[i]*x[i] for i in 1:S]
    @named sys = ODESystem(eqns, t, [x...], [p])
    return complete(sys)
end

function create_largeodesystem(S::Int)
    @variables (x(t))[1:S]
    @parameters p[1:S,1:S,1:S]

    eqns = [D(x[i]) ~ p[i,i,i]*x[i] for i in 1:S]
    @named sys = ODESystem(eqns, t, [x...], [p])
    return complete(sys)
end

and the benchmarks

julia> S = 100;

julia> sys = create_odesystem(S);

julia> lsys = create_largeodesystem(S);

julia> @btime create_odesystem($S);
  53.522 ms (623281 allocations: 30.24 MiB)

julia> @btime create_largeodesystem($S);
  7.394 s (76616496 allocations: 4.13 GiB)

I hope that these examples sketch my issue: creating the ODESystem is very slow for large systems where parameters are large. While this is to be expected, in many cases, I know that the variable is high-dimensional but incredibly sparse. As in the example, only the diagonal S elements are used, so all other elements can, in principle, be omitted.

Main question

That being said, say that you are given a sparse structure (either two-dim., using SparseArrays, or perhaps k-dimensional using SparseArrayKit), is there a way that I can programmatically create variables, or a variable array, that only creates the variables that are non-zero?

Other related questions

  1. Is there a way to efficiently change the value of the sparse parameters, most preferably using setp?
    In the above example, once the parameter values and init. conditions are chosen, one can use ModelingToolkit.setp(prob, p) to set the parameters, even when these are matrices and/or sparse arrays. Is this still possible when the variable is generated in another way?
  2. Can the for-loop over entries of a variable array be avoided when calling complete(sys)?

I hope that the problem and the question are clear. I have scoured the documentation of both ModelingToolkit and Symbolics, but I could not find anything related to large systems with sparse variable arrays. If it does exist, I would be very happy if you can point me towards this part of the docs.

2

Note that this is a limitation of the Symbolic code generation. I recommend using JuliaSimCompiler with its alternative IR-based code generation for sufficiently large systems.

Yeah we currently cannot specialize on that.

Open an issue in Symbolics.jl

Not entirely? There is a correctness check that does need to be employed. That said, JuliaSimCompiler has an approach that is much faster, but with a few trade-offs.

Currently I don’t think we have a model for sparse arrays like that. It’s a good point. Open an issue in Symbolics.jl.