Interpolation within ModelingToolkit framework

New to Julia
1

Hello,

Is there a way to interpolate within a ModelingToolkit’s (MTK) component model using the parameters/variables of the model? Interpolations.jl seems to struggle with Num types, and from the example in the documentation (Composing Ordinary Differential Equations · ModelingToolkit.jl) I understand that it is only possible to use pre-interpolated data that is then fed to the model as a time-variable function (?).

To illustrate this idea, I include a simple example using t as coordinate to obtain the interpolated variable.

using ModelingToolkit
using DifferentialEquations
using Interpolations


# Function to interpolate a time series
function interpolateTimeSeries(t, time, values)

	interpolator = CubicSplineInterpolation(time, values)
	interpolatedValue = interpolator(t)

end

@register interpolateTimeSeries(t, time, values)

# Source to produce the output of a timeseries table
function lookupTable(;name, time, values)
	
	@parameters t
	@variables output(t)

	equations = [output ~ interpolateTimeSeries(t, time, values)]

	ODESystem(equations, t, output, []; name=name)

end

############################## lookupTable source model #################
# Data
time = LinRange(1, 100, 100)
values = randn(100)

# Instantiate lookup-table source model
@named sourceLookupTable = lookupTable(time = time, values = values)
simpSource = structural_simplify(sourceLookupTable)
u0 = [sourceLookupTable => ]

# Define the problem
prob = ODEProblem(simpSource, u0, (0.0, 100.0))
sol = solve(prob, Tsit5())

This produces the following error:

ERROR: LoadError: TypeError: non-boolean (Num) used in boolean context
Stacktrace:
 [1] inbounds_index at $PATH\.julia\packages\Interpolations\qHlUr\src\extrapolation\extrapolation.jl:108 [inlined]
 [2] inbounds_position at $PATH\.julia\packages\Interpolations\qHlUr\src\extrapolation\extrapolation.jl:99 [inlined]
 [3] Extrapolation at $PATH\.julia\packages\Interpolations\qHlUr\src\extrapolation\extrapolation.jl:48 [inlined]
 [4] interpolateTimeSeries(::Num, ::LinRange{Float64}, ::Array{Float64,1}) at $PATH\SourcesSinks.jl:10
 [5] lookupTable(; name::Symbol, time::LinRange{Float64}, values::Array{Float64,1}) at $PATH\SourcesSinks.jl:22
 [6] top-level scope at $PATH\SourcesSinks.jl:34
 [7] include(::String) at .\client.jl:457
 [8] top-level scope at REPL[7]:1
in expression starting at $PATH\SourcesSinks.jl:34
2

That works IIRC? But even if it doesn’t, you can just @register your function.

4

Thanks for the quick reply @ChrisRackauckas. I did @register my function in the previous example (maybe not correctly?).

I also tried DataInterpolations.jl as you suggested, both with and without an external registered function, but I get the same error as in the example included in the first post.

Without external function

using ModelingToolkit
using DifferentialEquations
using Interpolations
using DataInterpolations


function lookupTable(;name, time, values)
	
	@parameters t
	@variables output(t)
 	interpolator = CubicSpline(values, time)

	equations = [output - interpolator(t)]

	ODESystem(equations, t, [output], []; name=name)

end

# Data
time = LinRange(1, 100, 100)
values = randn(100)

# Instantiate lookup-table source model
@named sourceLookupTable = lookupTable(time = time, values = values)
simpSource = structural_simplify(sourceLookupTable)

With external registered function

# Function to interpolate a time series
function interpolateTimeSeries(t, time, values)

	interpolator = CubicSpline(values, time)
	interpolatedValue = interpolator(t)

end

@register interpolateTimeSeries(t, time, values)

# Source to produce the output of a timeseries table
function lookupTable(;name, time, values)
	
	@parameters t
	@variables output(t)

	equations = [output ~ interpolateTimeSeries(t, time, values)]

	ODESystem(equations, t, [output], []; name=name)

end

############################## lookupTable source model #################
# Data
time = LinRange(1, 100, 100)
values = randn(100)

# Instantiate lookup-table source model
@named sourceLookupTable = lookupTable(time = time, values = values)
#@named sourceLookupTable = lookupTable()
simpSource = structural_simplify(sourceLookupTable)

In both cases, and the first example using Interpolations.jl, I get the same error:

ERROR: LoadError: TypeError: non-boolean (Num) used in boolean context

These errors originate from the use of @parameters as an input during the interpolation, e.g. interpolator(t) or interpolateTimeSeries(t, time, values) .

5

You just forgot to specify the non-symbolic parts arguments. The fix is:

@register interpolateTimeSeries(t, time::AbstractVector, values::AbstractVector)
1 Like
6

Many thanks for your help @ChrisRackauckas. The work of the MKT team in this package is fantastic.

Just some constructive feedback. I feel that the @register functionality might be useful to many people. However, a quick search of @register in the search docs bar only leads to Composing Ordinary Differential Equations - Specifying a time-variable forcing function, where it says:

" MTK allows to “register” arbitrary Julia functions, which are excluded from symbolic transformations but are just used as-is. So, you could, for example, interpolate a given time series using DataInterpolations.jl."

From this description, at least for a new user of MTK like me, it is a bit difficult to infer that I can specify some variables as non-symbolic and the remaining as symbolic. I took the liberty of using the example included in the documentation and include a function that interpolates some time-series data with a mix of non-symbolic and symbolic inputs. Here is the code:

using ModelingToolkit
using DifferentialEquations
using DataInterpolations 


@parameters tau
@variables t, x(t), f(t)
D = Differential(t)

function f_interpolate(t, table_t, table_u)

	interpolator = QuadraticInterpolation(table_u, table_t)
	output = interpolator(t)

end

@register f_interpolate(t, table_t::AbstractVector, table_u::AbstractVector)

table_t, table_u = LinRange(0.5, 10.5, 10), randn(10)

@named fol_external_f = ODESystem([f ~ f_interpolate(t, table_t, table_u), D(x) ~ (f - x) / tau])
prob = ODEProblem(structural_simplify(fol_external_f), [x => 0.0], (0.0, 10.0), [tau => 0.75])

sol = solve(prob)
plot(sol, vars=[x, f])
scatter!(table_t, table_u)

Capture3

If you feel this may complement the existing documentation, feel free to add it. Otherwise it is left here so it may help future new users.

7

That would be a good tutorial: Mixing Symbolic Models with Data, or something.

2 Likes
8

This is just a performance tip but I would have f_interpolate return the interpolating function. Right now every time f_interpolate is called it constructs a interpolation object. I would rewrite it as follows

f_interpolate = QuadraticInterpolation(table_u, table_t)
@register f_interpolate(t)

Assuming @register accepts DataInterpolations

9

It accepts any Julia function.

10

Agree. Thanks for the remark. The code above can become inefficient when one considers arrays.

11

I tried some of the above answers, but I too am running into errors trying to include CubicSpline interpolation from DataInterpolations.jl in my loss function. In the example below I am interpolating the model onto the same grid just as a test. Thank you in advance!

julia> using ModelingToolkit, DataInterpolations, Polynomials

julia> @variables c0 c1
2-element Vector{Num}:
 c0
 c1

julia> @parameters x[1:10] data[1:10]
2-element Vector{Symbolics.Arr{Num, 1}}:
 x[1:10]
 data[1:10]

julia> function f_interpolate(t, table_t, table_u)

               interpolator = CubicSpline(table_u, table_t)
               output = interpolator.(t)

       end
f_interpolate (generic function with 1 method)

julia> rms(x,y) = sqrt(sum((x .- y).^2) / length(x))
rms (generic function with 1 method)

julia> typeof(x)
Symbolics.Arr{Num, 1}

julia> @register f_interpolate(t::Vector{Symbolics.Arr{Num, 1}}, table_t::Vector{Symbolics.Arr{Num, 1}}, table_u::Vector{Symbolics.Arr{Num, 1}})

julia> loss = rms(data, f_interpolate(x, x, Polynomial([c0, c1]).(x)))
ERROR: MethodError: no constructors have been defined for Any
Stacktrace:
 [1] _broadcast_getindex_evalf
   @ ./broadcast.jl:670 [inlined]
 [2] _broadcast_getindex
   @ ./broadcast.jl:643 [inlined]
 [3] getindex
   @ ./broadcast.jl:597 [inlined]
 [4] copy
   @ ./broadcast.jl:899 [inlined]
 [5] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, Type{Any}, Tuple{Vector{Term{Any, Nothing}}}})
   @ Base.Broadcast ./broadcast.jl:860
 [6] munge_data(u::Symbolics.Arr{Any, 1}, t::Symbolics.Arr{Num, 1})
   @ DataInterpolations ~/.julia/packages/DataInterpolations/HJsu4/src/interpolation_utils.jl:44
 [7] CubicSpline(u::Symbolics.Arr{Any, 1}, t::Symbolics.Arr{Num, 1})
   @ DataInterpolations ~/.julia/packages/DataInterpolations/HJsu4/src/interpolation_caches.jl:161
 [8] f_interpolate(t::Symbolics.Arr{Num, 1}, table_t::Symbolics.Arr{Num, 1}, table_u::Symbolics.Arr{Any, 1})
   @ Main ./REPL[4]:3
 [9] top-level scope
   @ REPL[9]:1

julia>
12

Hi,
I’m getting the same error when I try to interpolate 2 (or more) dimensional data even if I register the function.

Hereafter a simplified version of my code for a 4D interpolation.

using ModelingToolkit
using Interpolations

@variables t

function interpolateMatrix(X1,X2,X3,X4,Data,bp1,bp2,bp3,bp4)

    interpolator = linear_interpolation((bp1, bp2, bp3, bp4), Data)
    output = interpolator(X1, X2, X3, X4)

end

@register interpolateMatrix(X1,X2,X3,X4,Data::AbstractVector,bp1::AbstractVector,bp2::AbstractVector,bp3::AbstractVector,bp4::AbstractVector)


function lookup_table_4D(Data,bp1,bp2,bp3,bp4;name)

    @variables X1(t), X2(t), X3(t), X4(t), output(t)

    eqs = [output ~ interpolateMatrix(X1,X2,X3,X4,Data,bp1,bp2,bp3,bp4)]

    ODESystem(eqs;name)
    
end

A = randn(9,9,9,4)
X1 = 1:9
X2 = 1:9
X3 = 1:9
X4 = 1:4

@named matA = lookup_table_4D(A, X1, X2, X3, X4)

The error I get is

ERROR: TypeError: non-boolean (Num) used in boolean context
Stacktrace:
 [1] inbounds_index
   @ $PATH\.julia\packages\Interpolations\nDwIa\src\extrapolation\extrapolation.jl:111 [inlined]
 [2] inbounds_position
   @ $PATH\.julia\packages\Interpolations\nDwIa\src\extrapolation\extrapolation.jl:102 [inlined]
 [3] Extrapolation
   @ $PATH\.julia\packages\Interpolations\nDwIa\src\extrapolation\extrapolation.jl:48 [inlined]
 [4] interpolateMatrix(X1::Num, X2::Num, X3::Num, X4::Num, Data::Array{Float64, 4}, bp1::UnitRange{Int64}, bp2::UnitRange{Int64}, bp3::UnitRange{Int64}, bp4::UnitRange{Int64})
   @ Main $PATH\4D_lut.jl:10
 [5] lookup_table_4D(Data::Array{Float64, 4}, bp1::UnitRange{Int64}, bp2::UnitRange{Int64}, bp3::UnitRange{Int64}, bp4::UnitRange{Int64}; name::Symbol)
   @ Main$PATH\4D_lut.jl:21
 [6] top-level scope
   @ $PATH\.julia\packages\ModelingToolkit\oE1IR\src\systems\abstractsystem.jl:923

Any idea on how to solve this problem ?

1 Like
13

Can you open an issue for that? DataInterpolations.jl should be fine, it has overloads and such. Interpolations.jl, not so much. It might need a different registration given it’s weird constructors.

14

I have encountered the same problem just yesterday also using Interpolations.jl, but I have overcome it by a combination of extrapolate and interpolate functions, together with @register, as it seemed that the error above indicates that there is a check on whether a symbol is within the bounds of the interpolating function. However, I have a complete beginner with the ModelingToolkit and Julia and I am not sure if my solution is correct.

I would also be willing to follow-up and raise the issue for that, should it be done at GitHub - JuliaMath/Interpolations.jl: Fast, continuous interpolation of discrete datasets in Julia?

15

The problem was that Data isn’t a vector. The correct syntax for the register is

@register interpolateMatrix(X1,X2,X3,X4,Data::AbstractArray,bp1::AbstractVector,bp2::AbstractVector,bp3::AbstractVector,bp4::AbstractVector)

The solution here : #2039

16

However, I get the same error when X4 is a vector instead of a scalar. I have tried to give the type ::Vector{Float64} or ::AbstractVector for this variable but none of them is working.

Here is the code:

using ModelingToolkit
using Interpolations
using Symbolics: scalarize

@variables t
δ = Differential(t)

function interpolateMatrix(X1,X2,X3,X4,Data,bp1,bp2,bp3,bp4)

    interpolator = linear_interpolation((bp1, bp2, bp3, bp4), Data)
    output = interpolator(X1, X2, X3, X4)
    
end

@register interpolateMatrix(X1,X2,X3,X4::Vector{Float64},Data::AbstractArray,bp1::AbstractVector,bp2::AbstractVector,bp3::AbstractVector,bp4::AbstractVector)::Vector{Float64}

function lookup_table_4D(Data,bp1,bp2,bp3,bp4;name)

    @parameters X4 = 1:4
    @variables X1(t) X2(t) X3(t) output(t)[1:4]
    
    eqs = [scalarize(output .~ interpolateMatrix(X1,X2,X3,X4,Data,bp1,bp2,bp3,bp4));]

    ODESystem(eqs;name)
end

A = randn(2,2,2,4)
bp1 = 1:2
bp2 = 1:2
bp3 = 1:2
bp4 = 1:4

@named matA = lookup_table_4D(A, bp1, bp2, bp3, bp4)

I have also noticed that the registered function interpolationMatrix has 19683 methods.

julia> interpolateMatrix
interpolateMatrix (generic function with 19683 methods)

I suppose that’s because the register generates a method for each combination of inputs’ type. However, if the types of inputs and parameters are known and fixed, only one of these methods is useful (ex: X#::Vector{Float64} ,Data::Array{Float64, 4} and bp#::Vector{Float64}). Is there a way to register only the required method?

17

I tried some of the codes above and they work, but the performance and number of allocations are very high how can I reduce it?

18

Hard to say without seeing your code, but this comment above is good advice.