Sparse Identification of Nonlinear Dynamics example in the tutorial fails

[hurak]

Sparse Identification of Nonlinear Dynamics example in the Nonlinear Systems tutorial fails (in Julia 1.6.3). Below I copy all the code and only show the REPL output for the last (failing) line:

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq
using Plots
using Random
using Symbolics: scalarize

Random.seed!(1111) # For the noise

# Create a nonlinear pendulum
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)

# Create the data with additional noise
X = sol[:,:] + 0.1 .* randn(size(sol))
DX = similar(sol[:,:])

for (i, xi) in enumerate(eachcol(sol[:,:]))
    DX[:,i] = pendulum(xi, [], sol.t[i])
end

ts = sol.t

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

@variables u[1:2] c[1:1]
@parameters w[1:2]
u = scalarize(u)
c = scalarize(c)
w = scalarize(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)

λs = exp10.(-10:0.1:-1)
opt = STLSQ(λs)

julia> res = solve(prob, basis, opt, progress = false, denoise = false, normalize = false, maxiter = 5000)
ERROR: BoundsError: attempt to access Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}} at index [6]
Stacktrace:
  [1] getindex(x::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}}, idx::Int64)
    @ Symbolics ~/.julia/packages/Symbolics/Cmx10/src/array-lib.jl:26
  [2] scalarize(term::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}}, idx::Tuple{Int64})
    @ Symbolics ~/.julia/packages/Symbolics/Cmx10/src/arrays.jl:497
  [3] scalarize(arr::Term{Real, Base.ImmutableDict{DataType, Any}})
    @ Symbolics ~/.julia/packages/Symbolics/Cmx10/src/arrays.jl:641
  [4] scalarize(arr::Num)
    @ Symbolics ~/.julia/packages/Symbolics/Cmx10/src/arrays.jl:643
  [5] (::Symbolics.var"#97#98"{Symbolics.Arr{Num, 1}})(i::Tuple{Int64})
    @ Symbolics ~/.julia/packages/Symbolics/Cmx10/src/arrays.jl:637
  [6] iterate
    @ ./generator.jl:47 [inlined]
  [7] collect_to!(dest::Vector{Num}, itr::Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}}}, Symbolics.var"#97#98"{Symbolics.Arr{Num, 1}}}, offs::Int64, st::Tuple{Tuple{Int64, Int64}})
    @ Base ./array.jl:728
  [8] collect_to_with_first!(dest::Vector{Num}, v1::Num, itr::Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}}}, Symbolics.var"#97#98"{Symbolics.Arr{Num, 1}}}, st::Tuple{Tuple{Int64, Int64}})
    @ Base ./array.jl:706
  [9] collect(itr::Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}}}, Symbolics.var"#97#98"{Symbolics.Arr{Num, 1}}})
    @ Base ./array.jl:687
 [10] map(f::Function, A::Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}}})
    @ Base ./abstractarray.jl:2323
 [11] scalarize(arr::Symbolics.Arr{Num, 1})
    @ Symbolics ~/.julia/packages/Symbolics/Cmx10/src/arrays.jl:636
 [12] build_parametrized_eqs(X::Matrix{Float64}, b::Basis)
    @ DataDrivenDiffEq ~/.julia/packages/DataDrivenDiffEq/AGBrC/src/solution.jl:133
 [13] build_solution(prob::DataDrivenProblem{Float64, false, DataDrivenDiffEq.Continuous}, Ξ::Matrix{Float64}, opt::STLSQ{Vector{Float64}}, b::Basis; eval_expression::Bool)
    @ DataDrivenDiffEq ~/.julia/packages/DataDrivenDiffEq/AGBrC/src/solution.jl:167
 [14] solve(p::DataDrivenProblem{Float64, false, DataDrivenDiffEq.Continuous}, b::Basis, opt::STLSQ{Vector{Float64}}; normalize::Bool, denoise::Bool, maxiter::Int64, round::Bool, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:progress,), Tuple{Bool}}})
    @ DataDrivenDiffEq ~/.julia/packages/DataDrivenDiffEq/AGBrC/src/solve/sindy.jl:53
 [15] top-level scope
    @ none:1

[anandijain]

if you want to make these issues productive, maybe just make a PR for jldoctests so that we know when doc breakages occur

[ChrisRackauckas]

Indeed this library should enable doctests

[hurak]

I am willing - in fact even eager - to contribute in this way (by formatting these examples as jldoctests), but I may need some little initial guidance. In particular, I do not know how to formulate the test when the correct result is currently not produced - against what shall I test then?

[AlCap23]

I am also new to doctests, but AFAIK we can test for the output of code blocks within the documentation. So something similar to res should be defined and a DataDrivenResult.

On the note of the failing docs example, I've boiled down a MWE for it. It seems related to the offset of the parameters I due to already present parameters in the Basis.

@parameters v[3:6]
collect(v)

Will result in

ERROR: BoundsError: attempt to access Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}} at index [5]
Stacktrace:
  [1] getindex(x::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}}, idx::Int64)
    @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/array-lib.jl:26
  [2] scalarize(term::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}}, idx::Tuple{Int64})
    @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/arrays.jl:497
  [3] scalarize(arr::Term{Real, Base.ImmutableDict{DataType, Any}})
    @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/arrays.jl:641
  [4] scalarize(arr::Num)
    @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/arrays.jl:643
  [5] (::Symbolics.var"#96#97"{Symbolics.Arr{Num, 1}})(i::Tuple{Int64})
    @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/arrays.jl:637
  [6] iterate
    @ ./generator.jl:47 [inlined]
  [7] collect(itr::Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}}}, Symbolics.var"#96#97"{Symbolics.Arr{Num, 1}}})
    @ Base ./array.jl:681
  [8] map(f::Function, A::Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}}})
    @ Base ./abstractarray.jl:2323
  [9] scalarize(arr::Symbolics.Arr{Num, 1})
    @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/arrays.jl:636
 [10] collect(x::Symbolics.Arr{Num, 1})
    @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/arrays.jl:653
 [11] top-level scope
    @ REPL[36]:1

I am using

     Status `~/.julia/environments/v1.6/Project.toml`
  [6e4b80f9] BenchmarkTools v1.2.0
  [2445eb08] DataDrivenDiffEq v0.6.6 `~/.julia/dev/DataDrivenDiffEq`
  [aae7a2af] DiffEqFlux v1.44.0
  [e30172f5] Documenter v0.27.10
  [5789e2e9] FileIO v1.11.2
  [587475ba] Flux v0.12.8
  [f6369f11] ForwardDiff v0.10.21
  [033835bb] JLD2 v0.4.15
  [23fbe1c1] Latexify v0.15.9
  [16fef848] LiveServer v0.7.0
  [bdcacae8] LoopVectorization v0.12.94
  [961ee093] ModelingToolkit v6.4.9
  [429524aa] Optim v1.5.0
  [1dea7af3] OrdinaryDiffEq v5.65.5
  [14b8a8f1] PkgTemplates v0.7.23
  [91a5bcdd] Plots v1.23.5
  [295af30f] Revise v3.1.20
  [0c5d862f] Symbolics v3.2.3

If I change it to use different variables, everything works out. I'll try out a newer Symbolics version, but It could be related to https://github.com/JuliaSymbolics/Symbolics.jl/issues/417

[AlCap23]

@shashi Updating to v3.5.1 did not improve. It seems that the offset is causing the array index to be kinda broken.

I've checked until here, and it seems that

v[1]

results in v[3] as an output. However, the Iterator cycles through 3..5. A quick fix would be to enumerate the idxs, but I don't think this will help in the global scope.

[AlCap23]

Closed via #304