Merge pull request #168 from SciML/MTK4_updates

MTK 4 Update

Author: AlCap23

Project.toml

@@ -8,6 +8,7 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
@@ -22,8 +23,9 @@ Compat = "2.2, 3.0"
 DSP = "0.6"
 DataInterpolations = "3.1"
 DiffEqBase = "6.45"
+DocStringExtensions = "0.7, 0.8"
 FiniteDifferences = "0.9.6, 0.10"
-ModelingToolkit = "3.17"
+ModelingToolkit = "4"
 ProximalOperators = "0.11, 0.12"
 QuadGK = "2.4"
 StatsBase = "0.32.0, 0.33"

README.md

@@ -54,8 +54,8 @@ end
 ################################################################
 
 @variables x y z
-u = Operation[x; y; z]
-polys = Operation[]
+u = [x; y; z]
+polys = Any[]
 for i ∈ 0:4
     for j ∈ 0:i
         for k ∈ 0:j

docs/src/basis.md

@@ -4,21 +4,23 @@ Many of the methods require the definition of a `Basis` on observables or
 functional forms. A `Basis` is generated via:
 
 ```julia
-Basis(h, u, parameters = [], iv = nothing)
+Basis(eqs::AbstractVector, states::AbstractVector; 
+    parameters::AbstractArray = [], iv = nothing,
+    simplify = false, linear_independent = false, name = gensym(:Basis), 
+    pins = [], observed = [], eval_expression = false,
+    kwargs...)
 ```
+where `eqs` is either a vector containing symbolic functions using 'ModelingToolkit.jl' or a general function with the typical DiffEq signature `h(u,p,t)`, which can be used with an `Num` or vector of `Num`. `states` are the dependent variables used to describe the Basis, and
+`parameters` are the optional parameters in the `Basis`. `iv` represents the independent variable of the system - in most cases the time. Additional arguments are `simplify`, which simplifies `eqs` before creating a `Basis`. `linear_dependent` breaks up `eqs` in linear independent elements which are unique. `name` is an optional name for the `Basis`, `pins` and `observed` can be using in accordance to ModelingToolkits documentation. `eval_expression` is used to generate a callable function from the eqs. If set to `false`, callable code will be returned. `true` will use `eval` on code returned from the function, which might cause worldage issues. 
 
-where `h` is either a vector of ModelingToolkit `Operation`s for the valid functional
-forms or a general function with the typical DiffEq signature `h(u,p,t)`, which can be used with an  `Operation` or vector of `Operation`. `u` are the ModelingToolkit `Variable`s used to describe the Basis, and
-`parameters` are the optional ModelingToolkit `Variable`s used to describe the
-parameters in the basis elements. `iv` represents the independent variable of the system - the time.
 
 ```@docs
 Basis
 ```
 
 ## Example
 
-We start by crearting some `Variables` and `Parameters` using `ModelingToolkit`.
+We start by crearting some variables and parameters using `ModelingToolkit`.
 ```@example basis
 using LinearAlgebra
 using DataDrivenDiffEq
@@ -29,14 +31,11 @@ using ModelingToolkit
 @parameters w[1:2]
 ```
 
-To define a basis, simply write down the equations you want to be included as a
-`Vector{Operation}`. Possible used parameters have to be given to the constructor.
+To define a basis, simply write down the equations you want to be included as a `Vector`. Possible used parameters have to be given to the constructor.
 ```@example basis
 h = [u[1]; u[2]; cos(w[1]*u[2]+w[2]*u[3])]
 b = Basis(h, u, parameters = w)
 ```
-What can a `Basis` do? Can it do stuff? Let's find out!
-
 `Basis` are callable with the signature of functions to be used in `DifferentialEquations`.
 So, the function value at a single point looks like:
 ```@example basis
@@ -70,7 +69,7 @@ push!(b, sin(u[1]))
 size(b)
 ```
 
-We can also define functions of time and add them
+We can also define functions of the independent variable and add them
 
 ```@example basis
 t = independent_variable(b)
@@ -116,11 +115,11 @@ b_f = Basis(f, u, parameters = w)
 println(b_f)
 ```
 
-This works for every function defined over `Operation`s. So to create a `Basis` from a `Flux` model, simply extend the activations used:
+This works for every function defined over `Num`s. So to create a `Basis` from a `Flux` model, simply extend the activations used:
 
 ```julia
 using Flux
-NNlib.σ(x::Operation) = 1 / (1+exp(-x))
+NNlib.σ(x::Num) = 1 / (1+exp(-x))
 
 c = Chain(Dense(3,2,σ), Dense(2, 1, σ))
 ps, re = Flux.destructure(c)
@@ -135,9 +134,6 @@ b = Basis(g, u, parameters = p)
 ## Functions
 
 ```@docs
-parameters
-variables
-DataDrivenDiffEq.independent_variable
 jacobian
 dynamics
 push!

docs/src/quickstart.md

@@ -238,7 +238,7 @@ So, let's create a bunch of basis functions for our problem first
 
 @variables u[1:2]
 
-h = Operation[u; u.^2; u.^3; sin.(u); cos.(u); 1]
+h = [u; u.^2; u.^3; sin.(u); cos.(u); 1]
 
 basis = Basis(h, u)
 nothing # hide

docs/src/sparse_identification/isindy.md

@@ -135,7 +135,7 @@ As before, we will also need a `Basis` to derive our equations from:
 
 ```@example iSINDy_2
 @variables u[1:4] t
-polys = Operation[]
+polys = Any[]
 for i ∈ 0:4
     if i == 0
         push!(polys, u[1]^0)
@@ -198,7 +198,7 @@ Alternatively, we can also use the input as an extended state `x`.
 
 ```@example iSINDy_2
 @variables u[1:4] t x
-polys = Operation[]
+polys = Any[]
 # Lots of basis functions -> sindy pi can handle more than ADM()
 for i ∈ 0:4
     if i == 0

docs/src/sparse_identification/sindy.md

@@ -59,8 +59,8 @@ To generate the symbolic equations, we need to define a ` Basis` over the variab
 
 ```@example SINDy_1
 @variables x y z
-u = Operation[x; y; z]
-polys = Operation[]
+u = [x; y; z]
+polys = Any[]
 for i ∈ 0:4
     for j ∈ 0:i
         for k ∈ 0:j

examples/Basis_Creation.jl

@@ -77,7 +77,7 @@ f(u, p, t) = [u[3]; u[2]*u[1]; sin(u[1])*u[2]]
 b = Basis(f, u)
 
 using Flux
-NNlib.σ(x::Operation) = 1 / (1+exp(-x))
+NNlib.σ(x::Num) = 1 / (1+exp(-x))
 # Build a fully
 c = Chain(Dense(3,5,σ), Dense(5, 2, σ))
 ps, re = Flux.destructure(c)

examples/DMDc_Examples.jl

@@ -21,4 +21,3 @@ eigvecs(sys)
 prob = DiscreteProblem(sys, X[:, 1], (0., 10.))
 sol_unforced = solve(prob,  FunctionMap())
 plot(sol_unforced)
-sol_unforced[:,:]

examples/ISInDy_Examples.jl

@@ -26,18 +26,7 @@ end
 
 # Create a basis
 @variables u[1:3]
-polys = Operation[]
-# Lots of basis functions
-for i ∈ 0:6
-    if i == 0
-        push!(polys, u[1]^0)
-    end
-    for ui in u
-        if i > 0
-            push!(polys, ui^i)
-        end
-    end
-end
+polys = [monomial_basis(u, 6);1]
 
 basis= Basis(polys, u)
 Ψ = ISINDy(X, DX, basis, ADM(1e-2), maxiter = 10, rtol = 0.1)
@@ -47,24 +36,14 @@ print_equations(Ψ)
 # 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()
-for i ∈ 0:10
-    if i == 0
-        push!(polys, u[1]^0)
-    end
-    for ui in u
-        if i > 0
-            push!(polys, ui^i)
-        end
-    end
-end
+polys = [monomial_basis(u, 10); 1]
 
 basis= Basis(polys, u)
 
 # Simply use any optimizer you would use for SINDy
 Ψ = ISINDy(X[:, :], DX[:, :], basis, STRRidge(1e-1), maxiter = 100, normalize = true)
-
+println(Ψ)
+print_equations(Ψ)
 
 # Transform into ODE System
 sys = ODESystem(Ψ)
@@ -80,4 +59,4 @@ plot(sol.t[:], sol[:,:]', color = :red, label = nothing)
 plot!(sol_.t, sol_[:, :]', color = :green, label = "Estimation")
 
 plot(sol.t, abs.(sol-sol_)')
-norm(sol[:,:]-sol_[:,:], 2) # approx 9e-7
+norm(sol[:,:]-sol_[:,:], 2) # approx 1.2e-13

examples/ISInDy_Examples2.jl

@@ -35,7 +35,7 @@ end
 plot(DX')
 
 @variables u[1:4] t
-polys = Operation[]
+polys = Any[]
 # Lots of basis functions -> sindy pi can handle more than ADM()
 for i ∈ 0:4
     if i == 0