Does it make sense to construct ODESystems with non-existant variables?

[Libbum]

I'm hitting up against this assertion error quite a lot.

Perhaps it's my miss-understanding of how to use SINDy correctly; but it's happened in more than one way now.

Current example: Playing around with a Rössler attractor using the standard parameter values. Test system looks like

function sys!(du, u, p, t)
    X, Y, Z = u
    a, b, c = p

    du[1] = - Y - Z
    du[2] = X + a*Y
    du[3] = b + Z*(X - c)
end

datasize = 30
tspan = (0.0f0, 3.0f0)
tsteps = range(tspan[1], tspan[2], length = datasize)
u0 = Float32[0.3, 0.5, 0.78]

p_ = Float32[0.2, 0.2, 5.7]

prob = ODEProblem(sys!, u0, tspan, p_)

From which I build a UODE case to 'learn' what the Y variable does in each equation

L = FastChain(FastDense(3, 32, tanh),FastDense(32, 32, tanh), FastDense(32, 3))
p = initial_params(L)

function dudt_(u, p, t)
    X, Y, Z = u
    a, b, c, = p_

    y = L(u,p)

    [y[1] - Z,
     X + y[2],
     b + Z*(X-c) + y[3]]
end

prob_nn = ODEProblem(dudt_,u0, tspan, p)

Which gives me a learned dataset

Xₙ = Float32[0.30009267 0.18289553 0.08949813 0.00889696 -0.06638512 -0.13437773 -0.20161553 -0.26835278 -0.33266342 -0.39145237 -0.44572976 -0.4981197 -0.54668456 -0.5909244 -0.628824 -0.65913796 -0.68393165 -0.6999598 -0.7114669 -0.71467036 -0.71099627 -0.6965863 -0.67767423 -0.64880353 -0.6104574 -0.5692104 -0.5193644 -0.4607151 -0.3988337 -0.33018124; 0.49992666 0.53692186 0.56088054 0.5779834 0.5844479 0.5897381 0.58248544 0.5705773 0.55083734 0.525924 0.49281082 0.4530428 0.4070489 0.3569066 0.30112132 0.23897965 0.17396559 0.104443274 0.032255013 -0.040149473 -0.11661649 -0.19195427 -0.26873833 -0.34198284 -0.41627797 -0.48611206 -0.55182093 -0.61592114 -0.67387533 -0.7259051; 0.7799829 0.45902893 0.27633953 0.16883376 0.10817303 0.07465119 0.055841275 0.047263786 0.040068798 0.038583606 0.03303134 0.034528494 0.034686353 0.03375564 0.032579035 0.031360127 0.031881377 0.030465402 0.030870909 0.029191166 0.032303184 0.03046857 0.031072713 0.03068368 0.031213444 0.030175872 0.03015326 0.03111261 0.030196927 0.032249533]
L̂ = Float32[-0.49630117 -0.5374406 -0.5616105 -0.57744056 -0.58307993 -0.5880346 -0.5810973 -0.5698522 -0.55059266 -0.52623427 -0.4931721 -0.45368993 -0.4076958 -0.3573648 -0.30131128 -0.2388681 -0.17372836 -0.10402171 -0.03196443 0.040302556 0.11614147 0.19122498 0.2672777 0.34007347 0.4139306 0.48327568 0.5487542 0.6127155 0.6707917 0.7230119; 0.09663023 0.10830986 0.1130031 0.11535166 0.11610068 0.11740245 0.11637088 0.11456688 0.11085699 0.10594397 0.09890834 0.090527125 0.08076144 0.07020234 0.058583632 0.045802265 0.032851975 0.019072365 0.005285439 -0.00857036 -0.02254451 -0.03704656 -0.05139447 -0.06569367 -0.080529824 -0.094920576 -0.109131694 -0.12360385 -0.13766418 -0.15090223; 0.0076180287 0.005571019 0.0036096796 0.0024604686 0.0017505623 0.0012024231 0.0008894168 0.00075354055 0.0008024089 0.0010130554 0.0013005063 0.0016464256 0.0019050539 0.0019929819 0.0018744282 0.0015147477 0.0008264594 -5.383417f-5 -0.0012840256 -0.0025322922 -0.0041975193 -0.005404722 -0.006792832 -0.007694889 -0.008312013 -0.008508753 -0.00819264 -0.0073981173 -0.005996678 -0.004254561]

Using the basis from the SINDy example documentation, I recover

opt = SR3(1e-2)
Ψ = SINDy(Xₙ, L̂, basis, opt; maxiter = 10000, normalize = true, denoise = true)

julia> println(Ψ)
Sparse Identification Result
No. of Parameters : 2
Active terms : 2
   Equation 1 : 1
   Equation 2 : 1
   Equation 3 : 0
Overall error (L2-Norm) : 0.07513197
   Equation 1 : 0.010420825
   Equation 2 : 0.039975792
   Equation 3 : 0.024735346
AICC :
   Equation 1 : 20.70024
   Equation 2 : 15.32237
   Equation 3 : 14.940946

SR3{Float64,UnionAll}(0.01, 1.0, ProximalOperators.NormL1) converged after 5473 iterations.

julia> print_equations(Ψ)
2 dimensional basis in ["x", "y", "z"]
f_1 = p₁ * y
f_2 = p₂ * y

Which is exactly what I expect. But since there is no result in z / f_3, then creating an ODESystem fails:

julia> unknown_sys = ODESystem(Ψ)
ERROR: AssertionError: length(b) == length(variables(b))
Stacktrace:
 [1] ODESystem(::Basis{Array{Operation,1},Array{Operation,1},Array{Operation,1},Operation}) at /home/tim/.julia/packages/DataDrivenDiffEq/l4Hqc/src/basis.jl:442
 [2] ODESystem(::SparseIdentificationResult) at /home/tim/.julia/packages/DataDrivenDiffEq/l4Hqc/src/sindy/results.jl:154
 [3] top-level scope at REPL[580]:1

This comes about so that eqs = dvs .~ b(vs, parameters(b), independent_variable(b)) can be a broadcasted operation. But is it possible/advisable to construct an f_3 = 0 or f_3 = p₃ case when a variable is not extant?

---

The other case where I come up against this problem can be discussed using the same example. Consider we know that our unknown does not affect the third equation, but don't know which of the three variables of our system cover the right-most term in the first and second equations. I would write something like:

L = FastChain(FastDense(3, 32, tanh),FastDense(32, 32, tanh), FastDense(32, 2))
p = initial_params(L)

function dudt_(u, p, t)
    X, Y, Z = u
    a, b, c, = p_

    y = L(u,p)

    [-Y + y[1],
     X + y[2],
     b + Z*(X-c)]
end

prob_nn = ODEProblem(dudt_,u0, tspan, p)

Notice here there are three input variables into the NN chain, but only two outputs. Since we're not even looking for a result for the third equation here, calling ODESystem on the result will again hit the assertion error.

If these problems don't really make sense, it'd be nice to discuss how to parse a problem like the above in the documentation, to avoid any future questions of this nature.

[AlCap23]

You can extract the rhs from the SINDy result and turn it into a separate equation manually via MTK. I guess if we automate this, we would run into an error creating the iip function, but I didn't tried this yet.

A workaround would be to allow users to modify the result via pushing eqs.

But I am open for discussion here ;)

[Libbum]

Thanks @AlCap23. For the moment I'd be happy to manually do it, if it seems like what I need here is an edge case.

What (where) is the iip function? I'm not familiar enough with the process or codebase to know what you mean by that. I guess it could be this section in the basis function?

If so, you're hunch about failing there is probably correct.

If I manually build an ODESystem:

@derivatives D'~independent_variable(b)

vs = similar(b.variables)
dvs = similar(b.variables)

for (i, vi) in enumerate(b.variables)
    vs[i] = ModelingToolkit.Operation(vi.op, [independent_variable(b)])
    dvs[i] = D(vs[i])
end

expr = b(vs, parameters(b), independent_variable(b))
# Hardcoded in this case, since I know z=0
eqs = dvs[1:2] .~ expr[1:2]
push!(eqs, dvs[3] ~ 0)

unknown_sys = ODESystem(eqs, independent_variable(b), variables(b), parameters(b))

Then, I cannot construct a basis from it:

unknown_eq = ODEFunction(unknown_sys)
julia> b = Basis((u, p, t)->unknown_eq(u, [1.; 1.], t), u)
ERROR: MethodError: no method matching Basis(::Array{Expression,1}, ::Array{Operation,1}; parameters=Operation[], iv=t)

Edit: Same issue comes up if I build completely from scratch:

@parameters t p[1:3]
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ p[1] * y,
       D(y) ~ p[2] * y,
       D(z) ~ p[3]]

unknown_sys = ODESystem(eqs)
unknown_eq = ODEFunction(unknown_sys)
b = Basis((u, p, t)->unknown_eq(u, [1.; 1.], t), u) # Error

Although if D(z) ~ p[3] * x (for example), there is no problem.

[AlCap23]

I'll try to have a look later on today.

[Libbum]

This may end up solving itself once updating to MTK 4.0. I think it has something to do with SciML/ModelingToolkit.jl#433, which all the breaking changes in 518 squashed. The differences between an Expression, Operation and Variable are now mostly removed. I think that Constant(0) is also no longer needed, which turns the set of Operations in the above case into an Expression.

---

Upon sleeping on this, of course if I want to do this manually I don't need to use a value in the third equation, so

@parameters t p[1:2]
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ p[1] * y,
       D(y) ~ p[2] * y]

recovered_sys = ODESystem(eqs)
recovered_eq = ODEFunction(recovered_sys)

function dudt(u, p, t)
    X, Y, Z = u
    a, b, c, = p_

    y = recovered_eq(u, p, t)
    #TODO: The result is flipped? Why?

    [y[2] - Z,
     X + y[1],
     b + Z*(X-c)]
end

works just fine. Apart from the oddity that the order of results end up reversed—no idea why.

[AlCap23]

Hey! Sorry for not being too responsive. I had a lot on my plate.

Soo, starting from your recent example with slight modifications

using ModelingToolkit

@parameters t p[1:2]
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ p[1] * y,
       D(y) ~ p[2] * y]

recovered_sys = ODESystem(eqs)
recovered_eq = ODEFunction(recovered_sys)

function dudt(u, p, t)
    X, Y, Z = u
    a, b, c = p

    y = recovered_eq(u, p, t)
    #TODO: The result is flipped? Why?

    [y[2] - Z,
     X + y[1],
     b + Z*(X-c)]
end

@parameters x[1:3] p[1:3] t # For human readable tracing...

recovered_eq(x, p, t)
dudt(x, p, t)

Results in

2-element Array{Operation,1}:
 p₁ * x₂
 p₂ * x₂

3-element Array{Operation,1}:
        p₂ * x₂ - x₃
        x₁ + p₁ * x₂
 p₂ + x₃ * (x₁ - p₃)

Which checks out. I have been using ModelingToolkit v3.21.0.

[AlCap23]

Your first example is indeed an issue. I assumed that a Basis should always be created with Array{Operation}. Since the output of the equations is an Expression, simply converting this works.

using DataDrivenDiffEq
@parameters t p[1:3]
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ p[1] * y,
       D(y) ~ p[2] * y,
       D(z) ~ p[3]]

unknown_sys = ODESystem(eqs)
unknown_eq = ODEFunction(unknown_sys)

@parameters u[1:3]

b = Basis((u, p, t)->Operation.(unknown_eq(u, [1.0; 1.0; 1.0], t)), u) # Error
b.basis

With an output of

2 dimensional basis in ["u₁", "u₂", "u₃"]

2-element Array{Expression,1}:
                            u₂
 ModelingToolkit.Constant(1.0)

Which is due to the duplicate equations p[1]*y = 1.0*y and `p[2]y = 1.0y which are consequently removed.

[AlCap23]

Circling back to the original question: I think we should allow something like a recover_basis(result) or similar to extract partly complete results. Maybe this work? :)

[Libbum]

Hey @AlCap23, thanks for getting back to me. I'm swamped myself, so will take a little to absorb this and respond. Certainly a low priority! :)