Register Array Symbolic Output Indexing

I’m looking to use an array function in a ModelingToolkit.jl model through the @register_array_symbolic macro, but I’m having trouble figuring out how to use it properly.

The following MWE shows what I’m trying to do:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

function simple_function(xdata)
    return xdata
end


@variables u(t)[1:2]
type1 = typeof(vcat(u))

@register_array_symbolic simple_function(xdata::type1) begin
    size = size(u,1)
    eltype = eltype(xdata)
end

eqs = [D(u[i]) ~ simple_function(vcat(u))[i] for i = 1:size(u,1)]

I want it to return a vector of equations that looks like this:

[D(u[1]) ~ simple_function(vcat(u))[1],
 D(u[2]) ~ simple_function(vcat(u))[2]]

but it returns the following error:

ERROR: BoundsError: attempt to access SymbolicUtils.BasicSymbolic{AbstractVector{Real}} at index [1]
Stacktrace:
 [1] getindex(x::SymbolicUtils.BasicSymbolic{AbstractVector{Real}}, idx::Int64)
   @ Symbolics \.julia\packages\Symbolics\EftHA\src\array-lib.jl:32
 [2] getindex(x::Symbolics.Arr{Num, 1}, idx::Int64)
   @ Symbolics \.julia\packages\Symbolics\EftHA\src\array-lib.jl:94
 [3] (::var"#11#12")(i::Int64)
   @ Main .\none:0
 [4] iterate
   @ .\generator.jl:48 [inlined]
 [5] collect(itr::Base.Generator{UnitRange{Int64}, var"#11#12"})
   @ Base .\array.jl:791
 [6] top-level scope
   @ Untitled-1:17

I’m not sure why a vector of reals cannot be indexed in this context, so any insight would be appreciated.

Thanks!

Hi! It seems you’re using the @register_array_symbolic macro incorrectly. It should express the properties (size, eltype, etc) of the output in terms of the input arguments. This is the correct syntax:

julia> @register_array_symbolic simple_function(xdata::Vector{Real}) begin
           size = (size(xdata,1),)
           eltype = eltype(xdata)
       end

Note how:

  • xdata::AbstractVector not ::type1. The type annotation should denote the type that the symbolic variable passed to simple_function should represent, not the actual type passed. This function takes a vector of real numbers. Symbolics is smart enough to figure out that you also want to be able to give it BasicSymbolic{Vector{Real}}, Symbolics.Arr{Num, 1}, Vector{Num}, etc. because all of those things represent a vector of reals.
  • size is a tuple
  • size is in terms of xdata not an arbitrary symbolic variable in the environment where the function is defined

This allows you to do:

julia> eqs = [D(u[i]) ~ simple_function(vcat(u))[i] for i = 1:size(u,1)]
2-element Vector{Equation}:
 Differential(t)((u(t))[1]) ~ (simple_function(SymbolicUtils.BasicSymbolic{Real}[(u(t))[1], (u(t))[2]]))[1]
 Differential(t)((u(t))[2]) ~ (simple_function(SymbolicUtils.BasicSymbolic{Real}[(u(t))[1], (u(t))[2]]))[2]

It’s also often useful to define ndims, so you could also define:

julia> @register_array_symbolic simple_function(xdata::Vector{Real}) begin
           size = (size(xdata,1),)
           eltype = eltype(xdata)
           ndims = 1
       end
1 Like

Thanks a lot! The input typing makes more sense to me now.

What is the purpose of ndims though? I’m not seeing anything about it in the docstring.

ndims is used for promote_symtype. It’s not something you should have to call, but it helps Symbolics figure out the number of dimensions of the output array without knowing the values passed to the registered function.

1 Like