Skip to content

Sparse Identification of Nonlinear Dynamics example in the tutorial fails #294

@hurak

Description

@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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions