Unable to reconstruct system from result

[ctessum]

When following this SINDy example, a logical final step would be to run a simulation with the discovered system of equations to inspect the result that it gives. One of the steps to do this would be to replace c[1] (the placeholder value for the forcing function) with the actual forcing function in the resulting equations. From a previous version of the documentation that seems like it may have disappeared, it seems that the way to do this would be something like:

t = get_iv(system)
subs_control = Dict(
    c[1] => exp(-((t - 5.0) / 5.0)^2)
)

eqs = map(equations(system)) do eq
    eq.lhs ~ substitute(eq.rhs, subs_control)
end

However, this gives the error ERROR: MethodError: <ₑ(::Num, ::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}) is ambiguous. Candidates: <ₑ(a::Number, b::SymbolicUtils.Symbolic) in SymbolicUtils at /.../.julia/packages/SymbolicUtils/qulQp/src/ordering.jl:9

Does anyone have any suggestions of how to fix this? Thanks!

[[deps.DataDrivenDiffEq]]
...
git-tree-sha1 = "52b8cdc6a05707d4385bba499653955a16466b86"
uuid = "2445eb08-9709-466a-b3fc-47e12bd697a2"
version = "1.0.0"

[[deps.Symbolics]]
...
git-tree-sha1 = "718328e81b547ef86ebf56fbc8716e6caea17b00"
uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
version = "4.13.0"

[AlCap23]

Yes! I'll work on that. I just ported the necessary stuff right now, but this is a good point. The missing documentation snippet can be found here, in the previous docs.

I think you need to unwrap the Num in this case, but I'll check out this evening.

[ctessum]

Thanks. I can confirm that the method in the previous doc worked with the previous version of the library, but it doesn't seem to work with v1.0.0.

[ctessum]

I was able to figure it out using Symbolics.unwrap, and also changing u and c to u(t) and c(t). Here is a full working example:

using DataDrivenDiffEq
using LinearAlgebra
using OrdinaryDiffEq
using DataDrivenSparse
using StableRNGs
using Symbolics
using Plots

rng = StableRNG(1337)

function pendulum(u, p, t)
    x = u[2]
    y = -9.81sin(u[1]) - 0.3u[2]^3 - 3.0 * cos(u[1]) - 10.0 * exp(-((t - 5.0) / 5.0)^2)
    return [x; y]
end

u0 = [0.99π; -1.0]
tspan = (0.0, 15.0)
prob = ODEProblem(pendulum, u0, tspan)
sol = solve(prob, Tsit5(), saveat = 0.01);

X = sol[:, :] + 0.2 .* randn(rng, size(sol));
ts = sol.t;

prob = ContinuousDataDrivenProblem(X, ts, GaussianKernel(),
                                   U = (u, p, t) -> [exp(-((t - 5.0) / 5.0)^2)],
                                   p = ones(2))

@parameters t
@variables u(t)[1:2] c(t)[1:1]
@parameters w[1:2]
u = collect(u)
c = collect(c)
w = collect(w)

h = Num[sin.(w[1] .* u[1]); cos.(w[2] .* u[1]); polynomial_basis(u, 5); c]

basis = Basis(h, u, parameters = w, controls = c);
println(basis) # hide

sampler = DataProcessing(split = 0.8, shuffle = true, batchsize = 30, rng = rng)
λs = exp10.(-10:0.1:0)
opt = STLSQ(λs)
res = solve(prob, basis, opt,
            options = DataDrivenCommonOptions(data_processing = sampler, digits = 1))

system = get_basis(res)
params = get_parameter_map(system)
equations(system)
println(system) # hide
println(params) # hide

t = Symbolics.unwrap(get_iv(system))
subs_control = Dict(
    c[1] => exp(-((t - 5.0) / 5.0)^2)
)

eqs = map(equations(system)) do eq
    eq.lhs ~ substitute(eq.rhs, subs_control)
end

@named sys = ODESystem(
    eqs,
    get_iv(system),
    states(system),
    parameters(system)
    );

x0 = [u[1] => u0[1], u[2] => u0[2]]
ps = get_parameter_map(system)

ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);
plot(estimate)