Skip to content

Commit

Permalink
Rename SInDy->SINDy
Browse files Browse the repository at this point in the history
  • Loading branch information
AlCap23 committed Aug 19, 2020
1 parent 0c02978 commit 20a14b3
Show file tree
Hide file tree
Showing 16 changed files with 103 additions and 103 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ end
basis = Basis(polys, u)

opt = STRRidge(0.1)
Ψ = SInDy(X, DX, basis, opt, maxiter = 100, normalize = true)
Ψ = SINDy(X, DX, basis, opt, maxiter = 100, normalize = true)
print_equations(Ψ)
get_error(Ψ)
```
Expand Down
14 changes: 7 additions & 7 deletions docs/src/quickstart.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -248,11 +248,11 @@ nothing # hide

```@example 3
opt = SR3(3e-1, 1.0)
Ψ = SInDy(X[:, 1:1000], DX[:, 1:1000], basis, opt, maxiter = 10000, 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
Expand All @@ -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)
52 changes: 26 additions & 26 deletions docs/src/sparse_identification/isindy.md
Original file line number Diff line number Diff line change
@@ -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}
Expand All @@ -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
Expand All @@ -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))
Expand All @@ -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)
Expand All @@ -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(Ψ)
```

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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(Ψ)
Expand All @@ -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(Ψ)
```

Expand All @@ -198,5 +198,5 @@ print_equations(Ψ)
## Functions

```@docs
ISInDy
ISINDy
```
2 changes: 1 addition & 1 deletion docs/src/sparse_identification/optimizers.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
22 changes: 11 additions & 11 deletions docs/src/sparse_identification/sindy.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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[]
Expand All @@ -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, opt, maxiter = 100, 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(Ψ)
```

Expand All @@ -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)
Expand Down Expand Up @@ -137,7 +137,7 @@ which resembles the papers results. Next, we could use [classical parameter esti
## Functions

```@docs
SInDy
SINDy
sparse_regression
```

Expand Down
2 changes: 1 addition & 1 deletion examples/Havok_Examples.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
10 changes: 5 additions & 5 deletions examples/ISInDy_Examples.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion examples/ISInDy_Examples2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
8 changes: 4 additions & 4 deletions examples/SInDy_Examples.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,13 +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, opt, p = [1.0; 1.0], maxiter = 100, 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, opt, p = [1.0; 1.0], maxiter = 5000, 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(Ψ)

Expand All @@ -65,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, opt, p = [0.5; 0.5], maxiter = 5000, 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)

Expand All @@ -74,7 +74,7 @@ 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, opt, p = [1.0; 1.0], maxiter = 15000)
Ψ = SINDy(sol[:,1:10], DX[:, 1:10], basis, λs, opt, p = [1.0; 1.0], maxiter = 15000)
println(Ψ)
print_equations(Ψ)

Expand Down
4 changes: 2 additions & 2 deletions examples/SInDy_Examples2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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, opt, maxiter = 100, normalize = false)
Ψ = SINDy(X_cropped, DX_sg, basis, opt, maxiter = 100, normalize = false)
print(Ψ)
print_equations(Ψ)
print_equations(Ψ, show_parameter = true)
Expand All @@ -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, opt, maxiter = 100, normalize = false)
Ψ = SINDy(X_smoothed, DX_sg, basis, opt, maxiter = 100, normalize = false)
print_equations(Ψ, show_parameter = true)
4 changes: 2 additions & 2 deletions src/DataDrivenDiffEq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,11 +60,11 @@ export print_equations
export get_coefficients, get_error, get_sparsity, get_aicc

include("./sindy/sindy.jl")
export SInDy
export SINDy
export sparse_regression, sparse_regression!

include("./sindy/isindy.jl")
export ISInDy
export ISINDy

include("./utils.jl")
export AIC, AICC, BIC
Expand Down

0 comments on commit 20a14b3

Please sign in to comment.