Merge pull request #138 from SciML/AlCap23-patch-2

Update docs and release version

Project.toml

@@ -1,7 +1,7 @@
 name = "DataDrivenDiffEq"
 uuid = "2445eb08-9709-466a-b3fc-47e12bd697a2"
 authors = ["Julius Martensen <julius.martensen@gmail.com>"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

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 rational functions `dx = f(x) / g(x)`. 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}
@@ -90,7 +90,7 @@ The parameters are off a little, but, as before, we can use `DiffEqFlux` to tune
 
 ## Example : Cart-Pole with Time-Dependent Control
 
-Implicit dynamics can also be reformulated as an explicit problem as stated in [this paper](https://arxiv.org/pdf/2004.02322.pdf). The algorithm tries to find the right 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:
+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
 
@@ -130,7 +130,8 @@ savefig("isindy_cartpole_data.png") # hide
 ```
 ![](isindy_cartpole_data.png)
 
-As always, we will also need a `Basis` to derive our equations from:
+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
 @variables u[1:4] t
@@ -161,7 +162,9 @@ basis= Basis(polys, u, iv = t)
 
 We added the time dependent input directly into the basis to account for its influence.
 
-Like for a normal sparse regression, we can use any `AbstractOptimizer` with a pareto front optimization over different thresholds. 
+*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. 
 
 ```@example isindy_2
 λ = exp10.(-4:0.1:-1)
@@ -177,14 +180,13 @@ ps = parameters(Ψ)
 estimator = ODEProblem(dudt, u0, tspan, ps)
 sol_ = solve(estimator, Tsit5(), saveat = dt)
 
-# Yeah! We got it right
 plot(solution.t[:], solution[:,:]', color = :red, label = nothing) # hide
 plot!(sol_.t, sol_[:, :]', color = :green, label = "Estimation") # hide
 savefig("isindy_cartpole_estimation.png") # hide
 ```
 ![](isindy_cartpole_estimation.png)
 
-Let's have a look at the equations recovered.
+Let's have a look at the equations recovered. They nearly match up.
 
 ```@example isindy_2
 print_equations(Ψ)

src/sindy/isindy.jl

@@ -14,7 +14,7 @@
 Performs an implicit sparse identification of nonlinear dynamics given the data matrices `X` and `Y` via the `AbstractBasis` `basis.`
 Keyworded arguments include the parameter (values) of the basis `p` and the timepoints `t`, which are passed in optionally.
 Tries to find a sparse vector inside the nullspace if `opt` is an `AbstractSubspaceOptimizer` or performs parallel implicit search if `opt` is a `AbstractOptimizer`.
-`maxiter` the maximum iterations to perform and `convergence_error` the
+`maxiter` denote the maximum iterations to perform and `convergence_error` the
 bound which causes the optimizer to stop. `denoise` defines if the matrix holding candidate trajectories should be thresholded via the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
 `normalize` normalizes the matrix holding candidate trajectories via the L2-Norm over each function.
 
@@ -145,8 +145,6 @@ function ImplicitSparseIdentificationResult(coeff::AbstractArray, equations::Bas
     b_, p_ = derive_implicit_parameterized_eqs(coeff, equations)
     ps = [p; p_]
 
-    println(sparsities)
-    println(b_)
     Ŷ = b_(X, ps, t)
     training_error = norm.(eachrow(Y-Ŷ), 2)
     aicc = similar(training_error)

test/isindy.jl

@@ -2,59 +2,59 @@
 @info "Starting implicit SINDy tests"
 @testset "ISInDy" begin
 
-    #@info "Nonlinear Implicit System"
-    ## Create a test problem
-    #function simple(u, p, t)
-    #    return [(2.0u[2]^2 - 3.0)/(1.0 + u[1]^2); -u[1]^2/(2.0 + u[2]^2); (1-u[2])/(1+u[3]^2)]
-    #end
-#
-    #u0 = [2.37; 1.58; -3.10]
-    #tspan = (0.0, 10.0)
-    #prob = ODEProblem(simple, u0, tspan)
-    #sol = solve(prob, Tsit5(), saveat = 0.1)
-#
-    ## Create the differential data
-    #X = sol[:,:]
-    #DX = similar(X)
-    #for (i, xi) in enumerate(eachcol(X))
-    #    DX[:, i] = simple(xi, [], 0.0)
-    #end
-#
-    ## Create a basis
-    #@variables u[1:3]
-    #polys = ModelingToolkit.Operation[]
-    ## Lots of basis functions
-    #for i ∈ 0:5
-    #    if i == 0
-    #        push!(polys, u[1]^0)
-    #    end
-    #    for ui in u
-    #        if i > 0
-    #            push!(polys, ui^i)
-    #        end
-    #    end
-    #end
-#
-    #basis= Basis(polys, u)
-#
-    #opt = ADM(1e-2)
-    #Ψ = ISInDy(X, DX, basis, opt = opt, maxiter = 100)
-#
-    ## Transform into ODE System
-    #sys = ODESystem(Ψ)
-    #dudt = ODEFunction(sys)
-    #ps = parameters(Ψ)
-    #print_equations(Ψ, show_parameter = true)
-#
-    #@test all(get_error(Ψ) .< 1e-6)
-    #@test length(ps) == 11
-    #@test get_sparsity(Ψ) == [4; 3; 4]
-    #
-    ## Simulate
-    #estimator = ODEProblem(dudt, u0, tspan, ps)
-    #sol_ = solve(estimator, Tsit5(), saveat = 0.1)
-    #@test sol[:,:] ≈ sol_[:,:]
-    #@test abs.(ps) ≈ abs.(Float64[-1/3 ; -1/3 ; -1.00 ; 2/3; 1.00 ;0.5 ;0.5 ; 1.0; 1.0; -1.0; 1.0])
+    @info "Nonlinear Implicit System"
+    # Create a test problem
+    function simple(u, p, t)
+        return [(2.0u[2]^2 - 3.0)/(1.0 + u[1]^2); -u[1]^2/(2.0 + u[2]^2); (1-u[2])/(1+u[3]^2)]
+    end
+
+    u0 = [2.37; 1.58; -3.10]
+    tspan = (0.0, 10.0)
+    prob = ODEProblem(simple, u0, tspan)
+    sol = solve(prob, Tsit5(), saveat = 0.1)
+
+    # Create the differential data
+    X = sol[:,:]
+    DX = similar(X)
+    for (i, xi) in enumerate(eachcol(X))
+        DX[:, i] = simple(xi, [], 0.0)
+    end
+
+    # Create a basis
+    @variables u[1:3]
+    polys = ModelingToolkit.Operation[]
+    # Lots of basis functions
+    for i ∈ 0:5
+        if i == 0
+            push!(polys, u[1]^0)
+        end
+        for ui in u
+            if i > 0
+                push!(polys, ui^i)
+            end
+        end
+    end
+
+    basis= Basis(polys, u)
+
+    opt = ADM(1e-2)
+    Ψ = ISInDy(X, DX, basis, opt, maxiter = 100)
+
+    # Transform into ODE System
+    sys = ODESystem(Ψ)
+    dudt = ODEFunction(sys)
+    ps = parameters(Ψ)
+    print_equations(Ψ, show_parameter = true)
+
+    @test all(get_error(Ψ) .< 1e-6)
+    @test length(ps) == 11
+    @test get_sparsity(Ψ) == [4; 3; 4]
+
+    # Simulate
+    estimator = ODEProblem(dudt, u0, tspan, ps)
+    sol_ = solve(estimator, Tsit5(), saveat = 0.1)
+    @test sol[:,:] ≈ sol_[:,:]
+    @test abs.(ps) ≈ abs.(Float64[-1/3 ; -1/3 ; -1.00 ; 2/3; 1.00 ;0.5 ;0.5 ; 1.0; 1.0; -1.0; 1.0])
 
 
     @info "Michaelis-Menten-Kinetics"