Merge pull request #140 from SciML/sindy_api

SINDY Signature and rename

Author: AlCap23

README.md

@@ -69,7 +69,7 @@ end
 basis = Basis(polys, u)
 
 opt = STRRidge(0.1)
-Ψ = SInDy(X, DX, basis, maxiter = 100, opt = opt, normalize = true)
+Ψ = SINDy(X, DX, basis, opt, maxiter = 100, normalize = true)
 print_equations(Ψ)
 get_error(Ψ)
 ```

docs/src/quickstart.md

@@ -191,9 +191,9 @@ Looking at the eigenvalues of the system, we see that the estimated eigenvalues
 
 ## Nonlinear Systems - Sparse Identification of Nonlinear Dynamics
 
-Okay, so far we can fit linear models via DMD and nonlinear models via EDMD. But what if we want to find a model of a nonlinear system *without moving to Koopman space*? Simple, we use [Sparse Identification of Nonlinear Dynamics](https://www.pnas.org/content/113/15/3932) or `SInDy`.
+Okay, so far we can fit linear models via DMD and nonlinear models via EDMD. But what if we want to find a model of a nonlinear system *without moving to Koopman space*? Simple, we use [Sparse Identification of Nonlinear Dynamics](https://www.pnas.org/content/113/15/3932) or `SINDy`.
 
-As the name suggests, `SInDy` finds the sparsest basis of functions which build the observed trajectory. Again, we will start with a nonlinear system
+As the name suggests, `SINDy` finds the sparsest basis of functions which build the observed trajectory. Again, we will start with a nonlinear system
 
 ```@example 3
 using DataDrivenDiffEq
@@ -224,7 +224,7 @@ savefig("nonlinear_pendulum.png") # hide
 
 which is the simple nonlinear pendulum with damping.
 
-Suppose we are like John and know nothing about the system, we have just the data in front of us. To apply `SInDy`, we need three ingredients:
+Suppose we are like John and know nothing about the system, we have just the data in front of us. To apply `SINDy`, we need three ingredients:
 
 + A `Basis` containing all possible candidate functions which might be in the model
 + An optimizer which is able to produce a sparse output
@@ -248,11 +248,11 @@ nothing # hide
 
 ```@example 3
 opt = SR3(3e-1, 1.0)
-Ψ = SInDy(X[:, 1:1000], DX[:, 1:1000], basis, maxiter = 10000, opt = opt, normalize = true)
+Ψ = SINDy(X[:, 1:1000], DX[:, 1:1000], basis, opt, maxiter = 10000, normalize = true)
 print_equations(Ψ) # hide
 ```
 
-We recovered the equations! Let's transform the `SInDyResult` into a performant piece of
+We recovered the equations! Let's transform the `SINDyResult` into a performant piece of
 Julia Code using `ODESystem`
 
 ```@example 3
@@ -267,6 +267,6 @@ estimation = solve(estimator, Tsit5(), saveat = solution.t)
 plot(solution.t[1:1000], solution[:,1:1000]', color = :red, line = :dot, label = nothing) # hide
 plot!(solution.t[1000:end], solution[:,1000:end]', color = :blue, line = :dot,label = nothing) # hide
 plot!(estimation, color = :green, label = "Estimation") # hide
-savefig("sindy_estimation.png") # hide
+savefig("SINDy_estimation.png") # hide
 ```
-![](sindy_estimation.png)
+![](SINDy_estimation.png)

docs/src/sparse_identification/isindy.md

@@ -1,6 +1,6 @@
 # Implicit Sparse Identification of Nonlinear Dynamics
 
-While `SInDy` works well for ODEs, some systems take the form of [DAE](https://diffeq.sciml.ai/stable/types/dae_types/)s. A common form is `f(x, p, t) - g(x, p, t)*dx = 0`. These can be inferred via `ISInDy`, which extends `SInDy` [for Implicit problems](https://ieeexplore.ieee.org/abstract/document/7809160). In particular, it solves
+While `SINDy` works well for ODEs, some systems take the form of [DAE](https://diffeq.sciml.ai/stable/types/dae_types/)s. A common form is `f(x, p, t) - g(x, p, t)*dx = 0`. These can be inferred via `ISINDy`, which extends `SINDy` [for Implicit problems](https://ieeexplore.ieee.org/abstract/document/7809160). In particular, it solves
 
 ```math
 \Xi = \min ~ \left\lVert \Theta(X, p, t)^{T} \Xi \right\rVert_{2} + \lambda ~ \left\lVert \Xi \right\rVert_{1}
@@ -13,7 +13,7 @@ where ``\Xi`` lies in the nullspace of ``\Theta``.
 Let's try to infer the [Michaelis-Menten Kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics), like in the corresponding paper. We start by generating the
 corresponding data.
 
-```@example isindy_1
+```@example iSINDy_1
 using DataDrivenDiffEq
 using ModelingToolkit
 using OrdinaryDiffEq
@@ -32,11 +32,11 @@ problem = ODEProblem(michaelis_menten, u0, tspan)
 solution = solve(problem, Tsit5(), saveat = 0.1, atol = 1e-7, rtol = 1e-7)
     
 plot(solution) # hide
-savefig("isindy_example.png")
+savefig("iSINDy_example.png")
 ```
-![](isindy_example.png)
+![](iSINDy_example.png)
 
-```@example isindy_1
+```@example iSINDy_1
 X = solution[:,:]
 DX = similar(X)
 for (i, xi) in enumerate(eachcol(X))
@@ -47,24 +47,24 @@ end
 basis= Basis([u^i for i in 0:4], [u])
 ```
 
-The signature of `ISInDy` is equal to `SInDy`, but requires an `AbstractSubspaceOptimizer`. Currently, `DataDrivenDiffEq` just implements `ADM()` based on [alternating directions](https://arxiv.org/pdf/1412.4659.pdf). `rtol` gets passed into the derivation of the `nullspace` via `LinearAlgebra`.
+The signature of `ISINDy` is equal to `SINDy`, but requires an `AbstractSubspaceOptimizer`. Currently, `DataDrivenDiffEq` just implements `ADM()` based on [alternating directions](https://arxiv.org/pdf/1412.4659.pdf). `rtol` gets passed into the derivation of the `nullspace` via `LinearAlgebra`.
 
 
-```@example isindy_1
+```@example iSINDy_1
 opt = ADM(1.1e-1)
 ```
 
-Since `ADM()` returns sparsified columns of the nullspace we need to find a pareto optimal solution. To achieve this, we provide a sufficient cost function `g` to `ISInDy`. This allows us to evaluate each individual column of the sparse matrix on its 0-norm (sparsity) and the 2-norm of the matrix vector product of ``\Theta^T \xi`` (nullspace). This is a default setting which can be changed by providing a function `f` which maps the coefficients and the library onto a feature space. Here, we want to set the focus on the the magnitude of the deviation from the nullspace.
+Since `ADM()` returns sparsified columns of the nullspace we need to find a pareto optimal solution. To achieve this, we provide a sufficient cost function `g` to `ISINDy`. This allows us to evaluate each individual column of the sparse matrix on its 0-norm (sparsity) and the 2-norm of the matrix vector product of ``\Theta^T \xi`` (nullspace). This is a default setting which can be changed by providing a function `f` which maps the coefficients and the library onto a feature space. Here, we want to set the focus on the the magnitude of the deviation from the nullspace.
 
-```@example isindy_1
-Ψ = ISInDy(X, DX, basis, opt, g = x->norm([1e-1*x[1]; x[2]]), maxiter = 1000, rtol = 0.99)
+```@example iSINDy_1
+Ψ = ISINDy(X, DX, basis, opt, g = x->norm([1e-1*x[1]; x[2]]), maxiter = 1000, rtol = 0.99)
 nothing #hide
 ```
 
 The function call returns a `SparseIdentificationResult`.
 As in [Sparse Identification of Nonlinear Dynamics](@ref), we can transform the `SparseIdentificationResult` into an `ODESystem`.
 
-```@example isindy_1
+```@example iSINDy_1
 # Transform into ODE System
 sys = ODESystem(Ψ)
 dudt = ODEFunction(sys)
@@ -75,13 +75,13 @@ estimation = solve(estimator, Tsit5(), saveat = 0.1)
 
 plot(solution, color = :red, label = "True") # hide
 plot!(estimation, color = :green, label = "Estimation") # hide
-savefig("isindy_example_final.png") # hide
+savefig("iSINDy_example_final.png") # hide
 ```
-![](isindy_example_final.png)
+![](iSINDy_example_final.png)
 
-The model recovered by `ISInDy` is  correct
+The model recovered by `ISINDy` is  correct
 
-```@example isindy_1
+```@example iSINDy_1
 print_equations(Ψ)
 ```
 
@@ -92,7 +92,7 @@ The parameters are off a little, but, as before, we can use `DiffEqFlux` to tune
 
 Implicit dynamics can also be reformulated as an explicit problem as stated in [this paper](https://arxiv.org/pdf/2004.02322.pdf). The algorithm searches the correct equations by trying out all candidate functions as a right hand side and performing a sparse regression onto the remaining set of candidates. Let's start by defining the problem and generate the data:
 
-```@example isindy_2
+```@example iSINDy_2
 
 using DataDrivenDiffEq
 using ModelingToolkit
@@ -125,15 +125,15 @@ for (i, xi) in enumerate(eachcol(X))
 end
 
 plot(solution) # hide
-savefig("isindy_cartpole_data.png") # hide
+savefig("iSINDy_cartpole_data.png") # hide
 
 ```
-![](isindy_cartpole_data.png)
+![](iSINDy_cartpole_data.png)
 
 We see that we include a forcing term `F` inside the model which is depending on `t`.
 As before, we will also need a `Basis` to derive our equations from:
 
-```@example isindy_2
+```@example iSINDy_2
 @variables u[1:4] t
 polys = Operation[]
 for i ∈ 0:4
@@ -164,12 +164,12 @@ We added the time dependent input directly into the basis to account for its inf
 
 *NOTE : Including input signals may change in future releases!*
 
-Like for a `SInDy`, we can use any `AbstractOptimizer` with a pareto front optimization over different thresholds. 
+Like for a `SINDy`, we can use any `AbstractOptimizer` with a pareto front optimization over different thresholds. 
 
-```@example isindy_2
+```@example iSINDy_2
 λ = exp10.(-4:0.1:-1)
 g(x) = norm([1e-3; 10.0] .* x, 2)
-Ψ = ISInDy(X[:,:], DX[:, :], basis, λ, STRRidge(), maxiter = 100, normalize = false, t = solution.t, g = g)
+Ψ = ISINDy(X[:,:], DX[:, :], basis, λ, STRRidge(), maxiter = 100, normalize = false, t = solution.t, g = g)
 
 # Transform into ODE System
 sys = ODESystem(Ψ)
@@ -182,13 +182,13 @@ sol_ = solve(estimator, Tsit5(), saveat = dt)
 
 plot(solution.t[:], solution[:,:]', color = :red, label = nothing) # hide
 plot!(sol_.t, sol_[:, :]', color = :green, label = "Estimation") # hide
-savefig("isindy_cartpole_estimation.png") # hide
+savefig("iSINDy_cartpole_estimation.png") # hide
 ```
-![](isindy_cartpole_estimation.png)
+![](iSINDy_cartpole_estimation.png)
 
 Let's have a look at the equations recovered. They nearly match up.
 
-```@example isindy_2
+```@example iSINDy_2
 print_equations(Ψ)
 ```
 
@@ -198,5 +198,5 @@ print_equations(Ψ)
 ## Functions
 
 ```@docs
-ISInDy
+ISINDy
 ```

docs/src/sparse_identification/optimizers.md

@@ -24,7 +24,7 @@ mutable struct MyOpt <: DataDrivenDiffEq.Optimize.AbstractOptimizer
 end
 ```
 
-To use `MyOpt` within `SInDy`, an `init!` function has to be implemented.
+To use `MyOpt` within `SINDy`, an `init!` function has to be implemented.
 
 ```julia
 function init!(X::AbstractArray, o::MyOpt, A::AbstractArray, Y::AbstractArray)

docs/src/sparse_identification/sindy.md

@@ -15,7 +15,7 @@ where, in most cases, ``Y``is the data matrix containing the derivatives of the
 As in the original paper, we will estimate the [Lorenz System](https://en.wikipedia.org/wiki/Lorenz_system).
 First, let's create the necessary data and have a look at the trajectory.
 
-```@example sindy_1
+```@example SINDy_1
 using DataDrivenDiffEq
 using ModelingToolkit
 using OrdinaryDiffEq
@@ -47,7 +47,7 @@ savefig("lorenz.png") #hide
 
 Additionally, we generate the *ideal* derivative data.
 
-```@example sindy_1
+```@example SINDy_1
 X = Array(solution)
 DX = similar(X)
 for (i, xi) in enumerate(eachcol(X))
@@ -57,7 +57,7 @@ end
 
 To generate the symbolic equations, we need to define a ` Basis` over the variables `x y z`. In this example, we will use all monomials up to degree of 4 and their products:
 
-```@example sindy_1
+```@example SINDy_1
 @variables x y z
 u = Operation[x; y; z]
 polys = Operation[]
@@ -79,25 +79,25 @@ nothing #hide
 
 To perform the sparse identification on our data, we need to define an `Optimizer`. Here, we will use `STRRidge`, which is described in the original paper. The threshold of the optimizer is set to `0.1`. An overview of the different optimizers can be found below.
 
-```@example sindy_1
+```@example SINDy_1
 opt = STRRidge(0.1)
-Ψ = SInDy(X, DX, basis, maxiter = 100, opt = opt, normalize = true)
+Ψ = SINDy(X, DX, basis, opt, maxiter = 100, normalize = true)
 ```
 
-`Ψ` is a `SInDyResult`, which stores some about the regression. As we can see, we have 7 active terms inside the model.
+`Ψ` is a `SINDyResult`, which stores some about the regression. As we can see, we have 7 active terms inside the model.
 To look at the equations, simply type
 
-```@example sindy_1
+```@example SINDy_1
 print_equations(Ψ)
 ```
 
 First, let's have a look at the ``L2``-Error and Akaikes Information Criterion of the result
 
-```@example sindy_1
+```@example SINDy_1
 get_error(Ψ)
 ```
 
-```@example sindy_1
+```@example SINDy_1
 get_aicc(Ψ)
 ```
 
@@ -106,7 +106,7 @@ We can also access the coefficient matrix ``\Xi`` directly via `get_coefficients
 To generate a numerical model usable in `DifferentialEquations`, we simply use the `ODESystem` function from `ModelingToolkit`.
 The resulting parameters used for the identification can be accessed via `parameters(Ψ)`. The returned vector also includes the parameters of the original `Basis` used to generate the result.
 
-```@example sindy_1
+```@example SINDy_1
 ps = parameters(Ψ)
 sys = ODESystem(Ψ)
 dudt = ODEFunction(sys)
@@ -137,7 +137,7 @@ which resembles the papers results. Next, we could use [classical parameter esti
 ## Functions
 
 ```@docs
-SInDy
+SINDy
 sparse_regression
 ```
 

examples/Havok_Examples.jl

@@ -60,7 +60,7 @@ plot(plot(z[1, :]),plot(z[end, :].^2),  layout = (2, 1))
 
 basis = Basis(u, u)
 opt = SR3(1e-1)
-b = SInDy(z[:, 1:end], dz[1:end-1, 1:end], basis, maxiter = 1000, opt = opt, normalize = true)
+b = SINDy(z[:, 1:end], dz[1:end-1, 1:end], basis, opt, maxiter = 1000, normalize = true)
 
 # Coincides with the paper results
 println(b)

examples/ISInDy_Examples.jl

@@ -40,15 +40,15 @@ for i ∈ 0:6
 end
 
 basis= Basis(polys, u)
-Ψ = ISInDy(X, DX, basis, ADM(1e-2), maxiter = 10, rtol = 0.1)
+Ψ = ISINDy(X, DX, basis, ADM(1e-2), maxiter = 10, rtol = 0.1)
 println(Ψ)
 print_equations(Ψ)
 
-# Lets use sindy pi for a parallel implicit implementation
+# Lets use SINDy pi for a parallel implicit implementation
 # Create a basis 
 @variables u[1:3]
 polys = Operation[]
-# Lots of basis functions -> sindy pi can handle more than ADM()
+# Lots of basis functions -> SINDy pi can handle more than ADM()
 for i ∈ 0:10
     if i == 0
         push!(polys, u[1]^0)
@@ -62,8 +62,8 @@ end
 
 basis= Basis(polys, u)
 
-# Simply use any optimizer you would use for sindy
-Ψ = ISInDy(X[:, :], DX[:, :], basis, STRRidge(1e-1), maxiter = 100, normalize = true)
+# Simply use any optimizer you would use for SINDy
+Ψ = ISINDy(X[:, :], DX[:, :], basis, STRRidge(1e-1), maxiter = 100, normalize = true)
 
 
 # Transform into ODE System

examples/ISInDy_Examples2.jl

@@ -63,7 +63,7 @@ basis= Basis(polys, u, iv = t)
 # Simply use any optimizer you would use for sindy
 λ = exp10.(-4:0.1:-1)
 g(x) = norm([1e-3; 10.0] .* x, 2)
-Ψ = ISInDy(X[:,:], DX[:, :], basis, λ, STRRidge(), maxiter = 100, normalize = false, t = solution.t, g = g)
+Ψ = ISINDy(X[:,:], DX[:, :], basis, λ, STRRidge(), maxiter = 100, normalize = false, t = solution.t, g = g)
 println(Ψ)
 print_equations(Ψ, show_parameter = true)
 

examples/SInDy_Examples.jl

@@ -50,12 +50,13 @@ println(basis)
 # than assumptions, converges fast but fails sometimes with too much noise
 opt = STRRidge(1e-2)
 # Enforce all 100 iterations
-Ψ = SInDy(sol[:,1:25], DX[:, 1:25], basis, p = [1.0; 1.0], maxiter = 100, opt = opt, convergence_error = 1e-5)
+Ψ = SINDy(sol[:,1:25], DX[:, 1:25], basis, opt, p = [1.0; 1.0], maxiter = 100, convergence_error = 1e-5)
 println(Ψ)
+print_equations(Ψ)
 
 # Lasso as ADMM, typically needs more information, more tuning
 opt = ADMM(1e-2, 1.0)
-Ψ = SInDy(sol[:,1:50], DX[:, 1:50], basis, p = [1.0; 1.0], maxiter = 5000, opt = opt, convergence_error = 1e-3)
+Ψ = SINDy(sol[:,1:50], DX[:, 1:50], basis, opt, p = [1.0; 1.0], maxiter = 5000, convergence_error = 1e-3)
 println(Ψ)
 print_equations(Ψ)
 
@@ -64,7 +65,7 @@ parameters(Ψ)
 
 # SR3, works good with lesser data and tuning
 opt = SR3(1e-2, 1.0)
-Ψ = SInDy(sol[:,1:end], DX[:, 1:end], basis, p = [0.5; 0.5], maxiter = 5000, opt = opt, convergence_error = 1e-5)
+Ψ = SINDy(sol[:,1:end], DX[:, 1:end], basis, opt, p = [0.5; 0.5], maxiter = 5000, convergence_error = 1e-5)
 println(Ψ)
 print_equations(Ψ, show_parameter = true)
 
@@ -73,8 +74,9 @@ print_equations(Ψ, show_parameter = true)
 λs = exp10.(-5:0.1:-1)
 # Use SR3 with high relaxation (allows the solution to diverge from LTSQ) and high iterations
 opt = SR3(1e-2, 5.0)
-Ψ = SInDy(sol[:,1:10], DX[:, 1:10], basis, λs, p = [1.0; 1.0], maxiter = 15000, opt = opt)
+Ψ = SINDy(sol[:,1:10], DX[:, 1:10], basis, λs, opt, p = [1.0; 1.0], maxiter = 15000)
 println(Ψ)
+print_equations(Ψ)
 
 # Transform into ODE System
 sys = ODESystem(Ψ)

examples/SInDy_Examples2.jl

@@ -71,7 +71,7 @@ basis = Basis(h, u)
 
 # Get the reduced basis via the sparse regression
 opt = STRRidge(0.1)
-Ψ = SInDy(X_cropped, DX_sg, basis, maxiter = 100, opt = opt, normalize = false)
+Ψ = SINDy(X_cropped, DX_sg, basis, opt, maxiter = 100, normalize = false)
 print(Ψ)
 print_equations(Ψ)
 print_equations(Ψ, show_parameter = true)
@@ -84,5 +84,5 @@ X_noisy = X + 0.01*randn(seed,size(X))
 X_noisy_cropped, DX_sg = savitzky_golay(X_noisy, windowSize, polyOrder, deriv=1, dt=dt)
 X_noisy_cropped, X_smoothed = savitzky_golay(X_noisy, windowSize, polyOrder, deriv=0, dt=dt)
 
-Ψ = SInDy(X_smoothed, DX_sg, basis, maxiter = 100, opt = opt, normalize = false)
-print_equations(Ψ)
+Ψ = SINDy(X_smoothed, DX_sg, basis, opt, maxiter = 100, normalize = false)
+print_equations(Ψ, show_parameter = true)