diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml deleted file mode 100644 index 3e91a6c97..000000000 --- a/.JuliaFormatter.toml +++ /dev/null @@ -1,4 +0,0 @@ -style = "sciml" -format_markdown = true -format_docstrings = true -annotate_untyped_fields_with_any = false diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index c240796cc..6762c6f3e 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -1,13 +1,19 @@ -name: "Format Check" +name: format-check on: push: branches: - 'master' + - 'main' + - 'release-' tags: '*' pull_request: jobs: - format-check: - name: "Format Check" - uses: "SciML/.github/.github/workflows/format-check.yml@v1" + runic: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: fredrikekre/runic-action@v1 + with: + version: '1' diff --git a/docs/make.jl b/docs/make.jl index 844bf8887..14c1c7dfd 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,7 +5,7 @@ using Pkg function dev_subpkg(subpkg) subpkg_path = abspath(joinpath(dirname(@__FILE__), "..", "lib", subpkg)) - Pkg.develop(PackageSpec(path = subpkg_path)) + return Pkg.develop(PackageSpec(path = subpkg_path)) end dev_subpkg("DataDrivenDMD") @@ -50,25 +50,34 @@ function create_tutorials(dirname, targetdir, excludes = []) mdpost(str) = replace(str, "@__CODE__" => code) Literate.markdown(ipath, targetdir) Literate.markdown(ipath, targetdir, execute = false, postprocess = mdpost) - push!(tutorials, - relpath(joinpath(targetdir, fname * ".md"), joinpath(@__DIR__, "src"))) + push!( + tutorials, + relpath(joinpath(targetdir, fname * ".md"), joinpath(@__DIR__, "src")) + ) end end return tutorials end -koopman_tutorial = create_tutorials(joinpath(@__DIR__, "src/libs/datadrivendmd/"), - joinpath(@__DIR__, "src/libs/datadrivendmd/examples")) -sparse_tutorial = create_tutorials(joinpath(@__DIR__, "src/libs/datadrivensparse/"), - joinpath(@__DIR__, "src/libs/datadrivensparse/examples")) -sr_tutorial = create_tutorials(joinpath(@__DIR__, "src/libs/datadrivensr/"), - joinpath(@__DIR__, "src/libs/datadrivensr/examples")) +koopman_tutorial = create_tutorials( + joinpath(@__DIR__, "src/libs/datadrivendmd/"), + joinpath(@__DIR__, "src/libs/datadrivendmd/examples") +) +sparse_tutorial = create_tutorials( + joinpath(@__DIR__, "src/libs/datadrivensparse/"), + joinpath(@__DIR__, "src/libs/datadrivensparse/examples") +) +sr_tutorial = create_tutorials( + joinpath(@__DIR__, "src/libs/datadrivensr/"), + joinpath(@__DIR__, "src/libs/datadrivensr/examples") +) # Must be after tutorials is created include("pages.jl") # Create the docs -makedocs(sitename = "DataDrivenDiffEq.jl", +makedocs( + sitename = "DataDrivenDiffEq.jl", authors = "Julius Martensen, Christopher Rackauckas, et al.", modules = [DataDrivenDiffEq, DataDrivenDMD, DataDrivenSparse, DataDrivenSR], clean = true, doctest = false, linkcheck = true, @@ -77,11 +86,16 @@ makedocs(sitename = "DataDrivenDiffEq.jl", "http://cwrowley.princeton.edu/papers/Hemati-2017a.pdf", "https://royalsocietypublishing.org/doi/10.1098/rspa.2020.0279", "https://www.pnas.org/doi/10.1073/pnas.1517384113", - "https://link.springer.com/article/10.1007/s00332-015-9258-5" + "https://link.springer.com/article/10.1007/s00332-015-9258-5", ], - format = Documenter.HTML(assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/DataDrivenDiffEq/stable/"), - pages = pages) + format = Documenter.HTML( + assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/DataDrivenDiffEq/stable/" + ), + pages = pages +) -deploydocs(repo = "github.com/SciML/DataDrivenDiffEq.jl.git"; - push_preview = true) +deploydocs( + repo = "github.com/SciML/DataDrivenDiffEq.jl.git"; + push_preview = true +) diff --git a/docs/pages.jl b/docs/pages.jl index db079454a..0bceea5fb 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -5,14 +5,17 @@ pages = [ "Basis" => "basis.md", "Solutions" => "solutions.md", "Utilities" => "utils.md", - "DataDrivenDMD" => ["Background" => "libs/datadrivendmd/koopman.md", - "Examples" => koopman_tutorial + "DataDrivenDMD" => [ + "Background" => "libs/datadrivendmd/koopman.md", + "Examples" => koopman_tutorial, ], - "DataDrivenSparse" => ["Background" => "libs/datadrivensparse/sparse_regression.md", - "Examples" => sparse_tutorial + "DataDrivenSparse" => [ + "Background" => "libs/datadrivensparse/sparse_regression.md", + "Examples" => sparse_tutorial, ], - "DataDrivenSR" => ["Background" => "libs/datadrivensr/symbolic_regression.md", - "Examples" => sr_tutorial + "DataDrivenSR" => [ + "Background" => "libs/datadrivensr/symbolic_regression.md", + "Examples" => sr_tutorial, ], - "Citing" => "citations.md" + "Citing" => "citations.md", ] diff --git a/docs/src/libs/datadrivendmd/example_01.jl b/docs/src/libs/datadrivendmd/example_01.jl index ee0c72be5..1f09ba2b1 100644 --- a/docs/src/libs/datadrivendmd/example_01.jl +++ b/docs/src/libs/datadrivendmd/example_01.jl @@ -47,6 +47,6 @@ get_basis(res) #md # ``` ## Test the result #src -@test rss(res) <= 1e-3 #src +@test rss(res) <= 1.0e-3 #src @test r2(res) >= 0.99 #src @test dof(res) == 3 #src diff --git a/docs/src/libs/datadrivendmd/example_04.jl b/docs/src/libs/datadrivendmd/example_04.jl index 154dec650..dbe8d0982 100644 --- a/docs/src/libs/datadrivendmd/example_04.jl +++ b/docs/src/libs/datadrivendmd/example_04.jl @@ -10,7 +10,7 @@ using DataDrivenDMD function slow_manifold(du, u, p, t) du[1] = p[1] * u[1] - du[2] = p[2] * (u[2] - u[1]^2) + return du[2] = p[2] * (u[2] - u[1]^2) end u0 = [3.0; -2.0] diff --git a/docs/src/libs/datadrivendmd/example_05.jl b/docs/src/libs/datadrivendmd/example_05.jl index 19e4bc660..3190cbae1 100644 --- a/docs/src/libs/datadrivendmd/example_05.jl +++ b/docs/src/libs/datadrivendmd/example_05.jl @@ -41,7 +41,7 @@ end # System dynamics function controlled_system!(du, u, p, t) ctrl = control_signal(t) - du .= A_true * u .+ B_true .* ctrl + return du .= A_true * u .+ B_true .* ctrl end # Generate "experimental" data diff --git a/docs/src/libs/datadrivensparse/example_01.jl b/docs/src/libs/datadrivensparse/example_01.jl index 40a136550..3c364a000 100644 --- a/docs/src/libs/datadrivensparse/example_01.jl +++ b/docs/src/libs/datadrivensparse/example_01.jl @@ -1,6 +1,6 @@ # # [Getting Started](@id getting_started) # -# The workflow for [DataDrivenDiffEq.jl](https://github.com/SciML/DataDrivenDiffEq.jl) is similar to other [SciML](https://sciml.ai/) packages. +# The workflow for [DataDrivenDiffEq.jl](https://github.com/SciML/DataDrivenDiffEq.jl) is similar to other [SciML](https://sciml.ai/) packages. # You start by defining a [`DataDrivenProblem`](@ref) and then dispatch on the [`solve`](@ref solve) command to return a [`DataDrivenSolution`](@ref). # Here is an outline of the required elements and choices: @@ -13,7 +13,7 @@ using ModelingToolkit using LinearAlgebra using DataDrivenSparse -# Generate a test problem +# Generate a test problem f(u) = u .^ 2 .+ 2.0u .- 1.0 X = randn(1, 100); diff --git a/docs/src/libs/datadrivensparse/example_02.jl b/docs/src/libs/datadrivensparse/example_02.jl index 838692cd7..ceed53f5b 100644 --- a/docs/src/libs/datadrivensparse/example_02.jl +++ b/docs/src/libs/datadrivensparse/example_02.jl @@ -38,9 +38,11 @@ ts = sol.t; # Using a [collocation method](@ref collocation), it automatically provides the derivative and smoothes the trajectory. Control signals can be passed # in as a function `(u,p,t)->control` or an array of measurements. -prob = ContinuousDataDrivenProblem(X, ts, GaussianKernel(), +prob = ContinuousDataDrivenProblem( + X, ts, GaussianKernel(), U = (u, p, t) -> [exp(-((t - 5.0) / 5.0)^2)], - p = ones(2)) + p = ones(2) +) #md plot(prob, size = (600,600)) @@ -65,8 +67,10 @@ println(basis) # hide sampler = DataProcessing(split = 0.8, shuffle = true, batchsize = 30, rng = rng) λs = exp10.(-10:0.1:0) opt = STLSQ(λs) -res = solve(prob, basis, opt, - options = DataDrivenCommonOptions(data_processing = sampler, digits = 1)) +res = solve( + prob, basis, opt, + options = DataDrivenCommonOptions(data_processing = sampler, digits = 1) +) #src println(res) #hide # !!! info diff --git a/docs/src/libs/datadrivensparse/example_03.jl b/docs/src/libs/datadrivensparse/example_03.jl index 027bf010a..c9cd9c721 100644 --- a/docs/src/libs/datadrivensparse/example_03.jl +++ b/docs/src/libs/datadrivensparse/example_03.jl @@ -13,7 +13,7 @@ using DataDrivenSparse using Test #src function michaelis_menten(u, p, t) - [0.6 - 1.5u[1] / (0.3 + u[1])] + return [0.6 - 1.5u[1] / (0.3 + u[1])] end u0 = [0.5] @@ -23,10 +23,14 @@ ode_problem = ODEProblem(michaelis_menten, u0, (0.0, 4.0)); # Since we have multiple trajectories at hand, we define a [`DataDrivenDataset`](@ref), which collects multiple problems but handles them as a unit # for the processing. -prob = DataDrivenDataset(map(1:2) do i - solve(remake(ode_problem, u0 = i * u0), - Tsit5(), saveat = 0.1, tspan = (0.0, 4.0)) -end...) +prob = DataDrivenDataset( + map(1:2) do i + solve( + remake(ode_problem, u0 = i * u0), + Tsit5(), saveat = 0.1, tspan = (0.0, 4.0) + ) + end... +) #md plot(prob) @@ -44,7 +48,7 @@ basis = Basis([h; h .* (D(u[1]))], u, implicits = D.(u), iv = t) # Next, we define the [`ImplicitOptimizer`](@ref) and `solve` the problem. It wraps a standard optimizer, by default [`STLSQ`](@ref), and performs # implicit sparse regression upon the selected basis. -opt = ImplicitOptimizer(1e-1:1e-1:5e-1) +opt = ImplicitOptimizer(1.0e-1:1.0e-1:5.0e-1) res = solve(prob, basis, opt) #md println(res) #hide diff --git a/docs/src/libs/datadrivensparse/example_04.jl b/docs/src/libs/datadrivensparse/example_04.jl index 34e80ec70..d11af3ab8 100644 --- a/docs/src/libs/datadrivensparse/example_04.jl +++ b/docs/src/libs/datadrivensparse/example_04.jl @@ -49,8 +49,10 @@ basis = Basis(eqs, x, independent_variable = t, implicits = D.(x)) # Next to varying over different sparsity penalties, we also want to batch our data using [`DataProcessing`](@ref). sampler = DataProcessing(split = 0.8, shuffle = true, batchsize = 30) -res = solve(dd_prob, basis, ImplicitOptimizer(STLSQ(1e-2:1e-2:1.0)), - options = DataDrivenCommonOptions(data_processing = sampler, digits = 2)) +res = solve( + dd_prob, basis, ImplicitOptimizer(STLSQ(1.0e-2:1.0e-2:1.0)), + options = DataDrivenCommonOptions(data_processing = sampler, digits = 2) +) #md println(res) #hide # And have a look at the resulting plot diff --git a/docs/src/libs/datadrivensparse/example_05.jl b/docs/src/libs/datadrivensparse/example_05.jl index 005360655..7c9228738 100644 --- a/docs/src/libs/datadrivensparse/example_05.jl +++ b/docs/src/libs/datadrivensparse/example_05.jl @@ -17,7 +17,7 @@ function cart_pole(u, p, t) du[1] = u[3] du[2] = u[4] du[3] = -(19.62 * sin(u[1]) + sin(u[1]) * cos(u[1]) * u[3]^2 + F * cos(u[1])) / - (2 - cos(u[1])^2) + (2 - cos(u[1])^2) du[4] = -(sin(u[1]) * u[3]^2 + 9.81 * sin(u[1]) * cos(u[1]) + F) / (2 - cos(u[1])^2) return du end @@ -35,15 +35,17 @@ for (i, xi) in enumerate(eachcol(X)) end t = solution.t -ddprob = ContinuousDataDrivenProblem(X, t, DX = DX[3:4, :], - U = (u, p, t) -> [-0.2 + 0.5 * sin(6 * t)]) +ddprob = ContinuousDataDrivenProblem( + X, t, DX = DX[3:4, :], + U = (u, p, t) -> [-0.2 + 0.5 * sin(6 * t)] +) #md plot(ddprob) -# Note that we just included the third and forth time derivative, assuming that we already know that the velocity `x[3:4]` is equal to the time +# Note that we just included the third and forth time derivative, assuming that we already know that the velocity `x[3:4]` is equal to the time # derivative of the position `x[1:2]`. -# Next, we define a sufficient [`Basis`](@ref). Again, we need to include `implicits` in the definition of -# our candidate functions and inform the [`Basis`](@ref) of it. +# Next, we define a sufficient [`Basis`](@ref). Again, we need to include `implicits` in the definition of +# our candidate functions and inform the [`Basis`](@ref) of it. @parameters t @variables u[1:4] du[1:2] x[1:1] @@ -70,8 +72,10 @@ basis = Basis(implicits, u, implicits = du, controls = x, iv = t); # We solve the problem by varying over a sufficient set of thresholds for the associated optimizer. -λ = [1e-4; 5e-4; 1e-3; 2e-3; 3e-3; 4e-3; 5e-3; 6e-3; 7e-3; 8e-3; 9e-3; 1e-2; 2e-2; 3e-2; - 4e-2; 5e-2] +λ = [ + 1.0e-4; 5.0e-4; 1.0e-3; 2.0e-3; 3.0e-3; 4.0e-3; 5.0e-3; 6.0e-3; 7.0e-3; 8.0e-3; 9.0e-3; 1.0e-2; 2.0e-2; 3.0e-2; + 4.0e-2; 5.0e-2 +] opt = ImplicitOptimizer(λ) res = solve(ddprob, basis, opt) @@ -81,8 +85,8 @@ res = solve(ddprob, basis, opt) system = get_basis(res) #md println(system) # hide -# We have recovered the correct equations of motion! -# Another visual check using the problem and the result yields +# We have recovered the correct equations of motion! +# Another visual check using the problem and the result yields #md plot( #md plot(ddprob), plot(res), layout = (1,2) diff --git a/docs/src/libs/datadrivensr/example_01.jl b/docs/src/libs/datadrivensr/example_01.jl index 310fd0e10..d0091303b 100644 --- a/docs/src/libs/datadrivensr/example_01.jl +++ b/docs/src/libs/datadrivensr/example_01.jl @@ -36,10 +36,12 @@ prob = ContinuousDataDrivenProblem(X, t, U = U) # To solve our problem, we will use [`EQSearch`](@ref), which provides a wrapper for the [symbolic regression interface](https://ai.damtp.cam.ac.uk/symbolicregression/stable/api/#Options). # We will stick to simple operations, use a `L1DistLoss`, and limit the verbosity of the algorithm. -eqsearch_options = SymbolicRegression.Options(binary_operators = [+, *], +eqsearch_options = SymbolicRegression.Options( + binary_operators = [+, *], loss = L1DistLoss(), verbosity = -1, progress = false, npop = 30, - timeout_in_seconds = 60.0) + timeout_in_seconds = 60.0 +) alg = EQSearch(eq_options = eqsearch_options) @@ -68,5 +70,5 @@ system = get_basis(res) #md # ``` ## Test #src -@test rss(res) .<= 5e-1 #src +@test rss(res) .<= 5.0e-1 #src @test r2(res) >= 0.95 #src diff --git a/docs/src/libs/datadrivensr/example_02.jl b/docs/src/libs/datadrivensr/example_02.jl index e190da11e..daf1c59c8 100644 --- a/docs/src/libs/datadrivensr/example_02.jl +++ b/docs/src/libs/datadrivensr/example_02.jl @@ -11,7 +11,7 @@ using DataDrivenSR function pendulum!(du, u, p, t) du[1] = u[2] - du[2] = -9.81 * sin(u[1]) + return du[2] = -9.81 * sin(u[1]) end u0 = [0.1, π / 2] @@ -35,10 +35,12 @@ u = collect(u) basis = Basis([polynomial_basis(u, 2); sin.(u)], u) -eqsearch_options = SymbolicRegression.Options(binary_operators = [+, *], +eqsearch_options = SymbolicRegression.Options( + binary_operators = [+, *], loss = L1DistLoss(), verbosity = -1, progress = false, npop = 30, - timeout_in_seconds = 60.0) + timeout_in_seconds = 60.0 +) alg = EQSearch(eq_options = eqsearch_options) @@ -63,5 +65,5 @@ system = get_basis(res) #md # ``` ## Test #src -@test rss(res) .<= 5e-2 #src +@test rss(res) .<= 5.0e-2 #src @test r2(res) >= 0.95 #src diff --git a/lib/DataDrivenDMD/src/DataDrivenDMD.jl b/lib/DataDrivenDMD/src/DataDrivenDMD.jl index 99c97c6b6..43af19fb3 100644 --- a/lib/DataDrivenDMD/src/DataDrivenDMD.jl +++ b/lib/DataDrivenDMD/src/DataDrivenDMD.jl @@ -20,7 +20,7 @@ using LinearAlgebra abstract type AbstractKoopmanAlgorithm <: AbstractDataDrivenAlgorithm end -# Results +# Results include("./result.jl") export KoopmanResult export get_operator, get_inputmap, get_outputmap, get_trainerror, get_testerror diff --git a/lib/DataDrivenDMD/src/algorithms.jl b/lib/DataDrivenDMD/src/algorithms.jl index a4390def7..035023ed9 100644 --- a/lib/DataDrivenDMD/src/algorithms.jl +++ b/lib/DataDrivenDMD/src/algorithms.jl @@ -22,14 +22,16 @@ end # General method with inputs function (x::AbstractKoopmanAlgorithm)( X::AbstractArray, Y::AbstractArray, U::AbstractArray, - B::AbstractArray) + B::AbstractArray + ) K, _ = x(X, Y - B * U) return (K, B) end function (x::AbstractKoopmanAlgorithm)( X::AbstractArray, Y::AbstractArray, U::AbstractArray, - ::Nothing) + ::Nothing + ) return x(X, Y, U) end @@ -117,8 +119,10 @@ function (x::DMDSVD{T})(X::AbstractArray, Y::AbstractArray) where {T <: Real} end # DMDc -function (x::DMDSVD{T})(X::AbstractArray, Y::AbstractArray, - U::AbstractArray) where {T <: Real} +function (x::DMDSVD{T})( + X::AbstractArray, Y::AbstractArray, + U::AbstractArray + ) where {T <: Real} isempty(U) && return x(X, Y) nx, m = size(X) nu, m = size(U) @@ -164,7 +168,7 @@ $(FIELDS) $(SIGNATURES) """ mutable struct TOTALDMD{R, A} <: - AbstractKoopmanAlgorithm where {R <: Number, A <: AbstractKoopmanAlgorithm} + AbstractKoopmanAlgorithm where {R <: Number, A <: AbstractKoopmanAlgorithm} truncation::R alg::A end @@ -182,8 +186,10 @@ function (x::TOTALDMD)(X::AbstractArray, Y::AbstractArray, U::AbstractArray) return x.alg(X * Q, Y * Q, U * Q) end -function (x::TOTALDMD)(X::AbstractArray, Y::AbstractArray, U::AbstractArray, - B::AbstractArray) +function (x::TOTALDMD)( + X::AbstractArray, Y::AbstractArray, U::AbstractArray, + B::AbstractArray + ) _, _, Q = truncated_svd([X; Y], x.truncation) K, _ = x.alg(X * Q, (Y - B * U) * Q) return (K, B) diff --git a/lib/DataDrivenDMD/src/result.jl b/lib/DataDrivenDMD/src/result.jl index 993537d43..c2b8e4499 100644 --- a/lib/DataDrivenDMD/src/result.jl +++ b/lib/DataDrivenDMD/src/result.jl @@ -24,8 +24,10 @@ struct KoopmanResult{K, B, C, Q, P, T} <: AbstractDataDrivenResult """Returncode""" retcode::DDReturnCode - function KoopmanResult(k_::K, b::B, c::C, q::Q, p::P, X::AbstractMatrix{T}, - Y::AbstractMatrix{T}, U::AbstractMatrix) where {K, B, C, Q, P, T} + function KoopmanResult( + k_::K, b::B, c::C, q::Q, p::P, X::AbstractMatrix{T}, + Y::AbstractMatrix{T}, U::AbstractMatrix + ) where {K, B, C, Q, P, T} k = Matrix(k_) rss = isempty(b) ? sum(abs2, Y .- c * k * X) : sum(abs2, Y .- c * (k * X .+ b * U)) dof = sum(!iszero, k) @@ -34,7 +36,7 @@ struct KoopmanResult{K, B, C, Q, P, T} <: AbstractDataDrivenResult ll = -nobs / 2 * log(rss / nobs) nll = -nobs / 2 * log(mean(abs2, Y .- vec(mean(Y, dims = 2)))) - new{K, B, C, Q, P, T}(k_, b, c, q, p, rss, ll, nll, dof, nobs, DDReturnCode(1)) + return new{K, B, C, Q, P, T}(k_, b, c, q, p, rss, ll, nll, dof, nobs, DDReturnCode(1)) end end diff --git a/lib/DataDrivenDMD/src/solve.jl b/lib/DataDrivenDMD/src/solve.jl index e5308524d..81dc4f954 100644 --- a/lib/DataDrivenDMD/src/solve.jl +++ b/lib/DataDrivenDMD/src/solve.jl @@ -1,13 +1,15 @@ # This will get called within init in DataDrivenDiffEq -function DataDrivenDiffEq.get_fit_targets(::A, prob::ABSTRACT_CONT_PROB, - basis::AbstractBasis) where { +function DataDrivenDiffEq.get_fit_targets( + ::A, prob::ABSTRACT_CONT_PROB, + basis::AbstractBasis + ) where { A <: - AbstractKoopmanAlgorithm -} + AbstractKoopmanAlgorithm, + } @unpack DX, X, p, t, U = prob - @assert size(DX, 1)==size(X, 1) "$(A) needs equal number of observed states and differentials for continuous problems!" + @assert size(DX, 1) == size(X, 1) "$(A) needs equal number of observed states and differentials for continuous problems!" Θ = basis(prob) n_x = size(X, 1) @@ -31,24 +33,28 @@ function DataDrivenDiffEq.get_fit_targets(::A, prob::ABSTRACT_CONT_PROB, return Θ, Ỹ, DX end -function DataDrivenDiffEq.get_fit_targets(::A, prob::ABSTRACT_DISCRETE_PROB, - basis::AbstractBasis) where { +function DataDrivenDiffEq.get_fit_targets( + ::A, prob::ABSTRACT_DISCRETE_PROB, + basis::AbstractBasis + ) where { A <: - AbstractKoopmanAlgorithm -} + AbstractKoopmanAlgorithm, + } # TODO Maybe we could, but this would require X[:, i+2] -> split in three here @assert !is_implicit(basis) "$(A) does not support implicit arguments in the basis for discrete problems!" @unpack X, p, t, U = prob - # Lift + # Lift Θ = basis(prob) n_b, m = size(Θ) Ỹ = zeros(eltype(Θ), n_b, m) if is_controlled(basis) foreach(1:m) do i - Ỹ[:, i] .= basis(X[:, i + 1], p, t[i + 1], - U[:, i + 1]) + Ỹ[:, i] .= basis( + X[:, i + 1], p, t[i + 1], + U[:, i + 1] + ) end else foreach(1:m) do i @@ -58,14 +64,14 @@ function DataDrivenDiffEq.get_fit_targets(::A, prob::ABSTRACT_DISCRETE_PROB, return Θ, Ỹ, X[:, 2:end] end -## Solve the Koopman +## Solve the Koopman function CommonSolve.solve!(prob::InternalDataDrivenProblem{A}) where { A <: - AbstractKoopmanAlgorithm -} + AbstractKoopmanAlgorithm, + } @unpack alg, basis, testdata, traindata, control_idx, options, problem, kwargs = prob @unpack selector = options - # Check for + # Check for results = alg(prob; kwargs...) # Get the best result based on selector @@ -74,7 +80,7 @@ function CommonSolve.solve!(prob::InternalDataDrivenProblem{A}) where { # Convert to basis new_basis = convert_to_basis(best_res, basis, problem, options, control_idx) # Build DataDrivenResult - DataDrivenSolution(new_basis, problem, alg, results, prob, best_res.retcode) + return DataDrivenSolution(new_basis, problem, alg, results, prob, best_res.retcode) end function convert_to_basis(res::KoopmanResult, basis::Basis, prob, options, control_idx) @@ -90,11 +96,13 @@ function convert_to_basis(res::KoopmanResult, basis::Basis, prob, options, contr Θ .= c * Matrix(k) end - DataDrivenDiffEq.__construct_basis(Θ, basis, prob, options) + return DataDrivenDiffEq.__construct_basis(Θ, basis, prob, options) end -function (algorithm::AbstractKoopmanAlgorithm)(prob::InternalDataDrivenProblem; - control_input = nothing, kwargs...) +function (algorithm::AbstractKoopmanAlgorithm)( + prob::InternalDataDrivenProblem; + control_input = nothing, kwargs... + ) (; traindata, testdata, control_idx, options) = prob (; abstol) = options # Preprocess control idx, indicates if any control is active in a single basis atom @@ -109,7 +117,7 @@ function (algorithm::AbstractKoopmanAlgorithm)(prob::InternalDataDrivenProblem; X̃, Ũ = X̃, DataDrivenDiffEq.__EMPTY_MATRIX end - map(traindata) do (X, Y, Z) + return map(traindata) do (X, Y, Z) if any(control_idx) X_, Y_, U_ = X[no_controls, :], Y[no_controls, :], X[control_idx, :] else diff --git a/lib/DataDrivenDMD/src/type.jl b/lib/DataDrivenDMD/src/type.jl index 94c70d850..cbaded562 100644 --- a/lib/DataDrivenDMD/src/type.jl +++ b/lib/DataDrivenDMD/src/type.jl @@ -149,11 +149,13 @@ the L2 error of the prediction exceeds the `threshold`. `p` and `t` are the parameters of the basis and the vector of timepoints, if necessary. """ -function update!(k::AbstractKoopman, +function update!( + k::AbstractKoopman, X::AbstractArray, Y::AbstractArray; p::AbstractArray = [], t::AbstractVector = [], U::AbstractArray = [], - threshold::T = eps()) where {T <: Real} + threshold::T = eps() + ) where {T <: Real} @assert updatable(k) "Linear Koopman is not updatable." Ψ₀ = k(X, p, t, U) diff --git a/lib/DataDrivenDMD/test/linear_autonomous.jl b/lib/DataDrivenDMD/test/linear_autonomous.jl index b303d1bec..8bfc0c711 100644 --- a/lib/DataDrivenDMD/test/linear_autonomous.jl +++ b/lib/DataDrivenDMD/test/linear_autonomous.jl @@ -22,8 +22,8 @@ rng = StableRNG(42) for alg in [DMDPINV(), DMDSVD(), TOTALDMD()] res = solve(prob, alg) - @test rss(res) <= 1e-2 - @test r2(res)≈0.95 atol=5e-1 + @test rss(res) <= 1.0e-2 + @test r2(res) ≈ 0.95 atol = 5.0e-1 @test dof(res) == 3 @test loglikelihood(res) >= 400.0 @@ -31,7 +31,7 @@ rng = StableRNG(42) @test Matrix(get_operator(operator_res)) ≈ A @test isempty(get_inputmap(operator_res)) @test get_outputmap(operator_res) ≈ I(2) - @test rss(operator_res) <= 1e-10 + @test rss(operator_res) <= 1.0e-10 end end end @@ -42,16 +42,16 @@ rng = StableRNG(42) for alg in [DMDPINV(), DMDSVD(), TOTALDMD()] res = solve(prob, alg) - @test rss(res) <= 1e-2 - @test r2(res)≈1.0 atol=1e-1 + @test rss(res) <= 1.0e-2 + @test r2(res) ≈ 1.0 atol = 1.0e-1 @test dof(res) == 4 @test loglikelihood(res) >= 85.0 foreach(get_results(res)) do operator_res - @test Matrix(get_operator(operator_res))≈A atol=1e-2 + @test Matrix(get_operator(operator_res)) ≈ A atol = 1.0e-2 @test isempty(get_inputmap(operator_res)) @test get_outputmap(operator_res) ≈ I(2) - @test rss(operator_res) <= 1e-2 + @test rss(operator_res) <= 1.0e-2 end end end @@ -69,16 +69,16 @@ end for alg in [DMDPINV(), DMDSVD(), TOTALDMD()] res = solve(prob, alg) - @test rss(res) <= 1e-2 - @test r2(res)≈0.95 atol=5e-1 + @test rss(res) <= 1.0e-2 + @test r2(res) ≈ 0.95 atol = 5.0e-1 @test dof(res) == 3 - @test loglikelihood(res) >= 400e3 + @test loglikelihood(res) >= 400.0e3 foreach(get_results(res)) do operator_res @test Matrix(get_operator(operator_res)) ≈ A @test isempty(get_inputmap(operator_res)) @test get_outputmap(operator_res) ≈ I(2) - @test rss(operator_res) <= 1e-10 + @test rss(operator_res) <= 1.0e-10 end end end @@ -92,9 +92,9 @@ end for alg in [DMDPINV(), DMDSVD(), TOTALDMD()] res = solve(prob, alg) @test rss(res) <= 2.0 - @test r2(res)≈1.0 atol=1e-2 + @test r2(res) ≈ 1.0 atol = 1.0e-2 @test dof(res) == 4 - @test loglikelihood(res) >= 85e3 + @test loglikelihood(res) >= 85.0e3 end end end @@ -108,21 +108,21 @@ end sol_ = solve(prob, Tsit5(), saveat = 0.01) # True Rank is 3 - X = Q * sol_[:, :] + 1e-3 * randn(rng, 20, 1001) - DX = Q * sol_(sol_.t, Val{1})[:, :] + 1e-3 * randn(rng, 20, 1001) + X = Q * sol_[:, :] + 1.0e-3 * randn(rng, 20, 1001) + DX = Q * sol_(sol_.t, Val{1})[:, :] + 1.0e-3 * randn(rng, 20, 1001) ddprob = ContinuousDataDrivenProblem(X, sol_.t, DX = DX) for alg in [TOTALDMD(3, DMDPINV()); TOTALDMD(0.01, DMDSVD(3))] res = solve(ddprob, alg, digits = 2) - @test rss(res) <= 1e-1 - @test r2(res)≈1.0 atol=1e-2 + @test rss(res) <= 1.0e-1 + @test r2(res) ≈ 1.0 atol = 1.0e-2 @test dof(res) == 400 - @test loglikelihood(res) >= 99e3 + @test loglikelihood(res) >= 99.0e3 foreach(get_results(res)) do operator_res K = Matrix(get_operator(operator_res)) - @test Q' * K * Q≈K̃ atol=1e-1 - @test Q * K̃ * Q'≈K atol=1e-1 + @test Q' * K * Q ≈ K̃ atol = 1.0e-1 + @test Q * K̃ * Q' ≈ K atol = 1.0e-1 end end end diff --git a/lib/DataDrivenDMD/test/linear_forced.jl b/lib/DataDrivenDMD/test/linear_forced.jl index 85d82cee2..da24655e5 100644 --- a/lib/DataDrivenDMD/test/linear_forced.jl +++ b/lib/DataDrivenDMD/test/linear_forced.jl @@ -44,7 +44,7 @@ end for alg in [DMDPINV(); DMDSVD(); TOTALDMD(2, DMDPINV())] res = solve(ddprob, alg, control_input = B) koopman_result = first(get_results(res)) - @test r2(res)≈0.95 atol=5e-1 + @test r2(res) ≈ 0.95 atol = 5.0e-1 @test dof(res) == 3 @test rss(res) <= eps() @test Matrix(get_operator(koopman_result)) ≈ [1.5 0; 0 0.1] diff --git a/lib/DataDrivenDMD/test/nonlinear_autonomous.jl b/lib/DataDrivenDMD/test/nonlinear_autonomous.jl index 8f6ebee02..ba859907d 100644 --- a/lib/DataDrivenDMD/test/nonlinear_autonomous.jl +++ b/lib/DataDrivenDMD/test/nonlinear_autonomous.jl @@ -57,9 +57,9 @@ end for alg in [DMDPINV(); DMDSVD(); TOTALDMD()] res = solve(prob, Ψ, alg, options = DataDrivenCommonOptions(digits = 2)) - @test get_parameter_values(res.basis)≈[-0.8, 0.63, -0.7] atol=5e-2 + @test get_parameter_values(res.basis) ≈ [-0.8, 0.63, -0.7] atol = 5.0e-2 @test loglikelihood(res) >= 50.0 @test r2(res) >= 0.95 - @test rss(res) <= 1e-1 + @test rss(res) <= 1.0e-1 end end diff --git a/lib/DataDrivenDMD/test/nonlinear_forced.jl b/lib/DataDrivenDMD/test/nonlinear_forced.jl index 09787fcf1..b9a4d2180 100644 --- a/lib/DataDrivenDMD/test/nonlinear_forced.jl +++ b/lib/DataDrivenDMD/test/nonlinear_forced.jl @@ -31,9 +31,9 @@ rng = StableRNG(42) for alg in [DMDPINV(); DMDSVD(); TOTALDMD()] res = solve(prob, Ψ, alg, options = DataDrivenCommonOptions(digits = 1)) koopman_res = first(get_results(res)) - @test r2(res)≈0.95 atol=5e-2 + @test r2(res) ≈ 0.95 atol = 5.0e-2 @test dof(res) == 4 - @test get_inputmap(koopman_res)≈[0; 1.0; 0.0;;] atol=1e-1 + @test get_inputmap(koopman_res) ≈ [0; 1.0; 0.0;;] atol = 1.0e-1 @test get_outputmap(koopman_res) ≈ [1.0 0 0; 0 1 0] @test get_parameter_values(res.basis) ≈ [-2.0; -0.4; 0.5; 0.9] end @@ -63,9 +63,9 @@ end for alg in [DMDPINV(); DMDSVD(); TOTALDMD()] res = solve(prob, Ψ, alg, options = DataDrivenCommonOptions(digits = 1)) koopman_res = first(get_results(res)) - @test r2(res)≈0.95 atol=5e-2 + @test r2(res) ≈ 0.95 atol = 5.0e-2 @test dof(res) == 4 - @test get_inputmap(koopman_res)≈[0; 0.0; 1.0;;] atol=1e-1 + @test get_inputmap(koopman_res) ≈ [0; 0.0; 1.0;;] atol = 1.0e-1 @test get_outputmap(koopman_res) ≈ [1.0 0 0; 0 1 1] end end diff --git a/lib/DataDrivenLux/src/DataDrivenLux.jl b/lib/DataDrivenLux/src/DataDrivenLux.jl index d82bbe5a6..f3d7fe70b 100644 --- a/lib/DataDrivenLux/src/DataDrivenLux.jl +++ b/lib/DataDrivenLux/src/DataDrivenLux.jl @@ -4,9 +4,9 @@ using DataDrivenDiffEq # Load specific (abstract) types using DataDrivenDiffEq: AbstractBasis, AbstractDataDrivenAlgorithm, - AbstractDataDrivenResult, AbstractDataDrivenProblem, DDReturnCode, - ABSTRACT_CONT_PROB, ABSTRACT_DISCRETE_PROB, - InternalDataDrivenProblem, is_implicit, is_controlled + AbstractDataDrivenResult, AbstractDataDrivenProblem, DDReturnCode, + ABSTRACT_CONT_PROB, ABSTRACT_DISCRETE_PROB, + InternalDataDrivenProblem, is_implicit, is_controlled using DocStringExtensions: DocStringExtensions, FIELDS, TYPEDEF, SIGNATURES using CommonSolve: CommonSolve, solve! @@ -23,7 +23,7 @@ using WeightInitializers: WeightInitializers, ones32, zeros32 using InverseFunctions: InverseFunctions, NoInverse using TransformVariables: TransformVariables, as, transform_logdensity using Distributions: Distributions, Distribution, Normal, Uniform, Univariate, dof, - loglikelihood, logpdf, mean, mode, quantile, scale, truncated + loglikelihood, logpdf, mean, mode, quantile, scale, truncated using DistributionsAD: DistributionsAD using StatsBase: StatsBase, aicc, nobs, nullloglikelihood, r2, rss, sum, weights @@ -58,7 +58,7 @@ end ## include("utils.jl") -## +## include("custom_priors.jl") export AdditiveError, MultiplicativeError export ObservedModel diff --git a/lib/DataDrivenLux/src/algorithms/crossentropy.jl b/lib/DataDrivenLux/src/algorithms/crossentropy.jl index 300303d7c..d240586f1 100644 --- a/lib/DataDrivenLux/src/algorithms/crossentropy.jl +++ b/lib/DataDrivenLux/src/algorithms/crossentropy.jl @@ -8,15 +8,20 @@ $(SIGNATURES) Uses the crossentropy method for discrete optimization to search the space of possible solutions. """ -function CrossEntropy(; populationsize = 100, functions = (sin, exp, cos, log, +, -, /, *), +function CrossEntropy(; + populationsize = 100, functions = (sin, exp, cos, log, +, -, /, *), arities = (1, 1, 1, 1, 2, 2, 2, 2), n_layers = 1, skip = true, loss = aicc, keep = 0.1, use_protected = true, distributed = false, threaded = false, rng = Random.default_rng(), optimizer = LBFGS(), optim_options = Optim.Options(), - observed = nothing, alpha = 0.999f0) - return CrossEntropy(CommonAlgOptions(; - populationsize, functions, arities, n_layers, skip, simplex = DirectSimplex(), loss, - keep, use_protected, distributed, threaded, rng, optimizer, - optim_options, optimiser = nothing, observed, alpha)) + observed = nothing, alpha = 0.999f0 + ) + return CrossEntropy( + CommonAlgOptions(; + populationsize, functions, arities, n_layers, skip, simplex = DirectSimplex(), loss, + keep, use_protected, distributed, threaded, rng, optimizer, + optim_options, optimiser = nothing, observed, alpha + ) + ) end Base.print(io::IO, ::CrossEntropy) = print(io, "CrossEntropy()") @@ -28,7 +33,7 @@ function init_model(x::CrossEntropy, basis::Basis, dataset::Dataset, intervals) # Get the parameter mapping variable_mask = map(enumerate(equations(basis))) do (i, eq) return any(ModelingToolkit.isvariable, ModelingToolkit.get_variables(eq.rhs)) && - IntervalArithmetic.iscommon(intervals[i]) + IntervalArithmetic.iscommon(intervals[i]) end variable_mask = Any[variable_mask...] @@ -37,14 +42,18 @@ function init_model(x::CrossEntropy, basis::Basis, dataset::Dataset, intervals) functions = map(convert_to_safe, functions) end - return LayeredDAG(length(basis), size(dataset.y, 1), n_layers, arities, functions; - skip, input_functions = variable_mask, x.options.simplex) + return LayeredDAG( + length(basis), size(dataset.y, 1), n_layers, arities, functions; + skip, input_functions = variable_mask, x.options.simplex + ) end function update_parameters!(cache::SearchCache{<:CrossEntropy}) - p̄ = mean(map(cache.candidates[cache.keeps]) do candidate - return ComponentVector(get_configuration(candidate.model.model, cache.p, candidate.st)) - end) + p̄ = mean( + map(cache.candidates[cache.keeps]) do candidate + return ComponentVector(get_configuration(candidate.model.model, cache.p, candidate.st)) + end + ) alpha = cache.alg.options.alpha @. cache.p = alpha * cache.p + (true - alpha) * p̄ return diff --git a/lib/DataDrivenLux/src/algorithms/randomsearch.jl b/lib/DataDrivenLux/src/algorithms/randomsearch.jl index 7f789246b..8d7e56de1 100644 --- a/lib/DataDrivenLux/src/algorithms/randomsearch.jl +++ b/lib/DataDrivenLux/src/algorithms/randomsearch.jl @@ -8,15 +8,20 @@ $(SIGNATURES) Performs a random search over the space of possible solutions to the symbolic regression problem. """ -function RandomSearch(; populationsize = 100, functions = (sin, exp, cos, log, +, -, /, *), +function RandomSearch(; + populationsize = 100, functions = (sin, exp, cos, log, +, -, /, *), arities = (1, 1, 1, 1, 2, 2, 2, 2), n_layers = 1, skip = true, loss = aicc, keep = 0.1, use_protected = true, distributed = false, threaded = false, rng = Random.default_rng(), optimizer = LBFGS(), optim_options = Optim.Options(), - observed = nothing, alpha = 0.999f0) - return RandomSearch(CommonAlgOptions(; - populationsize, functions, arities, n_layers, skip, simplex = Softmax(), loss, - keep, use_protected, distributed, threaded, rng, optimizer, - optim_options, optimiser = nothing, observed, alpha)) + observed = nothing, alpha = 0.999f0 + ) + return RandomSearch( + CommonAlgOptions(; + populationsize, functions, arities, n_layers, skip, simplex = Softmax(), loss, + keep, use_protected, distributed, threaded, rng, optimizer, + optim_options, optimiser = nothing, observed, alpha + ) + ) end Base.print(io::IO, ::RandomSearch) = print(io, "RandomSearch") diff --git a/lib/DataDrivenLux/src/algorithms/reinforce.jl b/lib/DataDrivenLux/src/algorithms/reinforce.jl index 66f53ccf5..7570417e1 100644 --- a/lib/DataDrivenLux/src/algorithms/reinforce.jl +++ b/lib/DataDrivenLux/src/algorithms/reinforce.jl @@ -10,18 +10,23 @@ $(SIGNATURES) Uses the REINFORCE algorithm to search over the space of possible solutions to the symbolic regression problem. """ -function Reinforce(; reward = RelativeReward(false), populationsize = 100, +function Reinforce(; + reward = RelativeReward(false), populationsize = 100, functions = (sin, exp, cos, log, +, -, /, *), arities = (1, 1, 1, 1, 2, 2, 2, 2), n_layers = 1, skip = true, loss = aicc, keep = 0.1, use_protected = true, distributed = false, threaded = false, rng = Random.default_rng(), optimizer = LBFGS(), optim_options = Optim.Options(), observed = nothing, - alpha = 0.999f0, optimiser = Adam(), ad_backend = AD.ForwardDiffBackend()) - return Reinforce(reward, + alpha = 0.999f0, optimiser = Adam(), ad_backend = AD.ForwardDiffBackend() + ) + return Reinforce( + reward, ad_backend, CommonAlgOptions(; populationsize, functions, arities, n_layers, skip, simplex = Softmax(), loss, keep, use_protected, distributed, threaded, rng, optimizer, - optim_options, optimiser, observed, alpha)) + optim_options, optimiser, observed, alpha + ) + ) end Base.print(io::IO, ::Reinforce) = print(io, "Reinforce") @@ -31,9 +36,11 @@ function reinforce_loss(candidates, p, alg) losses = map(alg.options.loss, candidates) rewards = alg.reward(losses) # ∇U(θ) = E[∇log(p)*R(t)] - return mean(map(enumerate(candidates)) do (i, candidate) - return rewards[i] * -candidate(p) - end) + return mean( + map(enumerate(candidates)) do (i, candidate) + return rewards[i] * -candidate(p) + end + ) end function update_parameters!(cache::SearchCache{<:Reinforce}) diff --git a/lib/DataDrivenLux/src/caches/cache.jl b/lib/DataDrivenLux/src/caches/cache.jl index b751b32b1..41522af0d 100644 --- a/lib/DataDrivenLux/src/caches/cache.jl +++ b/lib/DataDrivenLux/src/caches/cache.jl @@ -17,7 +17,7 @@ function init_model(x::AbstractDAGSRAlgorithm, basis::Basis, dataset::Dataset, i # Get the parameter mapping variable_mask = map(enumerate(equations(basis))) do (i, eq) return any(ModelingToolkit.isvariable, ModelingToolkit.get_variables(eq.rhs)) && - IntervalArithmetic.iscommon(intervals[i]) + IntervalArithmetic.iscommon(intervals[i]) end variable_mask = Any[variable_mask...] @@ -26,14 +26,20 @@ function init_model(x::AbstractDAGSRAlgorithm, basis::Basis, dataset::Dataset, i functions = map(convert_to_safe, functions) end - return LayeredDAG(length(basis), size(dataset.y, 1), n_layers, arities, functions; - skip = skip, input_functions = variable_mask, simplex = simplex) + return LayeredDAG( + length(basis), size(dataset.y, 1), n_layers, arities, functions; + skip = skip, input_functions = variable_mask, simplex = simplex + ) end -function init_cache(x::X where {X <: AbstractDAGSRAlgorithm}, - basis::Basis, problem::DataDrivenProblem; kwargs...) - (; rng, keep, observed, populationsize, optimizer, - optim_options, optimiser, loss) = x.options +function init_cache( + x::X where {X <: AbstractDAGSRAlgorithm}, + basis::Basis, problem::DataDrivenProblem; kwargs... + ) + (; + rng, keep, observed, populationsize, optimizer, + optim_options, optimiser, loss, + ) = x.options # Derive the model dataset = Dataset(problem) TData = eltype(dataset) @@ -41,7 +47,7 @@ function init_cache(x::X where {X <: AbstractDAGSRAlgorithm}, rng_ = Lux.replicate(rng) observed = isa(observed, ObservedModel) ? observed : - ObservedModel(dataset.y, fixed = true) + ObservedModel(dataset.y, fixed = true) parameters = ParameterDistributions(basis, TData) @@ -51,12 +57,15 @@ function init_cache(x::X where {X <: AbstractDAGSRAlgorithm}, ps = ComponentVector(Lux.initialparameters(rng_, model)) - # Derive the candidates + # Derive the candidates candidates = map(1:populationsize) do i - candidate = Candidate(rng_, model, basis, dataset; observed = observed, - parameterdist = parameters, ptype = TData) + candidate = Candidate( + rng_, model, basis, dataset; observed = observed, + parameterdist = parameters, ptype = TData + ) optimize_candidate!( - candidate, dataset; optimizer = optimizer, options = optim_options) + candidate, dataset; optimizer = optimizer, options = optim_options + ) return candidate end @@ -91,7 +100,8 @@ function init_cache(x::X where {X <: AbstractDAGSRAlgorithm}, optimiser_state = nothing end return SearchCache{typeof(x), ptype, typeof(optimiser_state)}( - x, candidates, ages, keeps, sorting, ps, dataset, optimiser_state) + x, candidates, ages, keeps, sorting, ps, dataset, optimiser_state + ) end function update_cache!(cache::SearchCache) @@ -123,7 +133,7 @@ end # Optimizes the cache and returns the loglikelihoods -# Serial +# Serial function optimize_cache!(cache::SearchCache{<:Any, __PROCESSUSE(1)}, p = cache.p) (; optimizer, optim_options) = cache.alg.options map(enumerate(cache.candidates)) do (i, candidate) @@ -132,7 +142,8 @@ function optimize_cache!(cache::SearchCache{<:Any, __PROCESSUSE(1)}, p = cache.p return true else optimize_candidate!( - candidate, cache.dataset, p; optimizer = optimizer, options = optim_options) + candidate, cache.dataset, p; optimizer = optimizer, options = optim_options + ) cache.ages[i] = 0 return true end @@ -143,13 +154,15 @@ end # Threaded function optimize_cache!(cache::SearchCache{<:Any, __PROCESSUSE(2)}, p = cache.p) (; optimizer, optim_options) = cache.alg.options - # Update all + # Update all Threads.@threads for i in 1:length(cache.keeps) if cache.keeps[i] cache.ages[i] += 1 else - optimize_candidate!(cache.candidates[i], cache.dataset, p; - optimizer = optimizer, options = optim_options) + optimize_candidate!( + cache.candidates[i], cache.dataset, p; + optimizer = optimizer, options = optim_options + ) cache.ages[i] = 0 end end @@ -165,8 +178,10 @@ function optimize_cache!(cache::SearchCache{<:Any, __PROCESSUSE(3)}, p = cache.p cache.ages[i] += 1 return true else - optimize_candidate!(cache.candidates[i], cache.dataset, p; - optimizer = optimizer, options = optim_options) + optimize_candidate!( + cache.candidates[i], cache.dataset, p; + optimizer = optimizer, options = optim_options + ) cache.ages[i] = 0 return true end diff --git a/lib/DataDrivenLux/src/caches/candidate.jl b/lib/DataDrivenLux/src/caches/candidate.jl index fbd553391..45edcf48a 100644 --- a/lib/DataDrivenLux/src/caches/candidate.jl +++ b/lib/DataDrivenLux/src/caches/candidate.jl @@ -7,7 +7,8 @@ mutable struct PathStatistics{T} <: StatsBase.StatisticalModel end function update_stats!( - stats::PathStatistics{T}, rss::T, ll::T, nullll::T, dof::Int) where {T} + stats::PathStatistics{T}, rss::T, ll::T, nullll::T, dof::Int + ) where {T} stats.dof = dof stats.loglikelihood = ll stats.nullloglikelihood = nullll @@ -27,8 +28,10 @@ StatsBase.r2(c::PathStatistics) = r2(c, :CoxSnell) model end -function (c::ComponentModel)(dataset::Dataset{T}, ps, st::NamedTuple, - p::AbstractVector{T}) where {T} +function (c::ComponentModel)( + dataset::Dataset{T}, ps, st::NamedTuple, + p::AbstractVector{T} + ) where {T} return first(c.model(c.basis(dataset, p), ps, st)) end function (c::ComponentModel)(ps, st::NamedTuple, paths::Vector{<:AbstractPathState}) @@ -91,12 +94,15 @@ get_scales(c::Candidate) = transform_scales(c.observed, c.scales) function Candidate( rng, model, basis, dataset::Dataset{T}; observed = ObservedModel(dataset.y), - parameterdist = ParameterDistributions(basis), ptype = Float32) where {T} + parameterdist = ParameterDistributions(basis), ptype = Float32 + ) where {T} # Create the initial state and path dataset_intervals = interval_eval(basis, dataset, get_interval(parameterdist)) - incoming_path = [PathState{ptype}(dataset_intervals[i], (), ((0, i),)) - for i in 1:length(basis)] + incoming_path = [ + PathState{ptype}(dataset_intervals[i], (), ((0, i),)) + for i in 1:length(basis) + ] ps, st = Lux.setup(rng, model) outgoing_path, st = sample(model, incoming_path, ps, st) @@ -119,9 +125,11 @@ function Candidate( stats = PathStatistics(rss, lls, null_ll, dof_, prod(size(dataset.y))) - return Candidate(Lux.replicate(rng), st, ComponentVector(ps), incoming_path, + return Candidate( + Lux.replicate(rng), st, ComponentVector(ps), incoming_path, outgoing_path, stats, observed, parameterdist, scales, parameters, - ComponentModel(basis, model)) + ComponentModel(basis, model) + ) end function update_values!(c::Candidate, ps, dataset) @@ -140,13 +148,16 @@ function update_values!(c::Candidate, ps, dataset) end @views function Distributions.logpdf( - c::Candidate, p::ComponentVector, dataset::Dataset{T}, ps = c.ps) where {T} + c::Candidate, p::ComponentVector, dataset::Dataset{T}, ps = c.ps + ) where {T} ŷ = c(dataset, ps, p.parameters) return logpdf(c, p, dataset.y, ŷ) end -function Distributions.logpdf(c::Candidate, p::AbstractVector, y::AbstractMatrix{T}, - ŷ::AbstractMatrix{T}) where {T} +function Distributions.logpdf( + c::Candidate, p::AbstractVector, y::AbstractMatrix{T}, + ŷ::AbstractMatrix{T} + ) where {T} return logpdf(c.observed, y, ŷ, p.scales) + logpdf(c.parameterdist, p.parameters) end @@ -154,7 +165,8 @@ initial_values(c::Candidate) = ComponentVector(; c.scales, c.parameters) function optimize_candidate!( c::Candidate, dataset::Dataset{T}, ps = c.ps; optimizer = Optim.LBFGS(), - options::Optim.Options = Optim.Options()) where {T} + options::Optim.Options = Optim.Options() + ) where {T} path, st = sample(c, ps) p_init = initial_values(c) @@ -199,7 +211,8 @@ end get_nodes(c::Candidate) = @ignore_derivatives get_nodes(c.outgoing_path) function convert_to_basis( - candidate::Candidate, ps = candidate.ps, options = DataDrivenCommonOptions()) + candidate::Candidate, ps = candidate.ps, options = DataDrivenCommonOptions() + ) (; basis, model) = candidate.model (; eval_expresssion) = options p_best = get_parameters(candidate) @@ -215,8 +228,10 @@ function convert_to_basis( eqs = collect(map(eq -> ModelingToolkit.substitute(eq, subs), eqs)) - return Basis(eqs, states(basis), parameters = p_new, iv = get_iv(basis), + return Basis( + eqs, states(basis), parameters = p_new, iv = get_iv(basis), controls = controls(basis), observed = observed(basis), implicits = implicit_variables(basis), - name = gensym(:Basis), eval_expression = eval_expresssion) + name = gensym(:Basis), eval_expression = eval_expresssion + ) end diff --git a/lib/DataDrivenLux/src/caches/dataset.jl b/lib/DataDrivenLux/src/caches/dataset.jl index 29bf894ca..e57c42c8c 100644 --- a/lib/DataDrivenLux/src/caches/dataset.jl +++ b/lib/DataDrivenLux/src/caches/dataset.jl @@ -11,9 +11,11 @@ end Base.eltype(::Dataset{T}) where {T} = T -function Dataset(X::AbstractMatrix, Y::AbstractMatrix, +function Dataset( + X::AbstractMatrix, Y::AbstractMatrix, U::AbstractMatrix = Array{eltype(X)}(undef, 0, 0), - t::AbstractVector = Array{eltype(X)}(undef, 0)) + t::AbstractVector = Array{eltype(X)}(undef, 0) + ) T = Base.promote_eltype(X, Y, U, t) X = convert.(T, X) Y = convert.(T, Y) diff --git a/lib/DataDrivenLux/src/custom_priors.jl b/lib/DataDrivenLux/src/custom_priors.jl index 9e7acee89..a51073a6d 100644 --- a/lib/DataDrivenLux/src/custom_priors.jl +++ b/lib/DataDrivenLux/src/custom_priors.jl @@ -7,8 +7,10 @@ An error following `ŷ ~ y + ϵ`. struct AdditiveError <: AbstractErrorModel end -function (x::AdditiveError)(d::D, y::T, ỹ::R, - scale::S = one(T)) where {D <: Type, T <: Number, S <: Number, R <: Number} +function (x::AdditiveError)( + d::D, y::T, ỹ::R, + scale::S = one(T) + ) where {D <: Type, T <: Number, S <: Number, R <: Number} return logpdf(d(y, scale), ỹ) end @@ -20,8 +22,10 @@ An error following `ŷ ~ y * (1+ϵ)`. struct MultiplicativeError <: AbstractErrorModel end -function (x::MultiplicativeError)(d::D, y::T, ỹ::R, - scale::S = one(T)) where {D <: Type, T <: Number, S <: Number, R <: Number} +function (x::MultiplicativeError)( + d::D, y::T, ỹ::R, + scale::S = one(T) + ) where {D <: Type, T <: Number, S <: Number, R <: Number} return logpdf(d(y, abs(y) * scale), ỹ) end @@ -34,9 +38,11 @@ end scale_transformation end -function ObservedDistribution(::Type{D}, errormodel::AbstractErrorModel; fixed = false, - transform = as(Real, 1e-5, TransformVariables.∞), - scale = 1.0) where {D <: Distributions.Distribution{Univariate, <:Any}} +function ObservedDistribution( + ::Type{D}, errormodel::AbstractErrorModel; fixed = false, + transform = as(Real, 1.0e-5, TransformVariables.∞), + scale = 1.0 + ) where {D <: Distributions.Distribution{Univariate, <:Any}} latent_scale = TransformVariables.inverse(transform, scale) return ObservedDistribution{fixed, D}(errormodel, latent_scale, transform) end @@ -54,35 +60,58 @@ get_dist(::ObservedDistribution{<:Any, D}) where {D} = D Base.show(io::IO, d::ObservedDistribution) = summary(io, d) function Distributions.logpdf( - d::ObservedDistribution{false}, x::X, x̂::Y, scale::Number) where {X, Y} - return sum(map( - xs -> d.errormodel( - get_dist(d), xs..., TransformVariables.transform(d.scale_transformation, scale)), - zip(x, x̂))) + d::ObservedDistribution{false}, x::X, x̂::Y, scale::Number + ) where {X, Y} + return sum( + map( + xs -> d.errormodel( + get_dist(d), xs..., TransformVariables.transform(d.scale_transformation, scale) + ), + zip(x, x̂) + ) + ) end function Distributions.logpdf( - d::ObservedDistribution{true}, x::X, x̂::Y, ::Number) where {X, Y} - return sum(map( - xs -> d.errormodel(get_dist(d), xs..., - TransformVariables.transform(d.scale_transformation, d.latent_scale)), - zip(x, x̂))) + d::ObservedDistribution{true}, x::X, x̂::Y, ::Number + ) where {X, Y} + return sum( + map( + xs -> d.errormodel( + get_dist(d), xs..., + TransformVariables.transform(d.scale_transformation, d.latent_scale) + ), + zip(x, x̂) + ) + ) end function Distributions.logpdf( - d::ObservedDistribution{false}, x::X, x̂::Number, scale::Number) where {X} - return sum(map( - xs -> d.errormodel(get_dist(d), xs, x̂, - TransformVariables.transform(d.scale_transformation, scale)), - x)) + d::ObservedDistribution{false}, x::X, x̂::Number, scale::Number + ) where {X} + return sum( + map( + xs -> d.errormodel( + get_dist(d), xs, x̂, + TransformVariables.transform(d.scale_transformation, scale) + ), + x + ) + ) end function Distributions.logpdf( - d::ObservedDistribution{true}, x::X, x̂::Number, ::Number) where {X} - return sum(map( - xs -> d.errormodel(get_dist(d), xs, x̂, - TransformVariables.transform(d.scale_transformation, d.latent_scale)), - x)) + d::ObservedDistribution{true}, x::X, x̂::Number, ::Number + ) where {X} + return sum( + map( + xs -> d.errormodel( + get_dist(d), xs, x̂, + TransformVariables.transform(d.scale_transformation, d.latent_scale) + ), + x + ) + ) end function transform_scales(d::ObservedDistribution, scale::Number) @@ -114,16 +143,22 @@ end Base.show(io::IO, o::ObservedModel) = summary(io, o) -function Distributions.logpdf(o::ObservedModel{M}, x::AbstractMatrix, x̂::AbstractMatrix, - scales::AbstractVector = ones(eltype(x̂), size(x, 1))) where {M} +function Distributions.logpdf( + o::ObservedModel{M}, x::AbstractMatrix, x̂::AbstractMatrix, + scales::AbstractVector = ones(eltype(x̂), size(x, 1)) + ) where {M} return sum(map(logpdf, o.observed_distributions, eachrow(x), eachrow(x̂), scales)) end -function Distributions.logpdf(o::ObservedModel{M}, x::AbstractMatrix, x̂::AbstractVector, - scales::AbstractVector = ones(eltype(x̂), size(x, 1))) where {M} - sum(map(axes(x, 1)) do i - return logpdf(o.observed_distributions[i], x[i, :], x̂[i], scales[i]) - end) +function Distributions.logpdf( + o::ObservedModel{M}, x::AbstractMatrix, x̂::AbstractVector, + scales::AbstractVector = ones(eltype(x̂), size(x, 1)) + ) where {M} + return sum( + map(axes(x, 1)) do i + return logpdf(o.observed_distributions[i], x[i, :], x̂[i], scales[i]) + end + ) end get_init(o::ObservedModel) = collect(map(get_init, o.observed_distributions)) @@ -141,7 +176,8 @@ end end function ParameterDistribution( - d::Distribution{Univariate}, init = mean(d), type::Type{T} = Float64) where {T} + d::Distribution{Univariate}, init = mean(d), type::Type{T} = Float64 + ) where {T} lower, upper = convert.(T, extrema(d)) lower_t = isinf(lower) ? -TransformVariables.∞ : lower upper_t = isinf(upper) ? TransformVariables.∞ : upper @@ -165,7 +201,7 @@ function Distributions.logpdf(p::ParameterDistribution, pval::T) where {T <: Num return transform_logdensity(p.transformation, Base.Fix1(logpdf, p.distribution), pval) end -# Parameters +# Parameters struct ParameterDistributions{T, N} distributions::NTuple{N, ParameterDistribution} diff --git a/lib/DataDrivenLux/src/lux/graph.jl b/lib/DataDrivenLux/src/lux/graph.jl index db721837f..355a1e259 100644 --- a/lib/DataDrivenLux/src/lux/graph.jl +++ b/lib/DataDrivenLux/src/lux/graph.jl @@ -12,15 +12,21 @@ $(FIELDS) layers end -function LayeredDAG(in_dimension::Int, out_dimension::Int, n_layers::Int, - fs::Vector{Pair{Function, Int}}; kwargs...) - return LayeredDAG(in_dimension, out_dimension, n_layers, - tuple(last.(fs)...), tuple(first.(fs)...); kwargs...) +function LayeredDAG( + in_dimension::Int, out_dimension::Int, n_layers::Int, + fs::Vector{Pair{Function, Int}}; kwargs... + ) + return LayeredDAG( + in_dimension, out_dimension, n_layers, + tuple(last.(fs)...), tuple(first.(fs)...); kwargs... + ) end -function LayeredDAG(in_dimension::Int, out_dimension::Int, n_layers::Int, +function LayeredDAG( + in_dimension::Int, out_dimension::Int, n_layers::Int, arities::Tuple, fs::Tuple; skip = false, eltype::Type{T} = Float32, - input_functions = Any[identity for i in 1:in_dimension], kwargs...) where {T} + input_functions = Any[identity for i in 1:in_dimension], kwargs... + ) where {T} n_inputs = in_dimension input_functions = copy(input_functions) @@ -33,8 +39,10 @@ function LayeredDAG(in_dimension::Int, out_dimension::Int, n_layers::Int, valid_idxs .= (arities .<= n_inputs) - layer = FunctionLayer(n_inputs, arities[valid_idxs], fs[valid_idxs]; skip = skip, - id_offset = i, input_functions = input_functions, kwargs...) + layer = FunctionLayer( + n_inputs, arities[valid_idxs], fs[valid_idxs]; skip = skip, + id_offset = i, input_functions = input_functions, kwargs... + ) if skip n_inputs = n_inputs + sum(valid_idxs) @@ -48,10 +56,14 @@ function LayeredDAG(in_dimension::Int, out_dimension::Int, n_layers::Int, return push!(layers, layer) end # The last layer is a decision node which uses an identity - push!(layers, - FunctionLayer(n_inputs, Tuple(1 for i in 1:out_dimension), + push!( + layers, + FunctionLayer( + n_inputs, Tuple(1 for i in 1:out_dimension), Tuple(identity for i in 1:out_dimension); skip = false, - input_functions = input_functions, id_offset = n_layers + 1, kwargs...)) + input_functions = input_functions, id_offset = n_layers + 1, kwargs... + ) + ) return LayeredDAG(Lux.Chain(layers...)) end @@ -66,11 +78,15 @@ end function get_loglikelihood(c::LayeredDAG, ps, st, paths::Vector{<:AbstractPathState}) lls = get_loglikelihood(c, ps, st) - sum(map(paths) do path - nodes = get_nodes(path) - sum(map(nodes) do (i, j) - i > 0 && return lls[i][j] - return 0.0f0 - end) - end) + return sum( + map(paths) do path + nodes = get_nodes(path) + sum( + map(nodes) do (i, j) + i > 0 && return lls[i][j] + return 0.0f0 + end + ) + end + ) end diff --git a/lib/DataDrivenLux/src/lux/layer.jl b/lib/DataDrivenLux/src/lux/layer.jl index 4b49dc7ec..d2985cbbe 100644 --- a/lib/DataDrivenLux/src/lux/layer.jl +++ b/lib/DataDrivenLux/src/lux/layer.jl @@ -15,21 +15,26 @@ end function FunctionLayer( in_dimension::Int, arities::Tuple, fs::Tuple; skip = false, id_offset = 1, - input_functions = Any[identity for i in 1:in_dimension], kwargs...) + input_functions = Any[identity for i in 1:in_dimension], kwargs... + ) nodes = map(eachindex(arities)) do i # We check if we have an inverse here - return FunctionNode(fs[i], arities[i], in_dimension, (id_offset, i); - input_functions, kwargs...) + return FunctionNode( + fs[i], arities[i], in_dimension, (id_offset, i); + input_functions, kwargs... + ) end inner_model = Lux.Chain(Lux.BranchLayer(nodes...), Lux.WrappedFunction(splat(vcat))) return FunctionLayer( - skip ? Lux.Parallel(vcat, inner_model, Lux.NoOpLayer()) : inner_model, skip) + skip ? Lux.Parallel(vcat, inner_model, Lux.NoOpLayer()) : inner_model, skip + ) end function get_loglikelihood(r::FunctionLayer, ps, st) if r.skip return _get_layer_loglikelihood( - r.nodes.layers[1].layers[1].layers, ps.layer_1.layer_1, st.layer_1.layer_1) + r.nodes.layers[1].layers[1].layers, ps.layer_1.layer_1, st.layer_1.layer_1 + ) else return _get_layer_loglikelihood(r.nodes.layers[1].layers, ps.layer_1, st.layer_1) end @@ -38,28 +43,41 @@ end function get_configuration(r::FunctionLayer, ps, st) if r.skip return _get_configuration( - r.nodes.layers[1].layers[1].layers, ps.layer_1.layer_1, st.layer_1.layer_1) + r.nodes.layers[1].layers[1].layers, ps.layer_1.layer_1, st.layer_1.layer_1 + ) else return _get_configuration(r.nodes.layers[1].layers, ps.layer_1, st.layer_1) end end @generated function _get_layer_loglikelihood( - layers::NamedTuple{fields}, ps, st::NamedTuple{fields}) where {fields} + layers::NamedTuple{fields}, ps, st::NamedTuple{fields} + ) where {fields} N = length(fields) st_symbols = [gensym() for _ in 1:N] - calls = [:($(st_symbols[i]) = get_loglikelihood( - layers.$(fields[i]), ps.$(fields[i]), st.$(fields[i]))) for i in 1:N] + calls = [ + :( + $(st_symbols[i]) = get_loglikelihood( + layers.$(fields[i]), ps.$(fields[i]), st.$(fields[i]) + ) + ) for i in 1:N + ] push!(calls, :(st = NamedTuple{$fields}((($(Tuple(st_symbols)...),))))) return Expr(:block, calls...) end @generated function _get_configuration( - layers::NamedTuple{fields}, ps, st::NamedTuple{fields}) where {fields} + layers::NamedTuple{fields}, ps, st::NamedTuple{fields} + ) where {fields} N = length(fields) st_symbols = [gensym() for _ in 1:N] - calls = [:($(st_symbols[i]) = get_configuration( - layers.$(fields[i]), ps.$(fields[i]), st.$(fields[i]))) for i in 1:N] + calls = [ + :( + $(st_symbols[i]) = get_configuration( + layers.$(fields[i]), ps.$(fields[i]), st.$(fields[i]) + ) + ) for i in 1:N + ] push!(calls, :(st = NamedTuple{$fields}((($(Tuple(st_symbols)...),))))) return Expr(:block, calls...) end diff --git a/lib/DataDrivenLux/src/lux/node.jl b/lib/DataDrivenLux/src/lux/node.jl index 12e505ea8..cec574e96 100644 --- a/lib/DataDrivenLux/src/lux/node.jl +++ b/lib/DataDrivenLux/src/lux/node.jl @@ -1,4 +1,3 @@ - """ $(TYPEDEF) @@ -39,13 +38,15 @@ function mask_inverse(::typeof(identity), arity::Int, in_f::AbstractVector) return ones(Bool, length(in_f)) end -function FunctionNode(f::F, arity::Int, input_dimension::Int, +function FunctionNode( + f::F, arity::Int, input_dimension::Int, id::Union{Int, NTuple{<:Any, Int}}; skip = false, simplex = Softmax(), - input_functions = [identity for i in 1:input_dimension], kwargs...) where {F} + input_functions = [identity for i in 1:input_dimension], kwargs... + ) where {F} input_mask = mask_inverse(f, arity, input_functions) - @assert sum(input_mask)>=1 "Input masks should enable at least one choice." - @assert length(input_mask)==input_dimension "Input dimension should be sized equally \ + @assert sum(input_mask) >= 1 "Input masks should enable at least one choice." + @assert length(input_mask) == input_dimension "Input dimension should be sized equally \ to input_mask" internal_node = InternalFunctionNode{id}(f, arity, input_dimension, simplex, input_mask) @@ -62,8 +63,10 @@ end function LuxCore.initialstates(rng::AbstractRNG, p::InternalFunctionNode) rand(rng) rng_ = LuxCore.replicate(rng) - return (; priors = init_weights(p.simplex, rng, sum(p.input_mask), p.arity), - active_inputs = zeros(Int, p.arity), temperature = 1.0f0, rng = rng_) + return (; + priors = init_weights(p.simplex, rng, sum(p.input_mask), p.arity), + active_inputs = zeros(Int, p.arity), temperature = 1.0f0, rng = rng_, + ) end @views function update_state(p::InternalFunctionNode, ps, st) @@ -99,9 +102,11 @@ end get_temperature(::FunctionNode, ps, st) = st.temperature function get_loglikelihood(::FunctionNode, ps, st) - return sum(map(enumerate(eachcol(ps.weights))) do (i, weight) - return logsoftmax(weight ./ st.temperature)[st.active_inputs[i]] - end) + return sum( + map(enumerate(eachcol(ps.weights))) do (i, weight) + return logsoftmax(weight ./ st.temperature)[st.active_inputs[i]] + end + ) end get_inputs(::FunctionNode, ps, st) = st.active_inputs diff --git a/lib/DataDrivenLux/src/lux/path_state.jl b/lib/DataDrivenLux/src/lux/path_state.jl index 377d735a8..a0669902a 100644 --- a/lib/DataDrivenLux/src/lux/path_state.jl +++ b/lib/DataDrivenLux/src/lux/path_state.jl @@ -9,11 +9,13 @@ struct PathState{T, PO <: Tuple, PI <: Tuple} <: AbstractPathState path_ids::PI function PathState{T}( - interval::Interval{T}, path_operators::PO, path_ids::PI) where {T, PO, PI} + interval::Interval{T}, path_operators::PO, path_ids::PI + ) where {T, PO, PI} return new{T, PO, PI}(interval, path_operators, path_ids) end function PathState{T}( - interval::Interval, path_operators::PO, path_ids::PI) where {T, PO, PI} + interval::Interval, path_operators::PO, path_ids::PI + ) where {T, PO, PI} return new{T, PO, PI}(Interval{T}(interval), path_operators, path_ids) end end @@ -32,21 +34,26 @@ get_nodes(state::PathState) = state.path_ids @inline tuplejoin(x, y, z...) = tuplejoin(tuplejoin(x, y), z...) function update_path( - f::F where {F <: Function}, id::Tuple{Int, Int}, state::PathState{T}) where {T} + f::F where {F <: Function}, id::Tuple{Int, Int}, state::PathState{T} + ) where {T} return PathState{T}( - f(get_interval(state)), (f, get_operators(state)...), (id, get_nodes(state)...)) + f(get_interval(state)), (f, get_operators(state)...), (id, get_nodes(state)...) + ) end function update_path(::Nothing, id::Tuple{Int, Int}, state::PathState{T}) where {T} return PathState{T}( - get_interval(state), (identity, get_operators(state)...), (id, get_nodes(state)...)) + get_interval(state), (identity, get_operators(state)...), (id, get_nodes(state)...) + ) end function update_path( - f::F where {F <: Function}, id::Tuple{Int, Int}, states::PathState{T}...) where {T} + f::F where {F <: Function}, id::Tuple{Int, Int}, states::PathState{T}... + ) where {T} return PathState{T}( f(get_interval.(states)...), (f, tuplejoin(map(get_operators, states)...)...), - (id, tuplejoin(map(get_nodes, states)...)...)) + (id, tuplejoin(map(get_nodes, states)...)...) + ) end # Compute the degrees of freedom diff --git a/lib/DataDrivenLux/src/lux/simplex.jl b/lib/DataDrivenLux/src/lux/simplex.jl index cb13aa1bc..a4ca8b258 100644 --- a/lib/DataDrivenLux/src/lux/simplex.jl +++ b/lib/DataDrivenLux/src/lux/simplex.jl @@ -9,7 +9,8 @@ on each row. struct Softmax <: AbstractSimplex end function (::Softmax)( - rng::AbstractRNG, x̂::AbstractVector, x::AbstractVector, κ = one(eltype(x))) + rng::AbstractRNG, x̂::AbstractVector, x::AbstractVector, κ = one(eltype(x)) + ) return softmax!(x̂, x ./ κ) end @@ -26,7 +27,8 @@ $(FIELDS) struct GumbelSoftmax <: AbstractSimplex end function (::GumbelSoftmax)( - rng::AbstractRNG, x̂::AbstractVector, x::AbstractVector, κ = one(eltype(x))) + rng::AbstractRNG, x̂::AbstractVector, x::AbstractVector, κ = one(eltype(x)) + ) z = -log.(-log.(rand(rng, size(x)...))) y = similar(x) foreach(axes(x, 2)) do i @@ -48,7 +50,8 @@ $(FIELDS) struct DirectSimplex <: AbstractSimplex end function (::DirectSimplex)( - rng::AbstractRNG, x̂::AbstractVector, x::AbstractVector, κ = one(eltype(x))) + rng::AbstractRNG, x̂::AbstractVector, x::AbstractVector, κ = one(eltype(x)) + ) return x̂ .= x end diff --git a/lib/DataDrivenLux/src/solve.jl b/lib/DataDrivenLux/src/solve.jl index 29826dc16..418a67847 100644 --- a/lib/DataDrivenLux/src/solve.jl +++ b/lib/DataDrivenLux/src/solve.jl @@ -1,5 +1,7 @@ -function DataDrivenDiffEq.get_fit_targets(::A, prob::AbstractDataDrivenProblem, - basis::Basis) where {A <: AbstractDAGSRAlgorithm} +function DataDrivenDiffEq.get_fit_targets( + ::A, prob::AbstractDataDrivenProblem, + basis::Basis + ) where {A <: AbstractDAGSRAlgorithm} return prob.X, DataDrivenDiffEq.get_implicit_data(prob) end @@ -8,8 +10,10 @@ end retcode <: DDReturnCode end -function CommonSolve.solve!(prob::InternalDataDrivenProblem{A}) where {A <: - AbstractDAGSRAlgorithm} +function CommonSolve.solve!(prob::InternalDataDrivenProblem{A}) where { + A <: + AbstractDAGSRAlgorithm, + } (; alg, basis, testdata, traindata, control_idx, options, problem, kwargs) = prob (; maxiters, progress, eval_expresssion, abstol) = options @@ -21,14 +25,18 @@ function CommonSolve.solve!(prob::InternalDataDrivenProblem{A}) where {A <: shows = min(5, sum(cache.keeps)) losses = map(alg.options.loss, cache.candidates[cache.keeps]) min_, max_ = extrema(losses) - [(:Iterations, iter), + [ + (:Iterations, iter), (:RSS, map(StatsBase.rss, cache.candidates[cache.keeps][1:shows])), (:Minimum, min_), (:Maximum, max_), (:Mode, mode(losses)), (:Mean, mean(losses)), - (:Probabilities, - map(x -> exp.(x(cache.p)), cache.candidates[cache.keeps][1:shows]))] + ( + :Probabilities, + map(x -> exp.(x(cache.p)), cache.candidates[cache.keeps][1:shows]), + ), + ] end end @@ -52,8 +60,10 @@ function CommonSolve.solve!(prob::InternalDataDrivenProblem{A}) where {A <: new_problem = DataDrivenDiffEq.remake_problem(problem, p = pnew) rss = sum( - abs2, new_basis(new_problem) .- DataDrivenDiffEq.get_implicit_data(new_problem)) + abs2, new_basis(new_problem) .- DataDrivenDiffEq.get_implicit_data(new_problem) + ) return DataDrivenSolution{typeof(rss)}( - new_basis, DDReturnCode(1), alg, [cache], new_problem, rss, length(pnew), prob) + new_basis, DDReturnCode(1), alg, [cache], new_problem, rss, length(pnew), prob + ) end diff --git a/lib/DataDrivenLux/test/cache.jl b/lib/DataDrivenLux/test/cache.jl index fa9f1dda9..4624fd962 100644 --- a/lib/DataDrivenLux/test/cache.jl +++ b/lib/DataDrivenLux/test/cache.jl @@ -14,9 +14,11 @@ Y = sin.(X[1:1, :]) dummy_basis = Basis(x, x) dummy_problem = DirectDataDrivenProblem(X, Y) -# We have 1 Choices in the first layer, 2 in the last -alg = RandomSearch(populationsize = 10, functions = (sin,), arities = (1,), - rng = rng, loss = rss, keep = 1, distributed = false) +# We have 1 Choices in the first layer, 2 in the last +alg = RandomSearch( + populationsize = 10, functions = (sin,), arities = (1,), + rng = rng, loss = rss, keep = 1, distributed = false +) cache = DataDrivenLux.init_cache(alg, dummy_basis, dummy_problem) rss_wrong = sum(abs2, Y .- X) diff --git a/lib/DataDrivenLux/test/crossentropy_solve.jl b/lib/DataDrivenLux/test/crossentropy_solve.jl index a66f317d1..166c9ab0e 100644 --- a/lib/DataDrivenLux/test/crossentropy_solve.jl +++ b/lib/DataDrivenLux/test/crossentropy_solve.jl @@ -33,15 +33,20 @@ dummy_dataset = DataDrivenLux.Dataset(dummy_problem) @parameters p [bounds = (-3.0, -1.0), dist = truncated(Normal(-2.0, 1.0), -3.0, -1.0)] b = Basis([x; exp.(x)], x) -# We have 1 Choices in the first layer, 2 in the last -alg = CrossEntropy(populationsize = 2_00, functions = (sin, exp, +), arities = (1, 1, 2), +# We have 1 Choices in the first layer, 2 in the last +alg = CrossEntropy( + populationsize = 2_00, functions = (sin, exp, +), arities = (1, 1, 2), rng = rng, n_layers = 3, use_protected = true, loss = bic, keep = 0.1, - threaded = true, optim_options = Optim.Options(time_limit = 0.2)) + threaded = true, optim_options = Optim.Options(time_limit = 0.2) +) -res = solve(dummy_problem, b, alg, +res = solve( + dummy_problem, b, alg, options = DataDrivenCommonOptions( - maxiters = 1_000, progress = parse(Bool, get(ENV, "CI", "false")), abstol = 0.0)) -@test rss(res) <= 1e-2 + maxiters = 1_000, progress = parse(Bool, get(ENV, "CI", "false")), abstol = 0.0 + ) +) +@test rss(res) <= 1.0e-2 @test aicc(res) <= -100.0 @test r2(res) >= 0.95 results = get_results(res) diff --git a/lib/DataDrivenLux/test/graphs.jl b/lib/DataDrivenLux/test/graphs.jl index 6bbbca029..2101cb5bf 100644 --- a/lib/DataDrivenLux/test/graphs.jl +++ b/lib/DataDrivenLux/test/graphs.jl @@ -23,8 +23,11 @@ X = randn(1, 10) Y, _ = dag(X, ps, new_st) @test y == [sin.(x[1]); sin.(x[1])] @test Y == [sin.(X[1:1, :]); sin.(X[1:1, :])] - @test exp(sum( - sum ∘ values, values(DataDrivenLux.get_loglikelihood(dag, ps, new_st)))) == 1.0f0 + @test exp( + sum( + sum ∘ values, values(DataDrivenLux.get_loglikelihood(dag, ps, new_st)) + ) + ) == 1.0f0 end @testset "Two Layer Skip" begin @@ -38,5 +41,5 @@ end @test y == [sin.(x[1]) .+ x[1]; x[1]] @test Y == [sin.(X[1:1, :]) .+ X[1:1, :]; X[1:1, :]] @test DataDrivenLux.get_loglikelihood(dag, ps, new_st, out_state) == - sum(Float32[-2.7725887, -1.3862944]) + sum(Float32[-2.7725887, -1.3862944]) end diff --git a/lib/DataDrivenLux/test/layers.jl b/lib/DataDrivenLux/test/layers.jl index 46173997b..313ec683e 100644 --- a/lib/DataDrivenLux/test/layers.jl +++ b/lib/DataDrivenLux/test/layers.jl @@ -18,8 +18,10 @@ using StableRNGs rng = StableRNG(43) ps, st = Lux.setup(rng, layer) layer_states, new_st = layer(states, ps, st) - @test all(exp.(values(DataDrivenLux.get_loglikelihood(layer, ps, new_st))) .≈ - (1 / 3, 1 / 9, 1 / 27)) + @test all( + exp.(values(DataDrivenLux.get_loglikelihood(layer, ps, new_st))) .≈ + (1 / 3, 1 / 9, 1 / 27) + ) intervals = map(DataDrivenLux.get_interval, layer_states) @test isequal_interval(intervals[1], interval(-1, 1)) diff --git a/lib/DataDrivenLux/test/randomsearch_solve.jl b/lib/DataDrivenLux/test/randomsearch_solve.jl index 71238c0ae..518f7bac6 100644 --- a/lib/DataDrivenLux/test/randomsearch_solve.jl +++ b/lib/DataDrivenLux/test/randomsearch_solve.jl @@ -28,20 +28,26 @@ dummy_dataset = DataDrivenLux.Dataset(dummy_problem) @test isempty(dummy_dataset.u_intervals) -for (data, _interval) in zip((X, Y, 1:size(X, 2)), - (dummy_dataset.x_intervals[1], dummy_dataset.y_intervals[1], dummy_dataset.t_interval)) +for (data, _interval) in zip( + (X, Y, 1:size(X, 2)), + (dummy_dataset.x_intervals[1], dummy_dataset.y_intervals[1], dummy_dataset.t_interval) + ) @test isequal_interval(_interval, interval(extrema(data))) end -# We have 1 Choices in the first layer, 2 in the last +# We have 1 Choices in the first layer, 2 in the last alg = RandomSearch(; populationsize = 10, functions = (sin, exp, *), arities = (1, 1, 2), rng, - n_layers = 2, loss = rss, keep = 2) + n_layers = 2, loss = rss, keep = 2 +) -res = solve(dummy_problem, alg, +res = solve( + dummy_problem, alg, options = DataDrivenCommonOptions( - maxiters = 50, progress = parse(Bool, get(ENV, "CI", "false")), abstol = 0.0)) -@test rss(res) <= 1e-2 + maxiters = 50, progress = parse(Bool, get(ENV, "CI", "false")), abstol = 0.0 + ) +) +@test rss(res) <= 1.0e-2 @test aicc(res) <= -100.0 @test r2(res) >= 0.95 results = get_results(res) diff --git a/lib/DataDrivenLux/test/reinforce_solve.jl b/lib/DataDrivenLux/test/reinforce_solve.jl index 9fb79e9e6..bf7b348ab 100644 --- a/lib/DataDrivenLux/test/reinforce_solve.jl +++ b/lib/DataDrivenLux/test/reinforce_solve.jl @@ -33,17 +33,21 @@ dummy_dataset = DataDrivenLux.Dataset(dummy_problem) @parameters p [bounds = (-3.0, -1.0), dist = truncated(Normal(-2.0, 1.0), -3.0, -1.0)] b = Basis([x; exp.(x)], x) -# We have 1 Choices in the first layer, 2 in the last +# We have 1 Choices in the first layer, 2 in the last alg = Reinforce(; populationsize = 200, functions = (sin, exp, +), arities = (1, 1, 2), rng, n_layers = 3, use_protected = true, loss = bic, keep = 10, threaded = true, - optim_options = Optim.Options(time_limit = 0.2), optimiser = AdamW(1e-2)) + optim_options = Optim.Options(time_limit = 0.2), optimiser = AdamW(1.0e-2) +) -res = solve(dummy_problem, b, alg, +res = solve( + dummy_problem, b, alg, options = DataDrivenCommonOptions( - maxiters = 1000, progress = parse(Bool, get(ENV, "CI", "false")), abstol = 0.0)) + maxiters = 1000, progress = parse(Bool, get(ENV, "CI", "false")), abstol = 0.0 + ) +) -@test rss(res) <= 1e-2 +@test rss(res) <= 1.0e-2 @test aicc(res) <= -100.0 @test r2(res) >= 0.95 results = get_results(res) diff --git a/lib/DataDrivenSR/src/DataDrivenSR.jl b/lib/DataDrivenSR/src/DataDrivenSR.jl index 7f5f6a35f..40d94e7a9 100644 --- a/lib/DataDrivenSR/src/DataDrivenSR.jl +++ b/lib/DataDrivenSR/src/DataDrivenSR.jl @@ -79,9 +79,11 @@ function SRResult(prob, hof, paretos) nobs = prod(size(y)) ll = iszero(rss) ? convert(eltype(rss), Inf) : -nobs / 2 * log(rss / nobs) ll0 = -nobs / 2 * log.(sum(abs2, y .- mean(y, dims = 2)[:, 1]) / nobs) - return SRResult(bs, hof, paretos, + return SRResult( + bs, hof, paretos, rss, ll, ll0, dof, nobs, - DDReturnCode(1)) + DDReturnCode(1) + ) end is_success(k::SRResult) = getfield(k, :retcode) == DDReturnCode(1) @@ -134,7 +136,7 @@ function convert_to_basis(paretofrontier, prob) @unpack alg, basis, problem, options = prob @unpack eq_options = alg @unpack maxiters, eval_expresssion, generate_symbolic_parameters, digits, - roundingmode = options + roundingmode = options eqs_ = map(paretofrontier) do dom node_to_symbolic(dom[end].tree, eq_options) @@ -143,8 +145,12 @@ function convert_to_basis(paretofrontier, prob) # Substitute with the basis elements atoms = map(xi -> xi.rhs, equations(basis)) - subs = Dict([SymbolicUtils.Sym{LiteralReal}(Symbol("x$(i)")) => x - for (i, x) in enumerate(atoms)]...) + subs = Dict( + [ + SymbolicUtils.Sym{LiteralReal}(Symbol("x$(i)")) => x + for (i, x) in enumerate(atoms) + ]... + ) eqs, ps = collect_numerical_parameters(eqs_) eqs = map(Base.Fix2(substitute, subs), eqs) @@ -173,12 +179,14 @@ function convert_to_basis(paretofrontier, prob) DataDrivenDiffEq._set_default_val(Num(ps_[i]), p[i]) end - Basis(eqs, states(basis), + return Basis( + eqs, states(basis), parameters = [p_new; ps], iv = get_iv(basis), controls = controls(basis), observed = observed(basis), implicits = implicit_variables(basis), name = gensym(:Basis), - eval_expression = eval_expresssion) + eval_expression = eval_expresssion + ) end # apply the algorithm on each dataset @@ -187,13 +195,15 @@ function (x::EQSearch)(ps::InternalDataDrivenProblem{EQSearch}, X, Y) @unpack maxiters, abstol = options @unpack weights, eq_options, numprocs, procs, parallelism, runtests = x - hofs = SymbolicRegression.equation_search(X, Y; + hofs = SymbolicRegression.equation_search( + X, Y; niterations = maxiters, weights = weights, options = eq_options, numprocs = numprocs, procs = procs, parallelism = parallelism, - runtests = runtests) + runtests = runtests + ) # We always want something which is a vector or tuple hofs = !isa(hofs, AbstractVector) ? [hofs] : hofs @@ -209,10 +219,10 @@ end function CommonSolve.solve!(ps::InternalDataDrivenProblem{EQSearch}) @unpack alg, basis, testdata, traindata, kwargs = ps @unpack weights, numprocs, procs, addprocs_function, parallelism, runtests, - eq_options = alg + eq_options = alg @unpack traindata, testdata, basis, options = ps @unpack maxiters, eval_expresssion, generate_symbolic_parameters, - digits, roundingmode, selector = options + digits, roundingmode, selector = options @unpack problem = ps results = map(traindata) do (X, Y) @@ -222,7 +232,7 @@ function CommonSolve.solve!(ps::InternalDataDrivenProblem{EQSearch}) idx = argmin(map(selector, results)) best_res = results[idx] - DataDrivenSolution(best_res.basis, problem, alg, results, ps, best_res.retcode) + return DataDrivenSolution(best_res.basis, problem, alg, results, ps, best_res.retcode) end export EQSearch diff --git a/lib/DataDrivenSR/test/runtests.jl b/lib/DataDrivenSR/test/runtests.jl index 320ec990e..be923bb8b 100644 --- a/lib/DataDrivenSR/test/runtests.jl +++ b/lib/DataDrivenSR/test/runtests.jl @@ -7,36 +7,48 @@ rng = StableRNG(42) X = rand(rng, 2, 50) @testset "Simple" begin - alg = DataDrivenSR.EQSearch(eq_options = Options(unary_operators = [sin, exp], - binary_operators = [*], maxdepth = 1, - seed = 42, - verbosity = -1, progress = false)) + alg = DataDrivenSR.EQSearch( + eq_options = Options( + unary_operators = [sin, exp], + binary_operators = [*], maxdepth = 1, + seed = 42, + verbosity = -1, progress = false + ) + ) f(x) = [sin(x[1]); exp(x[2])] Y = hcat(map(f, eachcol(X))...) prob = DirectDataDrivenProblem(X, Y) res = solve(prob, alg) @test r2(res) >= 0.95 - @test rss(res) <= 1e-5 + @test rss(res) <= 1.0e-5 end @testset "Univariate" begin - alg = DataDrivenSR.EQSearch(eq_options = Options(unary_operators = [sin, exp], - binary_operators = [*], maxdepth = 1, - seed = 42, - verbosity = -1, progress = false)) + alg = DataDrivenSR.EQSearch( + eq_options = Options( + unary_operators = [sin, exp], + binary_operators = [*], maxdepth = 1, + seed = 42, + verbosity = -1, progress = false + ) + ) Y = sin.(X[1:1, :]) prob = DirectDataDrivenProblem(X, Y) res = solve(prob, alg) @test r2(res) >= 0.95 - @test rss(res) <= 1e-5 + @test rss(res) <= 1.0e-5 end @testset "Lifted" begin - alg = DataDrivenSR.EQSearch(eq_options = Options(unary_operators = [sin, exp], - binary_operators = [+], maxdepth = 1, - seed = 42, - verbosity = -1, progress = false)) + alg = DataDrivenSR.EQSearch( + eq_options = Options( + unary_operators = [sin, exp], + binary_operators = [+], maxdepth = 1, + seed = 42, + verbosity = -1, progress = false + ) + ) f(x) = [sin(x[1] .^ 2); exp(x[2] * x[1])] Y = hcat(map(f, eachcol(X))...) @@ -46,5 +58,5 @@ end prob = DirectDataDrivenProblem(X, Y) res = solve(prob, basis, alg) @test r2(res) >= 0.95 - @test rss(res) <= 1e-5 + @test rss(res) <= 1.0e-5 end diff --git a/lib/DataDrivenSparse/src/DataDrivenSparse.jl b/lib/DataDrivenSparse/src/DataDrivenSparse.jl index 5fc7f7351..22d9ee13e 100644 --- a/lib/DataDrivenSparse/src/DataDrivenSparse.jl +++ b/lib/DataDrivenSparse/src/DataDrivenSparse.jl @@ -70,13 +70,13 @@ StatsBase.nobs(x::AbstractSparseRegressionCache) = begin end function StatsBase.loglikelihood(x::AbstractSparseRegressionCache) - begin + return begin -nobs(x) / 2 * log(rss(x) / nobs(x)) end end function StatsBase.nullloglikelihood(x::AbstractSparseRegressionCache) - begin + return begin @unpack B̃ = x -nobs(x) / 2 * log(mean(abs2, B̃ .- mean(vec(B̃)))) end @@ -96,9 +96,11 @@ get_proximal(x::AbstractSparseRegressionAlgorithm) = SoftThreshold() include("solver.jl") export SparseLinearSolver -function (x::X where {X <: AbstractSparseRegressionAlgorithm})(X, Y; +function (x::X where {X <: AbstractSparseRegressionAlgorithm})( + X, Y; options::DataDrivenCommonOptions = DataDrivenCommonOptions(), - kwargs...) + kwargs... + ) solver = SparseLinearSolver(x, options = options) results = solver(X, Y) # Keep this here for now diff --git a/lib/DataDrivenSparse/src/algorithms/ADMM.jl b/lib/DataDrivenSparse/src/algorithms/ADMM.jl index dbffda85f..8183692bb 100644 --- a/lib/DataDrivenSparse/src/algorithms/ADMM.jl +++ b/lib/DataDrivenSparse/src/algorithms/ADMM.jl @@ -25,9 +25,9 @@ mutable struct ADMM{T, R <: Number} <: AbstractSparseRegressionAlgorithm """Augmented Lagrangian parameter""" rho::R - function ADMM(threshold::T = 1e-1, ρ::R = 1.0) where {T, R} + function ADMM(threshold::T = 1.0e-1, ρ::R = 1.0) where {T, R} @assert all(threshold .> zero(eltype(threshold))) "Threshold must be positive definite" - @assert zero(R)<ρ "Augmented Lagrangian parameter should be positive definite" + @assert zero(R) < ρ "Augmented Lagrangian parameter should be positive definite" return new{T, R}(threshold, ρ) end end @@ -52,13 +52,13 @@ struct ADMMCache{fat, C, A, AT, BT, T, ATT, BTT} <: AbstractSparseRegressionCach end function init_cache(alg::ADMM, A::AbstractMatrix, b::AbstractVector) - init_cache(alg, A, permutedims(b)) + return init_cache(alg, A, permutedims(b)) end function init_cache(alg::ADMM, A::AbstractMatrix, B::AbstractMatrix) n_x, m_x = size(A) - @assert size(B, 1)==1 "Caches only hold single targets!" + @assert size(B, 1) == 1 "Caches only hold single targets!" λ = minimum(get_thresholds(alg)) @@ -81,12 +81,16 @@ function init_cache(alg::ADMM, A::AbstractMatrix, B::AbstractMatrix) active_set!(idx, proximal, coefficients, λ / rho) - return ADMMCache{fat, typeof(coefficients), typeof(idx), typeof(X), typeof(Y), - typeof(rho), typeof(A), typeof(B)}(coefficients, zero(coefficients), + return ADMMCache{ + fat, typeof(coefficients), typeof(idx), typeof(X), typeof(Y), + typeof(rho), typeof(A), typeof(B), + }( + coefficients, zero(coefficients), idx, proximal, zero(coefficients), zero(coefficients), - X, Y, rho, A, B) + X, Y, rho, A, B + ) end # Fat regression diff --git a/lib/DataDrivenSparse/src/algorithms/Implicit.jl b/lib/DataDrivenSparse/src/algorithms/Implicit.jl index 3f72605a4..0e6848d81 100644 --- a/lib/DataDrivenSparse/src/algorithms/Implicit.jl +++ b/lib/DataDrivenSparse/src/algorithms/Implicit.jl @@ -19,11 +19,11 @@ ImplicitOptimizer(0.1f0, ADMM) ``` """ mutable struct ImplicitOptimizer{T <: AbstractSparseRegressionAlgorithm} <: - AbstractSparseRegressionAlgorithm + AbstractSparseRegressionAlgorithm """Explicit Optimizer""" optimizer::T - function ImplicitOptimizer(threshold = 1e-1, opt = STLSQ) + function ImplicitOptimizer(threshold = 1.0e-1, opt = STLSQ) optimizer = opt(threshold) return new{typeof(optimizer)}(optimizer) end @@ -37,10 +37,12 @@ Base.summary(opt::ImplicitOptimizer) = "Implicit Optimizer using " * summary(opt get_threshold(opt::ImplicitOptimizer) = get_threshold(opt.optimizer) -function (x::ImplicitOptimizer)(X, Y; +function (x::ImplicitOptimizer)( + X, Y; options::DataDrivenCommonOptions = DataDrivenCommonOptions(), necessary_idx = ones(Bool, size(X, 1)), - kwargs...) + kwargs... + ) @unpack optimizer = x @unpack verbose = options diff --git a/lib/DataDrivenSparse/src/algorithms/SR3.jl b/lib/DataDrivenSparse/src/algorithms/SR3.jl index 191b3c765..9c99ea51d 100644 --- a/lib/DataDrivenSparse/src/algorithms/SR3.jl +++ b/lib/DataDrivenSparse/src/algorithms/SR3.jl @@ -39,11 +39,15 @@ mutable struct SR3{T, V, P <: AbstractProximalOperator} <: AbstractSparseRegress """Proximal operator""" proximal::P - function SR3(threshold::T = 1e-1, nu::V = 1.0, - R::P = HardThreshold()) where {T, V <: Number, - P <: AbstractProximalOperator} + function SR3( + threshold::T = 1.0e-1, nu::V = 1.0, + R::P = HardThreshold() + ) where { + T, V <: Number, + P <: AbstractProximalOperator, + } @assert all(threshold .> zero(eltype(threshold))) "Threshold must be positive definite" - @assert nu>zero(V) "Relaxation must be positive definite" + @assert nu > zero(V) "Relaxation must be positive definite" λ = isa(R, HardThreshold) ? threshold .^ 2 / 2 : threshold return new{typeof(λ), V, P}(λ, nu, R) @@ -60,7 +64,7 @@ end Base.summary(::SR3) = "SR3" struct SR3Cache{C, A, P <: AbstractProximalOperator, AT, BT, T, ATT, BTT} <: - AbstractSparseRegressionCache + AbstractSparseRegressionCache X::C X_prev::C active_set::A @@ -77,13 +81,13 @@ struct SR3Cache{C, A, P <: AbstractProximalOperator, AT, BT, T, ATT, BTT} <: end function init_cache(alg::SR3, A::AbstractMatrix, b::AbstractVector) - init_cache(alg, A, permutedims(b)) + return init_cache(alg, A, permutedims(b)) end function init_cache(alg::SR3, A::AbstractMatrix, B::AbstractMatrix) n_x, m_x = size(A) - @assert size(B, 1)==1 "Caches only hold single targets!" + @assert size(B, 1) == 1 "Caches only hold single targets!" λ = minimum(get_thresholds(alg)) @@ -99,12 +103,16 @@ function init_cache(alg::SR3, A::AbstractMatrix, B::AbstractMatrix) active_set!(idx, proximal, coefficients, λ) - return SR3Cache{typeof(coefficients), typeof(idx), typeof(proximal), typeof(X), - typeof(Y), typeof(nu), typeof(A), typeof(B)}(coefficients, + return SR3Cache{ + typeof(coefficients), typeof(idx), typeof(proximal), typeof(X), + typeof(Y), typeof(nu), typeof(A), typeof(B), + }( + coefficients, copy(coefficients), idx, proximal, zero(coefficients), - X, Y, nu, A, B) + X, Y, nu, A, B + ) end function step!(cache::SR3Cache, λ::T) where {T <: Number} diff --git a/lib/DataDrivenSparse/src/algorithms/STLSQ.jl b/lib/DataDrivenSparse/src/algorithms/STLSQ.jl index 0a5d02a13..9df849803 100644 --- a/lib/DataDrivenSparse/src/algorithms/STLSQ.jl +++ b/lib/DataDrivenSparse/src/algorithms/STLSQ.jl @@ -36,15 +36,15 @@ opt = STLSQ(Float32[1e-2; 1e-1]) This was formally `STRRidge` and has been renamed. """ struct STLSQ{T <: Union{Number, AbstractVector}, R <: Number} <: - AbstractSparseRegressionAlgorithm + AbstractSparseRegressionAlgorithm """Sparsity threshold""" thresholds::T """Ridge regression parameter""" rho::R - function STLSQ(threshold::T = 1e-1, rho::R = zero(eltype(T))) where {T, R <: Number} + function STLSQ(threshold::T = 1.0e-1, rho::R = zero(eltype(T))) where {T, R <: Number} @assert all(threshold .> zero(eltype(threshold))) "Threshold must be positive definite" - @assert rho>=zero(R) "Ridge regression parameter must be positive definite!" + @assert rho >= zero(R) "Ridge regression parameter must be positive definite!" return new{T, R}(threshold, rho) end end @@ -52,7 +52,7 @@ end Base.summary(::STLSQ) = "STLSQ" struct STLSQCache{usenormal, C <: AbstractArray, A <: BitArray, AT, BT, ATT, BTT} <: - AbstractSparseRegressionCache + AbstractSparseRegressionCache X::C X_prev::C active_set::A @@ -65,12 +65,12 @@ struct STLSQCache{usenormal, C <: AbstractArray, A <: BitArray, AT, BT, ATT, BTT end function init_cache(alg::STLSQ, A::AbstractMatrix, b::AbstractVector) - init_cache(alg, A, permutedims(b)) + return init_cache(alg, A, permutedims(b)) end function init_cache(alg::STLSQ, A::AbstractMatrix, B::AbstractMatrix) n_x, m_x = size(A) - @assert size(B, 1)==1 "Caches only hold single targets!" + @assert size(B, 1) == 1 "Caches only hold single targets!" @unpack rho = alg λ = minimum(get_thresholds(alg)) @@ -94,10 +94,14 @@ function init_cache(alg::STLSQ, A::AbstractMatrix, B::AbstractMatrix) active_set!(active_set, proximal, coefficients, λ) - return STLSQCache{usenormal, typeof(coefficients), typeof(active_set), typeof(X), - typeof(Y), typeof(A), typeof(B)}(coefficients, prev_coefficients, + return STLSQCache{ + usenormal, typeof(coefficients), typeof(active_set), typeof(X), + typeof(Y), typeof(A), typeof(B), + }( + coefficients, prev_coefficients, active_set, get_proximal(alg), - X, Y, A, B) + X, Y, A, B + ) end function step!(cache::STLSQCache, λ::T) where {T} diff --git a/lib/DataDrivenSparse/src/algorithms/proximals.jl b/lib/DataDrivenSparse/src/algorithms/proximals.jl index 5d6ea5ca8..f3099324a 100644 --- a/lib/DataDrivenSparse/src/algorithms/proximals.jl +++ b/lib/DataDrivenSparse/src/algorithms/proximals.jl @@ -1,6 +1,8 @@ -@inline function (s::AbstractProximalOperator)(x::AbstractArray{<:Number}, +@inline function (s::AbstractProximalOperator)( + x::AbstractArray{<:Number}, y::AbstractArray{<:Bool}, - λ::T) where {T <: Real} + λ::T + ) where {T <: Real} @assert size(y) == size(x) active_set!(y, s, x, λ) for i in eachindex(x) @@ -21,8 +23,10 @@ See [by Zheng et al., 2018](https://ieeexplore.ieee.org/document/8573778). """ struct SoftThreshold <: AbstractProximalOperator end; -@inline function active_set!(idx::BitArray, ::SoftThreshold, x::AbstractArray{T}, - λ::T) where {T} +@inline function active_set!( + idx::BitArray, ::SoftThreshold, x::AbstractArray{T}, + λ::T + ) where {T} @assert size(idx) == size(x) @inbounds foreach(eachindex(x)) do i idx[i] = abs(x[i]) > λ @@ -37,8 +41,10 @@ end return end -@inline function (s::SoftThreshold)(y::AbstractArray{N}, x::AbstractArray{N}, - λ::T) where {N <: Number, T <: Real} +@inline function (s::SoftThreshold)( + y::AbstractArray{N}, x::AbstractArray{N}, + λ::T + ) where {N <: Number, T <: Real} @assert size(y) == size(x) for i in eachindex(x) y[i] = sign(x[i]) * max(abs(x[i]) - λ, zero(eltype(x))) @@ -58,8 +64,10 @@ See [by Zheng et al., 2018](https://ieeexplore.ieee.org/document/8573778). """ struct HardThreshold <: AbstractProximalOperator end; -@inline function active_set!(idx::BitArray, ::HardThreshold, x::AbstractArray, - λ::T) where {T} +@inline function active_set!( + idx::BitArray, ::HardThreshold, x::AbstractArray, + λ::T + ) where {T} @assert size(idx) == size(x) @inbounds foreach(eachindex(x)) do i idx[i] = abs(x[i]) > sqrt(2 * λ) @@ -74,8 +82,10 @@ end return end -@inline function (s::HardThreshold)(y::AbstractArray{N}, x::AbstractArray{N}, - λ::T) where {N <: Number, T <: Real} +@inline function (s::HardThreshold)( + y::AbstractArray{N}, x::AbstractArray{N}, + λ::T + ) where {N <: Number, T <: Real} @assert all(size(y) .== size(x)) for i in eachindex(x) y[i] = abs(x[i]) > sqrt(2 * λ) ? x[i] : zero(eltype(x)) @@ -112,8 +122,10 @@ end ClippedAbsoluteDeviation() = ClippedAbsoluteDeviation(NaN) -@inline function active_set!(idx::BitArray, h::ClippedAbsoluteDeviation, - x::AbstractArray, λ::T) where {T} +@inline function active_set!( + idx::BitArray, h::ClippedAbsoluteDeviation, + x::AbstractArray, λ::T + ) where {T} @assert size(idx) == size(x) @unpack ρ = h ρ = isnan(ρ) ? convert(T, 5) * λ : convert(T, ρ) @@ -132,8 +144,10 @@ function (s::ClippedAbsoluteDeviation)(x::AbstractArray, λ::T) where {T <: Real return end -function (s::ClippedAbsoluteDeviation)(y::AbstractArray{N}, x::AbstractArray{N}, - λ::T) where {N <: Number, T <: Real} +function (s::ClippedAbsoluteDeviation)( + y::AbstractArray{N}, x::AbstractArray{N}, + λ::T + ) where {N <: Number, T <: Real} @assert all(size(y) .== size(x)) ρ = isnan(s.ρ) ? convert(eltype(x), 5) * λ : convert(eltype(x), s.ρ) for i in eachindex(x) diff --git a/lib/DataDrivenSparse/src/commonsolve.jl b/lib/DataDrivenSparse/src/commonsolve.jl index b701ddb52..2553bd60d 100644 --- a/lib/DataDrivenSparse/src/commonsolve.jl +++ b/lib/DataDrivenSparse/src/commonsolve.jl @@ -1,6 +1,8 @@ -function CommonSolve.solve!(ps::InternalDataDrivenProblem{ - <:AbstractSparseRegressionAlgorithm -}) +function CommonSolve.solve!( + ps::InternalDataDrivenProblem{ + <:AbstractSparseRegressionAlgorithm, + } + ) @unpack alg, basis, testdata, traindata, problem, options, transform = ps results = map(traindata) do (X, Y) @@ -18,15 +20,16 @@ function CommonSolve.solve!(ps::InternalDataDrivenProblem{ coefficients = permutedims(StatsBase.transform(transform, coefficients)) new_basis = DataDrivenDiffEq.__construct_basis(coefficients, basis, problem, options) - DataDrivenSolution(new_basis, problem, alg, results, ps, best_res.retcode) + return DataDrivenSolution(new_basis, problem, alg, results, ps, best_res.retcode) end function __sparse_regression( ps::InternalDataDrivenProblem{ - <:AbstractSparseRegressionAlgorithm + <:AbstractSparseRegressionAlgorithm, }, X::AbstractArray, - Y::AbstractArray) + Y::AbstractArray + ) @unpack alg, testdata, options, transform = ps coefficients, optimal_thresholds, optimal_iterations = alg(X, Y, options = options) @@ -45,13 +48,17 @@ function __sparse_regression( dof = sum(abs.(coefficients) .> 0.0) - SparseRegressionResult(coefficients, dof, optimal_thresholds, + return SparseRegressionResult( + coefficients, dof, optimal_thresholds, optimal_iterations, testerror, trainerror, - retcode) + retcode + ) end -function __sparse_regression(ps::InternalDataDrivenProblem{<:ImplicitOptimizer}, - X::AbstractArray, Y::AbstractArray) +function __sparse_regression( + ps::InternalDataDrivenProblem{<:ImplicitOptimizer}, + X::AbstractArray, Y::AbstractArray + ) @unpack alg, testdata, options, transform, basis, problem, implicit_idx = ps @assert DataDrivenDiffEq.is_implicit(basis) "The provided `Basis` does not have implicit variables!" @@ -66,16 +73,20 @@ function __sparse_regression(ps::InternalDataDrivenProblem{<:ImplicitOptimizer}, candidate_matrix[i, j] = implicit_idx[i, j] || sum(implicit_idx[i, idx]) == 0 end - opt_coefficients = zeros(eltype(problem), size(candidate_matrix, 2), - size(candidate_matrix, 1)) + opt_coefficients = zeros( + eltype(problem), size(candidate_matrix, 2), + size(candidate_matrix, 1) + ) opt_thresholds = [] opt_iterations = [] foreach(enumerate(eachcol(candidate_matrix))) do (i, idx) # We enforce that one of the implicit variables is necessary for success coeff, thresholds, - iters = alg(X[idx, :], Y, options = options, - necessary_idx = implicit_idx[idx, i]) + iters = alg( + X[idx, :], Y, options = options, + necessary_idx = implicit_idx[idx, i] + ) opt_coefficients[i:i, idx] .= coeff push!(opt_thresholds, thresholds) push!(opt_iterations, iters) @@ -95,7 +106,9 @@ function __sparse_regression(ps::InternalDataDrivenProblem{<:ImplicitOptimizer}, dof = sum(abs.(opt_coefficients) .> 0.0) - SparseRegressionResult(opt_coefficients, dof, opt_thresholds, + return SparseRegressionResult( + opt_coefficients, dof, opt_thresholds, opt_iterations, testerror, trainerror, - retcode) + retcode + ) end diff --git a/lib/DataDrivenSparse/src/result.jl b/lib/DataDrivenSparse/src/result.jl index 92f1f4103..b6d6a779e 100644 --- a/lib/DataDrivenSparse/src/result.jl +++ b/lib/DataDrivenSparse/src/result.jl @@ -1,5 +1,5 @@ struct SparseRegressionResult{X <: AbstractArray, L, IT, T, TE, R} <: - AbstractDataDrivenResult + AbstractDataDrivenResult "Coefficient matrix" coefficients::X "Number of nonzeros coefficients" @@ -20,7 +20,7 @@ is_success(k::SparseRegressionResult) = getfield(k, :retcode) == DDReturnCode(1) l2error(k::SparseRegressionResult) = is_success(k) ? getfield(k, :testerror) : Inf function l2error(k::SparseRegressionResult{<:Any, <:Any, Nothing}) - is_success(k) ? getfield(k, :traineerror) : Inf + return is_success(k) ? getfield(k, :traineerror) : Inf end get_coefficients(k::SparseRegressionResult) = getfield(k, :coefficients) diff --git a/lib/DataDrivenSparse/src/solver.jl b/lib/DataDrivenSparse/src/solver.jl index a4c71c0b9..fa6aa9cc0 100644 --- a/lib/DataDrivenSparse/src/solver.jl +++ b/lib/DataDrivenSparse/src/solver.jl @@ -8,21 +8,25 @@ struct SparseLinearSolver{A <: AbstractSparseRegressionAlgorithm, T <: Number} selector::Function end -function SparseLinearSolver(x::A; - options = DataDrivenCommonOptions()) where { +function SparseLinearSolver( + x::A; + options = DataDrivenCommonOptions() + ) where { A <: - AbstractSparseRegressionAlgorithm -} - return SparseLinearSolver(x, + AbstractSparseRegressionAlgorithm, + } + return SparseLinearSolver( + x, options.abstol, options.reltol, options.maxiters, - options.verbose, options.progress, options.selector) + options.verbose, options.progress, options.selector + ) end init_cache(alg::SparseLinearSolver, X, Y) = init_cache(alg.algorithm, X, Y) function (alg::SparseLinearSolver)(X::AbstractMatrix, Y::AbstractMatrix) @unpack verbose = alg - map(axes(Y, 1)) do i + return map(axes(Y, 1)) do i if verbose if i > 1 @printf "\n" diff --git a/lib/DataDrivenSparse/test/cartpole.jl b/lib/DataDrivenSparse/test/cartpole.jl index 4a89ab0b7..92b33e7af 100644 --- a/lib/DataDrivenSparse/test/cartpole.jl +++ b/lib/DataDrivenSparse/test/cartpole.jl @@ -28,8 +28,10 @@ for (i, xi) in enumerate(eachcol(X)) end t = solution.t -ddprob = ContinuousDataDrivenProblem(X, t, DX = DX[3:4, :], - U = (u, p, t) -> [-0.2 + 0.5 * sin(6 * t)]) +ddprob = ContinuousDataDrivenProblem( + X, t, DX = DX[3:4, :], + U = (u, p, t) -> [-0.2 + 0.5 * sin(6 * t)] +) @variables u[1:4] x[1:1] t du = [Symbolics.variable("du", i) for i in 3:4] @@ -54,13 +56,17 @@ push!(implicits, x[1] * sin(u[1])) basis = Basis(implicits, u, controls = x, iv = t, implicits = du) -λ = [1e-4; 5e-4; 1e-3; 2e-3; 3e-3; 4e-3; 5e-3; 6e-3; 7e-3; 8e-3; 9e-3; 1e-2; 2e-2; 3e-2; - 4e-2; 5e-2] +λ = [ + 1.0e-4; 5.0e-4; 1.0e-3; 2.0e-3; 3.0e-3; 4.0e-3; 5.0e-3; 6.0e-3; 7.0e-3; 8.0e-3; 9.0e-3; 1.0e-2; 2.0e-2; 3.0e-2; + 4.0e-2; 5.0e-2 +] -res = solve(ddprob, basis, ImplicitOptimizer(λ), - options = DataDrivenCommonOptions(verbose = false, digits = 3)) +res = solve( + ddprob, basis, ImplicitOptimizer(λ), + options = DataDrivenCommonOptions(verbose = false, digits = 3) +) @test r2(res) >= 0.95 -@test rss(res) <= 1e-2 +@test rss(res) <= 1.0e-2 @test dof(res) == 10 @test get_parameter_values(res.basis) ≈ - [-0.101, 0.05, -1.0, -0.05, -0.05, -0.203, 0.101, -0.101, -1.0, -0.101] + [-0.101, 0.05, -1.0, -0.05, -0.05, -0.203, 0.101, -0.101, -1.0, -0.101] diff --git a/lib/DataDrivenSparse/test/michaelis_menten.jl b/lib/DataDrivenSparse/test/michaelis_menten.jl index 8d2f5df7a..83a7314be 100644 --- a/lib/DataDrivenSparse/test/michaelis_menten.jl +++ b/lib/DataDrivenSparse/test/michaelis_menten.jl @@ -6,7 +6,7 @@ using Test using StatsBase function michaelis_menten(u, p, t) - [0.6 - 1.5u[1] / (0.3 + u[1])] # Should be 0.6*0.3+0.6u[1] - 1.5u[1] = u[2]*u[1]-0.3*u[2] + return [0.6 - 1.5u[1] / (0.3 + u[1])] # Should be 0.6*0.3+0.6u[1] - 1.5u[1] = u[2]*u[1]-0.3*u[2] end u0 = [0.5] @@ -28,30 +28,34 @@ prob = DataDrivenDataset(DataDrivenProblem(solution_1), DataDrivenProblem(soluti @testset "Groundtruth" begin prob = DataDrivenDataset(DataDrivenProblem(solution_1), DataDrivenProblem(solution_2)) - opts = [ImplicitOptimizer(STLSQ(5e-2, 1e-7)); - ImplicitOptimizer(STLSQ(1e-2:1e-2:1e-1, 1e-7))] + opts = [ + ImplicitOptimizer(STLSQ(5.0e-2, 1.0e-7)); + ImplicitOptimizer(STLSQ(1.0e-2:1.0e-2:1.0e-1, 1.0e-7)) + ] for opt in opts res = solve(prob, basis, opt) @test r2(res) >= 0.9 - @test rss(res) < 1e-3 + @test rss(res) < 1.0e-3 @test dof(res) == 4 end end @testset "Noise" begin rng = StableRNG(1111) - prob = DataDrivenDataset(map((solution_1, solution_2)) do sol - X = Array(sol) - X .+= 0.01 * randn(rng, size(X)) - t = sol.t - ContinuousDataDrivenProblem(X, t, GaussianKernel()) - end...) - - opts = [ImplicitOptimizer(ADMM(1e-2:1e-4:1e-1))] + prob = DataDrivenDataset( + map((solution_1, solution_2)) do sol + X = Array(sol) + X .+= 0.01 * randn(rng, size(X)) + t = sol.t + ContinuousDataDrivenProblem(X, t, GaussianKernel()) + end... + ) + + opts = [ImplicitOptimizer(ADMM(1.0e-2:1.0e-4:1.0e-1))] for opt in opts res = solve(prob, basis, opt, options = DataDrivenCommonOptions()) @test r2(res) >= 0.9 - @test rss(res) <= 2e-1 + @test rss(res) <= 2.0e-1 @test dof(res) == 3 end end diff --git a/lib/DataDrivenSparse/test/pendulum.jl b/lib/DataDrivenSparse/test/pendulum.jl index 51f486218..86f991254 100644 --- a/lib/DataDrivenSparse/test/pendulum.jl +++ b/lib/DataDrivenSparse/test/pendulum.jl @@ -25,16 +25,18 @@ sol = solve(prob, Tsit5(), saveat = dt) @testset "Groundtruth" begin dd_prob = DataDrivenProblem(sol) for opt in [ - STLSQ(1e-1), - STLSQ(1e-2:1e-2:1e-1, 0.0001), - ADMM(1e-2), - SR3(1e-2, SoftThreshold()), - SR3(1e-1, ClippedAbsoluteDeviation()), - SR3(5e-1) - ] - res = solve(dd_prob, basis, opt, - options = DataDrivenCommonOptions(maxiters = 10_000, digits = 1)) - @test r2(res)≈0.9 atol=1e-1 + STLSQ(1.0e-1), + STLSQ(1.0e-2:1.0e-2:1.0e-1, 0.0001), + ADMM(1.0e-2), + SR3(1.0e-2, SoftThreshold()), + SR3(1.0e-1, ClippedAbsoluteDeviation()), + SR3(5.0e-1), + ] + res = solve( + dd_prob, basis, opt, + options = DataDrivenCommonOptions(maxiters = 10_000, digits = 1) + ) + @test r2(res) ≈ 0.9 atol = 1.0e-1 @test rss(res) <= 500.0 @test loglikelihood(res) >= 100.0 @test 2 <= dof(res) <= 4 @@ -46,23 +48,26 @@ end t = sol.t rng = StableRNG(21) - X_n = X .+ 1e-1 * randn(rng, size(X)) + X_n = X .+ 1.0e-1 * randn(rng, size(X)) dd_prob = ContinuousDataDrivenProblem(X_n, t, GaussianKernel()) for opt in [ - STLSQ(0.5), - STLSQ(0.5, 0.001), - ADMM(1e-2), - SR3(1e-2, SoftThreshold()), - SR3(1e-1, ClippedAbsoluteDeviation()), - SR3(5e-1) - ] - res = solve(dd_prob, basis, opt, + STLSQ(0.5), + STLSQ(0.5, 0.001), + ADMM(1.0e-2), + SR3(1.0e-2, SoftThreshold()), + SR3(1.0e-1, ClippedAbsoluteDeviation()), + SR3(5.0e-1), + ] + res = solve( + dd_prob, basis, opt, options = DataDrivenCommonOptions( normalize = DataNormalization(ZScoreTransform), denoise = true, - maxiters = 10_000, digits = 1)) - @test r2(res)≈0.9 atol=1e-1 + maxiters = 10_000, digits = 1 + ) + ) + @test r2(res) ≈ 0.9 atol = 1.0e-1 @test rss(res) <= 100.0 @test loglikelihood(res) >= 100.0 @test 2 <= dof(res) <= 4 diff --git a/lib/DataDrivenSparse/test/sparse_linear_solve.jl b/lib/DataDrivenSparse/test/sparse_linear_solve.jl index 6cf5e3b61..dbe073e45 100644 --- a/lib/DataDrivenSparse/test/sparse_linear_solve.jl +++ b/lib/DataDrivenSparse/test/sparse_linear_solve.jl @@ -10,9 +10,15 @@ using StableRNGs rng = StableRNG(42) # Generate data t = 0.0:0.1:10.0 - X = permutedims(reduce(hcat, - (sin.(0.1 .* t), cos.(0.5 .* t), sin.(2.0 .* t .^ 2), - cos.(0.5 .* t .^ 2), exp.(-t)))) + X = permutedims( + reduce( + hcat, + ( + sin.(0.1 .* t), cos.(0.5 .* t), sin.(2.0 .* t .^ 2), + cos.(0.5 .* t .^ 2), exp.(-t), + ) + ) + ) A = [0.68 0.0 0.0 0.0 -1.2] Ỹ = A * X Y = Ỹ + 0.01 * randn(rng, size(Ỹ)) @@ -20,15 +26,19 @@ using StableRNGs true_dof = 2 for alg in [STLSQ, ADMM, SR3] alg_ = alg(LinRange(0.5 * first(λ), 1.5 * last(λ), 20)) - solver = SparseLinearSolver(alg_, - options = DataDrivenCommonOptions(verbose = false, - maxiters = 10_000)) + solver = SparseLinearSolver( + alg_, + options = DataDrivenCommonOptions( + verbose = false, + maxiters = 10_000 + ) + ) res, _... = solver(X, Y) res = first(res) @test rss(res) <= 1.2 @test aicc(res) <= -400.0 @test true_dof == dof(res) - @test r2(res)≈1.0 atol=6e-2 + @test r2(res) ≈ 1.0 atol = 6.0e-2 end end @@ -36,38 +46,54 @@ end rng = StableRNG(52) # Generate data t = 0.0:0.5:2.0 - X = permutedims(reduce(hcat, - (sin.(0.5 .* t), cos.(0.5 .* t), sin.(2.0 .* t .^ 2), - cos.(0.5 .* t .^ 2), exp.(-t), randn(rng, length(t))))) + X = permutedims( + reduce( + hcat, + ( + sin.(0.5 .* t), cos.(0.5 .* t), sin.(2.0 .* t .^ 2), + cos.(0.5 .* t .^ 2), exp.(-t), randn(rng, length(t)), + ) + ) + ) A = [0.68 0.0 0.0 0.0 -1.2 0.0] Y = A * X λ = extrema(abs.(A)[abs.(A) .> 0.0]) true_dof = 2 for alg in [STLSQ, ADMM, SR3] alg_ = alg(LinRange(0.1, 1.6, 15)) - solver = SparseLinearSolver(alg_, - options = DataDrivenCommonOptions(verbose = false, - maxiters = 10_000)) + solver = SparseLinearSolver( + alg_, + options = DataDrivenCommonOptions( + verbose = false, + maxiters = 10_000 + ) + ) res, _... = solver(X, Y) res = first(res) @test rss(res) <= 1.5e-1 @test aicc(res) <= -5.0 @test true_dof == dof(res) - @test r2(res)≈1.0 atol=1e-1 + @test r2(res) ≈ 1.0 atol = 1.0e-1 end end @testset "Implicit Optimizer" begin t = 0.0:0.1:10.0 - X = permutedims(reduce(hcat, - (sin.(0.5 .* t .+ 0.1), cos.(0.5 .* t), sin.(2.0 .* t .^ 2), - cos.(0.5 .* t .^ 2 .- 0.1), exp.(-t)))) + X = permutedims( + reduce( + hcat, + ( + sin.(0.5 .* t .+ 0.1), cos.(0.5 .* t), sin.(2.0 .* t .^ 2), + cos.(0.5 .* t .^ 2 .- 0.1), exp.(-t), + ) + ) + ) Y = permutedims(0.5 * X[1, :] + 0.22 * X[4, :] - 2.0 * X[3, :]) X = vcat(X, Y) for alg in [STLSQ(0.1, 1.0), ADMM(), SR3()] opt = ImplicitOptimizer(alg) rescoeff, _... = opt(X, Y, options = DataDrivenCommonOptions(maxiters = 2000)) - @test vec(rescoeff)≈[0.25; 0.0; -1.0; 0.11; 0.0; -0.5] atol=5e-2 + @test vec(rescoeff) ≈ [0.25; 0.0; -1.0; 0.11; 0.0; -0.5] atol = 5.0e-2 end end @@ -78,15 +104,15 @@ end # Test case 1: Very small data values that lead to zero coefficients after regularization N = 3 - X̂ = randn(rng, N, 50) * 1e-10 - Ŷ = randn(rng, 1, 50) * 1e-10 + X̂ = randn(rng, N, 50) * 1.0e-10 + Ŷ = randn(rng, 1, 50) * 1.0e-10 @variables u[1:N] b = polynomial_basis(u, 2) basis = Basis(b, u) problem = DirectDataDrivenProblem(X̂, Ŷ) - λ = 1e-1 + λ = 1.0e-1 opt = ADMM(λ) options = DataDrivenCommonOptions() @@ -100,7 +126,7 @@ end Ŷ2 = randn(rng, 1, 50) problem2 = DirectDataDrivenProblem(X̂2, Ŷ2) - λ_high = 1e10 # Very high regularization + λ_high = 1.0e10 # Very high regularization opt_high = ADMM(λ_high) result2 = @test_nowarn solve(problem2, basis, opt_high, options = options) diff --git a/src/DataDrivenDiffEq.jl b/src/DataDrivenDiffEq.jl index 917192575..36b786b37 100644 --- a/src/DataDrivenDiffEq.jl +++ b/src/DataDrivenDiffEq.jl @@ -19,7 +19,7 @@ using SymbolicUtils: operation, arguments, iscall, issym using Symbolics using Symbolics: scalarize, variable, value @reexport using ModelingToolkit: unknowns, parameters, independent_variable, observed, - get_iv, get_observed + get_iv, get_observed # Local Difference operator (removed from Symbolics v7) include("./difference.jl") @@ -30,13 +30,13 @@ using QuadGK using Statistics using StatsBase @reexport using StatsBase: rss, r2, aic, aicc, bic, summarystats, loglikelihood, - nullloglikelihood, nobs, dof + nullloglikelihood, nobs, dof using DataInterpolations @reexport using DataInterpolations: ConstantInterpolation, LinearInterpolation, - QuadraticInterpolation, LagrangeInterpolation, - QuadraticSpline, CubicSpline, BSplineInterpolation, - BSplineApprox, Curvefit + QuadraticInterpolation, LagrangeInterpolation, + QuadraticSpline, CubicSpline, BSplineInterpolation, + BSplineApprox, Curvefit @reexport using MLUtils: splitobs, DataLoader @reexport using StatsBase: ZScoreTransform, UnitRangeTransform diff --git a/src/basis/build_function.jl b/src/basis/build_function.jl index a0f099297..186e1407e 100644 --- a/src/basis/build_function.jl +++ b/src/basis/build_function.jl @@ -3,31 +3,38 @@ struct DataDrivenFunction{IMPL, CTRLS, F1, F2} <: AbstractDataDrivenFunction{IMP f_iip::F2 end -function DataDrivenFunction(rhs, implicits, states, parameters, iv, controls, - eval_expression = false) +function DataDrivenFunction( + rhs, implicits, states, parameters, iv, controls, + eval_expression = false + ) _is_implicit = !isempty(implicits) _is_controlled = !isempty(controls) if !eval_expression f_oop, - f_iip = build_function(rhs, + f_iip = build_function( + rhs, value.(implicits), value.(states), value.(parameters), [value(iv)], value.(controls), - expression = Val{false}) + expression = Val{false} + ) else ex_oop, - ex_iip = build_function(rhs, + ex_iip = build_function( + rhs, value.(implicits), value.(states), value.(parameters), [value(iv)], value.(controls), - expression = Val{true}) + expression = Val{true} + ) f_oop = eval(ex_oop) f_iip = eval(ex_iip) end return DataDrivenFunction{_is_implicit, _is_controlled, typeof(f_oop), typeof(f_iip)}( f_oop, - f_iip) + f_iip + ) end _apply_function(f::DataDrivenFunction, du, u, p, t, c) = begin @@ -37,88 +44,121 @@ end function _apply_function!(f::DataDrivenFunction, res, du, u, p, t, c) (; f_iip) = f - f_iip(res, du, u, p, t, c) + return f_iip(res, du, u, p, t, c) end # Dispatch -# OOP +# OOP # Without controls or implicits -function (f::DataDrivenFunction{false, false})(u::AbstractVector, p::P, - t::Number) where { +function (f::DataDrivenFunction{false, false})( + u::AbstractVector, p::P, + t::Number + ) where { P <: - Union{AbstractArray, Tuple -}} - _apply_function(f, __EMPTY_VECTOR, u, p, t, __EMPTY_VECTOR) + Union{ + AbstractArray, Tuple, + }, + } + return _apply_function(f, __EMPTY_VECTOR, u, p, t, __EMPTY_VECTOR) end # Without implicits, with controls -function (f::DataDrivenFunction{false, true})(u::AbstractVector, p::P, t::Number, - c::AbstractVector) where { +function (f::DataDrivenFunction{false, true})( + u::AbstractVector, p::P, t::Number, + c::AbstractVector + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_function(f, __EMPTY_VECTOR, u, p, t, c) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_function(f, __EMPTY_VECTOR, u, p, t, c) end # With implicit, without controls -function (f::DataDrivenFunction{true, false})(du::AbstractVector, u::AbstractVector, p::P, - t::Number) where { +function (f::DataDrivenFunction{true, false})( + du::AbstractVector, u::AbstractVector, p::P, + t::Number + ) where { P <: - Union{AbstractArray, Tuple}} - _apply_function(f, du, u, p, t, __EMPTY_VECTOR) + Union{AbstractArray, Tuple}, + } + return _apply_function(f, du, u, p, t, __EMPTY_VECTOR) end # With implicit and controls -function (f::DataDrivenFunction{true, true})(du::AbstractVector, u::AbstractVector, p::P, +function (f::DataDrivenFunction{true, true})( + du::AbstractVector, u::AbstractVector, p::P, t::Number, - c::AbstractVector) where { + c::AbstractVector + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_function(f, du, u, p, t, c) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_function(f, du, u, p, t, c) end -# IIP +# IIP # Without controls or implicits function (f::DataDrivenFunction{false, false})( res::AbstractVector, u::AbstractVector, p::P, - t::Number) where { + t::Number + ) where { P <: - Union{AbstractArray, Tuple -}} - _apply_function!(f, res, __EMPTY_VECTOR, u, p, t, __EMPTY_VECTOR) + Union{ + AbstractArray, Tuple, + }, + } + return _apply_function!(f, res, __EMPTY_VECTOR, u, p, t, __EMPTY_VECTOR) end # Without implicits, with controls -function (f::DataDrivenFunction{false, true})(res::AbstractVector, u::AbstractVector, p::P, +function (f::DataDrivenFunction{false, true})( + res::AbstractVector, u::AbstractVector, p::P, t::Number, - c::AbstractVector) where { + c::AbstractVector + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_function!(f, res, __EMPTY_VECTOR, u, p, t, c) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_function!(f, res, __EMPTY_VECTOR, u, p, t, c) end # With implicit, without controls -function (f::DataDrivenFunction{true, false})(res::AbstractVector, du::AbstractVector, +function (f::DataDrivenFunction{true, false})( + res::AbstractVector, du::AbstractVector, u::AbstractVector, p::P, - t::Number) where { + t::Number + ) where { P <: - Union{AbstractArray, Tuple}} - _apply_function!(f, res, du, u, p, t, __EMPTY_VECTOR) + Union{AbstractArray, Tuple}, + } + return _apply_function!(f, res, du, u, p, t, __EMPTY_VECTOR) end # With implicit and controls -function (f::DataDrivenFunction{true, true})(res::AbstractVector, du::AbstractVector, +function (f::DataDrivenFunction{true, true})( + res::AbstractVector, du::AbstractVector, u::AbstractVector, p::P, t::Number, - c::AbstractVector) where { + c::AbstractVector + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_function!(f, res, du, u, p, t, c) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_function!(f, res, du, u, p, t, c) end ## @@ -131,97 +171,134 @@ maybeview(x, id) = x function _check_array_inputs(res, du, u, p, t, c) # Collect the keys here n_obs = size(u, 2) - @assert n_obs==length(t) "Number of observations $(n_obs) does not match timepoints $(length(t))" + @assert n_obs == length(t) "Number of observations $(n_obs) does not match timepoints $(length(t))" isempty(du) || - @assert n_obs==size(du, 2) "Number of observations $(n_obs) does not match implicits $(size(du, 2))" + @assert n_obs == size(du, 2) "Number of observations $(n_obs) does not match implicits $(size(du, 2))" isempty(c) || - @assert n_obs==size(c, 2) "Number of observations $(n_obs) does not match controls $(size(c, 2))" + @assert n_obs == size(c, 2) "Number of observations $(n_obs) does not match controls $(size(c, 2))" - isempty(res) || - @assert n_obs==size(res, 2) "Number of observations $(n_obs) does not match residuals $(size(res, 2))" + return isempty(res) || + @assert n_obs == size(res, 2) "Number of observations $(n_obs) does not match residuals $(size(res, 2))" end -function _apply_vec_function(f::DataDrivenFunction, du::AbstractMatrix, u::AbstractMatrix, - p::AbstractVector, t::AbstractVector, c::AbstractMatrix) +function _apply_vec_function( + f::DataDrivenFunction, du::AbstractMatrix, u::AbstractMatrix, + p::AbstractVector, t::AbstractVector, c::AbstractMatrix + ) _check_array_inputs(__EMPTY_MATRIX, du, u, p, t, c) - reduce(hcat, + return reduce( + hcat, map(axes(u, 2)) do i - _apply_function(f, + _apply_function( + f, maybeview(du, i), maybeview(u, i), view(p, :), maybeview(t, i), - maybeview(c, i)) - end) + maybeview(c, i) + ) + end + ) end -function _apply_vec_function!(f::DataDrivenFunction, res::AbstractMatrix, +function _apply_vec_function!( + f::DataDrivenFunction, res::AbstractMatrix, du::AbstractMatrix, u::AbstractMatrix, p::AbstractVector, - t::AbstractVector, c::AbstractMatrix) + t::AbstractVector, c::AbstractMatrix + ) _check_array_inputs(res, du, u, p, t, c) - foreach(axes(u, 2)) do i - _apply_function!(f, + return foreach(axes(u, 2)) do i + _apply_function!( + f, maybeview(res, i), maybeview(du, i), maybeview(u, i), view(p, :), - maybeview(t, i), maybeview(c, i)) + maybeview(t, i), maybeview(c, i) + ) end end -## OOP +## OOP -function (f::DataDrivenFunction{false, false})(u::AbstractMatrix, p::P, - t::AbstractVector) where { +function (f::DataDrivenFunction{false, false})( + u::AbstractMatrix, p::P, + t::AbstractVector + ) where { P <: Union{ - AbstractArray, - Tuple}} - _apply_vec_function(f, __EMPTY_MATRIX, u, p, t, __EMPTY_MATRIX) + AbstractArray, + Tuple, + }, + } + return _apply_vec_function(f, __EMPTY_MATRIX, u, p, t, __EMPTY_MATRIX) end -function (f::DataDrivenFunction{false, true})(u::AbstractMatrix, p::P, t::AbstractVector, - c::AbstractMatrix) where { +function (f::DataDrivenFunction{false, true})( + u::AbstractMatrix, p::P, t::AbstractVector, + c::AbstractMatrix + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_vec_function(f, __EMPTY_MATRIX, u, p, t, c) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_vec_function(f, __EMPTY_MATRIX, u, p, t, c) end -function (f::DataDrivenFunction{true, false})(du::AbstractMatrix, u::AbstractMatrix, p::P, - t::AbstractVector) where { +function (f::DataDrivenFunction{true, false})( + du::AbstractMatrix, u::AbstractMatrix, p::P, + t::AbstractVector + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_vec_function(f, du, u, p, t, __EMPTY_MATRIX) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_vec_function(f, du, u, p, t, __EMPTY_MATRIX) end -## IIP +## IIP function (f::DataDrivenFunction{false, false})( res::AbstractMatrix, u::AbstractMatrix, p::P, - t::AbstractVector) where { + t::AbstractVector + ) where { P <: Union{ - AbstractArray, - Tuple}} - _apply_vec_function!(f, res, __EMPTY_MATRIX, u, p, t, __EMPTY_MATRIX) + AbstractArray, + Tuple, + }, + } + return _apply_vec_function!(f, res, __EMPTY_MATRIX, u, p, t, __EMPTY_MATRIX) end -function (f::DataDrivenFunction{false, true})(res::AbstractMatrix, u::AbstractMatrix, p::P, +function (f::DataDrivenFunction{false, true})( + res::AbstractMatrix, u::AbstractMatrix, p::P, t::AbstractVector, - c::AbstractMatrix) where { + c::AbstractMatrix + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_vec_function!(f, res, __EMPTY_MATRIX, u, p, t, c) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_vec_function!(f, res, __EMPTY_MATRIX, u, p, t, c) end -function (f::DataDrivenFunction{true, false})(res::AbstractMatrix, du::AbstractMatrix, +function (f::DataDrivenFunction{true, false})( + res::AbstractMatrix, du::AbstractMatrix, u::AbstractMatrix, p::P, - t::AbstractVector) where { + t::AbstractVector + ) where { P <: - Union{AbstractArray, - Tuple}} - _apply_vec_function!(f, res, du, u, p, t, __EMPTY_MATRIX) + Union{ + AbstractArray, + Tuple, + }, + } + return _apply_vec_function!(f, res, du, u, p, t, __EMPTY_MATRIX) end diff --git a/src/basis/type.jl b/src/basis/type.jl index 342c416cd..fae520f55 100644 --- a/src/basis/type.jl +++ b/src/basis/type.jl @@ -69,8 +69,10 @@ struct Basis{IMPL, CTRLS} <: AbstractBasis """Internal systems""" systems::Vector{AbstractBasis} - function Basis(eqs, states, ctrls, ps, observed, iv, implicit, f, name, systems; - checks::Bool = true) + function Basis( + eqs, states, ctrls, ps, observed, iv, implicit, f, name, systems; + checks::Bool = true + ) if checks # Currently do nothing here #check_variables(dvs, iv) @@ -82,12 +84,14 @@ struct Basis{IMPL, CTRLS} <: AbstractBasis imp_ = !isempty(implicit) ctrls_ = !isempty(ctrls) - new{imp_, ctrls_}(eqs, states, ctrls, ps, observed, iv, implicit, f, name, systems) + return new{imp_, ctrls_}(eqs, states, ctrls, ps, observed, iv, implicit, f, name, systems) end end -function __preprocess_basis(eqs, states, ctrls, ps, observed, iv, implicit, name, systems, - simplify, linear_independent, eval_expression) +function __preprocess_basis( + eqs, states, ctrls, ps, observed, iv, implicit, name, systems, + simplify, linear_independent, eval_expression + ) # Check for iv iv === nothing && (iv = Symbolics.variable(:t)) iv = value(iv) @@ -95,58 +99,68 @@ function __preprocess_basis(eqs, states, ctrls, ps, observed, iv, implicit, name eqs = Symbolics.scalarize(eqs) lhs = isa(eqs, AbstractVector{Equation}) ? - map(Base.Fix2(getfield, :lhs), eqs) : - map(Base.Fix1(Symbolics.variable, :φ), 1:length(eqs)) + map(Base.Fix2(getfield, :lhs), eqs) : + map(Base.Fix1(Symbolics.variable, :φ), 1:length(eqs)) rhs = isa(eqs, AbstractVector{Equation}) ? - map(Base.Fix2(getfield, :rhs), eqs) : eqs + map(Base.Fix2(getfield, :rhs), eqs) : eqs rhs = Num.(rhs) lhs = Num.(lhs) # Scalarize all variables states, controls, - parameters, - implicits, - observed = value.(collect(states)), - value.(collect(ctrls)), - value.(collect(ps)), - value.(collect(implicit)), - value.(collect(observed)) + parameters, + implicits, + observed = value.(collect(states)), + value.(collect(ctrls)), + value.(collect(ps)), + value.(collect(implicit)), + value.(collect(observed)) # Filter out zeros rhs = [eq for eq in rhs if ~isequal(Num(eq), zero(Num))] rhs = linear_independent ? create_linear_independent_eqs(rhs, false) : rhs unique!(rhs, simplify) - f = DataDrivenFunction(rhs, + f = DataDrivenFunction( + rhs, implicits, states, parameters, iv, - controls, eval_expression) + controls, eval_expression + ) eqs = reduce(vcat, map(Symbolics.Equation, lhs, rhs); init = Equation[]) eqs = isa(eqs, AbstractVector) ? collect(eqs) : [collect(eqs)] return eqs, states, controls, parameters, observed, iv, implicits, f, name, - systems + systems end ## Constructors -function Basis(eqs::AbstractVector, states::AbstractVector; +function Basis( + eqs::AbstractVector, states::AbstractVector; parameters::AbstractVector = [], iv = nothing, controls::AbstractVector = [], implicits = [], observed::AbstractVector = [], name = gensym(:Basis), simplify = false, linear_independent = false, eval_expression = false, - kwargs...) - return Basis(__preprocess_basis(eqs, states, controls, parameters, observed, iv, - implicits, name, AbstractBasis[], simplify, - linear_independent, eval_expression)...) -end - -function Basis(f::Function, states::AbstractVector; parameters::AbstractVector = [], + kwargs... + ) + return Basis( + __preprocess_basis( + eqs, states, controls, parameters, observed, iv, + implicits, name, AbstractBasis[], simplify, + linear_independent, eval_expression + )... + ) +end + +function Basis( + f::Function, states::AbstractVector; parameters::AbstractVector = [], controls::AbstractVector = [], implicits::AbstractVector = [], - iv = nothing, kwargs...) + iv = nothing, kwargs... + ) isnothing(iv) && (iv = Num(Variable(:t))) try @@ -159,8 +173,10 @@ function Basis(f::Function, states::AbstractVector; parameters::AbstractVector = else eqs = f(implicits, states, parameters, iv, controls) end - return Basis(eqs, states, parameters = parameters, iv = iv, controls = controls, - implicits = implicits; kwargs...) + return Basis( + eqs, states, parameters = parameters, iv = iv, controls = controls, + implicits = implicits; kwargs... + ) catch e rethrow(e) end @@ -203,6 +219,7 @@ end continue end end + return end @inline function Base.print(io::IO, x::AbstractBasis, fullview::Bool) @@ -236,6 +253,7 @@ end for (i, eq) in enumerate(equations(x)) println(io, "$i : $(eq.lhs) = $(eq.rhs)") end + return end ## Getters @@ -287,100 +305,134 @@ get_f(b::AbstractBasis) = getfield(b, :f) # OOP # Without controls or implicits -function (b::Basis{false, false})(u::AbstractVector, p::P = parameters(b), - t::Number = get_iv(b)) where { +function (b::Basis{false, false})( + u::AbstractVector, p::P = parameters(b), + t::Number = get_iv(b) + ) where { P <: - Union{AbstractArray, Tuple}} + Union{AbstractArray, Tuple}, + } f = get_f(b) - f(u, p, t) + return f(u, p, t) end # Without implicits, with controls -function (b::Basis{false, true})(u::AbstractVector, +function (b::Basis{false, true})( + u::AbstractVector, p::P = parameters(b), t::Number = get_iv(b), - c::AbstractVector = controls(b)) where { + c::AbstractVector = controls(b) + ) where { P <: Union{ - AbstractArray, - Tuple}} + AbstractArray, + Tuple, + }, + } f = get_f(b) - f(u, p, t, c) + return f(u, p, t, c) end # With implicit, without controls -function (b::Basis{true, false})(du::AbstractVector, u::AbstractVector, +function (b::Basis{true, false})( + du::AbstractVector, u::AbstractVector, p::P = parameters(b), - t::Number = get_iv(b)) where { + t::Number = get_iv(b) + ) where { P <: - Union{AbstractArray, Tuple}} + Union{AbstractArray, Tuple}, + } f = get_f(b) - f(du, u, p, t) + return f(du, u, p, t) end # With implicit and controls -function (b::Basis{true, true})(du::AbstractVector, u::AbstractVector, +function (b::Basis{true, true})( + du::AbstractVector, u::AbstractVector, p::P = parameters(b), t::Number = get_iv(b), - c::AbstractVector = controls(b)) where { + c::AbstractVector = controls(b) + ) where { P <: - Union{AbstractArray, - Tuple}} + Union{ + AbstractArray, + Tuple, + }, + } f = get_f(b) - f(du, u, p, t, c) + return f(du, u, p, t, c) end # Array -function (b::Basis{false, false})(u::AbstractMatrix, p::P, - t::AbstractVector) where {P <: - Union{AbstractArray, Tuple}} +function (b::Basis{false, false})( + u::AbstractMatrix, p::P, + t::AbstractVector + ) where { + P <: + Union{AbstractArray, Tuple}, + } f = get_f(b) - f(u, p, t) + return f(u, p, t) end -function (b::Basis{true, false})(du::AbstractMatrix, u::AbstractMatrix, p::P, - t::AbstractVector) where {P <: Union{AbstractArray, Tuple}} +function (b::Basis{true, false})( + du::AbstractMatrix, u::AbstractMatrix, p::P, + t::AbstractVector + ) where {P <: Union{AbstractArray, Tuple}} f = get_f(b) - f(du, u, p, t) + return f(du, u, p, t) end -function (b::Basis{false, true})(u::AbstractMatrix, p::P, t::AbstractVector, - c::AbstractMatrix) where {P <: Union{AbstractArray, Tuple}} +function (b::Basis{false, true})( + u::AbstractMatrix, p::P, t::AbstractVector, + c::AbstractMatrix + ) where {P <: Union{AbstractArray, Tuple}} f = get_f(b) - f(u, p, t, c) + return f(u, p, t, c) end -function (b::Basis{true, true})(du::AbstractMatrix, u::AbstractMatrix, p::P, +function (b::Basis{true, true})( + du::AbstractMatrix, u::AbstractMatrix, p::P, t::AbstractVector, - c::AbstractMatrix) where {P <: Union{AbstractArray, Tuple}} + c::AbstractMatrix + ) where {P <: Union{AbstractArray, Tuple}} f = get_f(b) - f(du, u, p, t, c) + return f(du, u, p, t, c) end -function (b::Basis{false, false})(res::AbstractMatrix, u::AbstractMatrix, p::P, - t::AbstractVector) where {P <: - Union{AbstractArray, Tuple}} +function (b::Basis{false, false})( + res::AbstractMatrix, u::AbstractMatrix, p::P, + t::AbstractVector + ) where { + P <: + Union{AbstractArray, Tuple}, + } f = get_f(b) - f(res, u, p, t) + return f(res, u, p, t) end function (b::Basis{true, false})( res::AbstractMatrix, du::AbstractMatrix, u::AbstractMatrix, p::P, - t::AbstractVector) where {P <: Union{AbstractArray, Tuple}} + t::AbstractVector + ) where {P <: Union{AbstractArray, Tuple}} f = get_f(b) - f(res, du, u, p, t) + return f(res, du, u, p, t) end -function (b::Basis{false, true})(res::AbstractMatrix, u::AbstractMatrix, p::P, +function (b::Basis{false, true})( + res::AbstractMatrix, u::AbstractMatrix, p::P, t::AbstractVector, - c::AbstractMatrix) where {P <: Union{AbstractArray, Tuple}} + c::AbstractMatrix + ) where {P <: Union{AbstractArray, Tuple}} f = get_f(b) - f(res, u, p, t, c) + return f(res, u, p, t, c) end -function (b::Basis{true, true})(res::AbstractMatrix, du::AbstractMatrix, u::AbstractMatrix, +function (b::Basis{true, true})( + res::AbstractMatrix, du::AbstractMatrix, u::AbstractMatrix, p::P, t::AbstractVector, - c::AbstractMatrix) where {P <: Union{AbstractArray, Tuple}} + c::AbstractMatrix + ) where {P <: Union{AbstractArray, Tuple}} f = get_f(b) - f(res, du, u, p, t, c) + return f(res, du, u, p, t, c) end ## Information and Iteration @@ -396,9 +448,11 @@ Base.iterate(x::B, id) where {B <: AbstractBasis} = iterate(equations(x), id) ## Internal update function __update!(b::AbstractBasis, eval_expression = false) - ff = DataDrivenFunction([bi.rhs for bi in collect(equations(b))], + ff = DataDrivenFunction( + [bi.rhs for bi in collect(equations(b))], implicit_variables(b), unknowns(b), parameters(b), [get_iv(b)], - controls(b), eval_expression) + controls(b), eval_expression + ) @set! b.f = ff return end @@ -421,21 +475,24 @@ If control variables are defined, the function can also be called by `f(u,p,t,co If the Jacobian with respect to other variables is needed, it can be passed via a second argument. """ function jacobian(x::Basis, eval_expression::Bool = false) - jacobian( - x, unknowns(x), eval_expression) + return jacobian( + x, unknowns(x), eval_expression + ) end function jacobian(x::Basis, s, eval_expression::Bool = false) j = Symbolics.jacobian([xi.rhs for xi in equations(x)], s) - return DataDrivenFunction(j, + return DataDrivenFunction( + j, implicit_variables(x), unknowns(x), parameters(x), [get_iv(x)], - controls(x), eval_expression) + controls(x), eval_expression + ) end ## Utilities function Base.deleteat!(b::Symbolics.Arr{T, N}, idxs) where {T, N} - deleteat!(Symbolics.unwrap(b), idxs) + return deleteat!(Symbolics.unwrap(b), idxs) end ## Interfacing && merging @@ -465,7 +522,7 @@ function Base.unique!(b::Basis, simplify_eqs = false; eval_expression = false) end deleteat!(equations(b), idx) simplify_eqs && map(ModelingToolkit.simplify, equations(b)) - __update!(b, eval_expression) + return __update!(b, eval_expression) end function Base.deleteat!(b::Basis, inds; eval_expression = false) @@ -474,8 +531,10 @@ function Base.deleteat!(b::Basis, inds; eval_expression = false) return end -function Base.push!(b::Basis, eqs::AbstractArray, simplify_eqs = true; - eval_expression = false) +function Base.push!( + b::Basis, eqs::AbstractArray, simplify_eqs = true; + eval_expression = false + ) @inbounds for eq in eqs push!(b, eq, false) end @@ -517,7 +576,7 @@ function Base.isequal(x::Basis, y::Basis) length(x) == length(y) || return false yrhs = [yi.rhs for yi in equations(y)] xrhs = [xi.rhs for xi in equations(x)] - isequal(yrhs, xrhs) + return isequal(yrhs, xrhs) end """ @@ -538,7 +597,7 @@ function get_parameter_values(x::Basis) if isempty(ps) return Float64[] end - map(ps) do p + return map(ps) do p # In Symbolics v7, hasmetadata check for VariableDefaultValue may not work # Use try-catch to handle getdefaultval which throws if no default exists val = try @@ -563,7 +622,7 @@ This extends `getmetadata` in a way that all parameters have a numeric value. Values are unwrapped from symbolic wrappers to ensure compatibility with ODEProblem. """ function get_parameter_map(x::Basis) - map(parameters(x)) do p + return map(parameters(x)) do p # In Symbolics v7, hasmetadata check for VariableDefaultValue may not work # Use try-catch to handle getdefaultval which throws if no default exists val = try @@ -619,7 +678,7 @@ end # Override show to avoid ModelingToolkit's display method that needs namespacing function Base.show(io::IO, ::MIME"text/plain", b::AbstractBasis) - Base.print(io, b) + return Base.print(io, b) end # SciMLBase interface - declare whether this is an in-place function diff --git a/src/basis/utils.jl b/src/basis/utils.jl index a44b8e965..656188778 100644 --- a/src/basis/utils.jl +++ b/src/basis/utils.jl @@ -5,7 +5,7 @@ _is_constant(x::Number) = true function _is_constant(x) # In SymbolicUtils v4+, constants are wrapped in Const type # Check using isconst if available - SymbolicUtils.isconst(x) + return SymbolicUtils.isconst(x) end count_operation(x::Number, op::Function, nested::Bool = true) = 0 @@ -29,7 +29,7 @@ function count_operation(x::SymbolicUtils.BasicSymbolic, op::Function, nested::B end function count_operation(x::Num, op::Function, nested::Bool = true) - count_operation(value(x), op, nested) + return count_operation(value(x), op, nested) end function count_operation(x, ops::AbstractArray, nested::Bool = true) @@ -37,7 +37,7 @@ function count_operation(x, ops::AbstractArray, nested::Bool = true) end function count_operation(x::AbstractArray, op::Function, nested::Bool = true) - sum([count_operation(xi, op, nested) for xi in x]) + return sum([count_operation(xi, op, nested) for xi in x]) end function count_operation(x::AbstractArray, ops::AbstractArray, nested::Bool = true) @@ -46,7 +46,7 @@ function count_operation(x::AbstractArray, ops::AbstractArray, nested::Bool = tr counter += count_operation(xi, op, nested) end - counter + return counter end function split_term!(x::AbstractArray, o, ops::AbstractArray = [+]) @@ -72,7 +72,7 @@ function split_term!(x::AbstractArray, o, ops::AbstractArray = [+]) end function split_term!(x::AbstractArray, o::Num, ops::AbstractArray = [+]) - split_term!(x, value(o), ops) + return split_term!(x, value(o), ops) end remove_constant_factor(x::Num) = remove_constant_factor(value(x)) @@ -121,28 +121,28 @@ function is_dependent(x::SymbolicUtils.BasicSymbolic, y::SymbolicUtils.BasicSymb # In SymbolicUtils v4, occursin was removed. Use get_variables instead. # Check if y appears in the variables of x vars = Symbolics.get_variables(x) - y in vars + return y in vars end function is_dependent(x::Any, y::SymbolicUtils.BasicSymbolic) - false + return false end function is_dependent(x::SymbolicUtils.BasicSymbolic, y::Any) - false + return false end function is_dependent(x::Num, y::Num) - is_dependent(y.val, x.val) + return is_dependent(y.val, x.val) end function is_dependent(x::Num, y::AbstractVector{Num}) - map(yi -> is_dependent(x, yi), y) + return map(yi -> is_dependent(x, yi), y) end function is_dependent(x::AbstractVector{Num}, y::AbstractVector{Num}) inds = reduce(hcat, map(xi -> is_dependent(xi, y), x)) - inds = reshape(inds, length(y), length(x)) + return inds = reshape(inds, length(y), length(x)) end is_not_dependent(x, y) = .!is_dependent(x, y) @@ -157,7 +157,7 @@ end function is_binary(f::Function, t::Type = Number) f ∈ [+, -, *, /, ^] && return true - !is_unary(f, t) + return !is_unary(f, t) end function ariety(f::Function, t::Type = Number) diff --git a/src/commonsolve.jl b/src/commonsolve.jl index af0101050..f390cfe3d 100644 --- a/src/commonsolve.jl +++ b/src/commonsolve.jl @@ -1,11 +1,12 @@ - ## INTERNAL USE ONLY # This is a way to create a datadriven problem relatively efficient and handle all algorithms -struct InternalDataDrivenProblem{A <: AbstractDataDrivenAlgorithm, B <: AbstractBasis, TD, - T <: DataLoader, F, CI, VI, PI, SI, - O <: DataDrivenCommonOptions, - P <: AbstractDataDrivenProblem, K} +struct InternalDataDrivenProblem{ + A <: AbstractDataDrivenAlgorithm, B <: AbstractBasis, TD, + T <: DataLoader, F, CI, VI, PI, SI, + O <: DataDrivenCommonOptions, + P <: AbstractDataDrivenProblem, K, + } # The Algorithm alg::A # Data and Normalization @@ -33,8 +34,10 @@ end # This is a preprocess step, which commonly returns the implicit data. # For Koopman Algorithms this is not true -function get_fit_targets(::AbstractDataDrivenAlgorithm, prob::AbstractDataDrivenProblem, - basis::AbstractBasis) +function get_fit_targets( + ::AbstractDataDrivenAlgorithm, prob::AbstractDataDrivenProblem, + basis::AbstractBasis + ) Y = get_implicit_data(prob) X = basis(prob) return X, Y @@ -44,14 +47,17 @@ end function CommonSolve.init( prob::AbstractDataDrivenProblem, alg::AbstractDataDrivenAlgorithm; options::DataDrivenCommonOptions = DataDrivenCommonOptions(), - kwargs...) - init(prob, unit_basis(prob), alg; options = options, kwargs...) + kwargs... + ) + return init(prob, unit_basis(prob), alg; options = options, kwargs...) end -function CommonSolve.init(prob::AbstractDataDrivenProblem, basis::AbstractBasis, +function CommonSolve.init( + prob::AbstractDataDrivenProblem, basis::AbstractBasis, alg::AbstractDataDrivenAlgorithm = ZeroDataDrivenAlgorithm(); options::DataDrivenCommonOptions = DataDrivenCommonOptions(), - kwargs...) + kwargs... + ) @unpack denoise, normalize, data_processing = options # This function handles preprocessing of the variables @@ -88,9 +94,11 @@ function CommonSolve.init(prob::AbstractDataDrivenProblem, basis::AbstractBasis, test, loader = data_processing(data) - return InternalDataDrivenProblem(alg, test, loader, dt, control_idx, implicit_idx, + return InternalDataDrivenProblem( + alg, test, loader, dt, control_idx, implicit_idx, parameter_idx, state_idx, - options, basis, prob, kwargs) + options, basis, prob, kwargs + ) end function CommonSolve.solve!(::InternalDataDrivenProblem{ZeroDataDrivenAlgorithm}) diff --git a/src/difference.jl b/src/difference.jl index 8f53f14ab..c6fdbed06 100644 --- a/src/difference.jl +++ b/src/difference.jl @@ -37,11 +37,11 @@ end SymbolicUtils.promote_symtype(::Difference, ::Type{T}) where {T} = T function Base.show(io::IO, D::Difference) - print(io, "Difference(", D.t, "; dt=", D.dt, ", update=", D.update, ")") + return print(io, "Difference(", D.t, "; dt=", D.dt, ", update=", D.update, ")") end Base.nameof(::Difference) = :Difference function Base.:(==)(D1::Difference, D2::Difference) - isequal(D1.t, D2.t) && isequal(D1.dt, D2.dt) && isequal(D1.update, D2.update) + return isequal(D1.t, D2.t) && isequal(D1.dt, D2.dt) && isequal(D1.update, D2.update) end Base.hash(D::Difference, u::UInt) = hash(D.dt, hash(D.t, xor(u, 0x055640d6d952f101))) diff --git a/src/problem/set.jl b/src/problem/set.jl index 1f8124c02..fdc2fbfa4 100644 --- a/src/problem/set.jl +++ b/src/problem/set.jl @@ -29,13 +29,17 @@ end # Constructor -function DataDrivenDataset(probs::Vararg{T, N}; name = gensym(:DDSet), - kwargs...) where {T <: AbstractDataDrivenProblem, N} +function DataDrivenDataset( + probs::Vararg{T, N}; name = gensym(:DDSet), + kwargs... + ) where {T <: AbstractDataDrivenProblem, N} return DataDrivenDataset(name, probs, map(length, probs)) end -function DataDrivenDataset(solutions::Vararg{T, N}; name = gensym(:DDSet), - kwargs...) where {T <: DiffEqBase.DESolution, N} +function DataDrivenDataset( + solutions::Vararg{T, N}; name = gensym(:DDSet), + kwargs... + ) where {T <: DiffEqBase.DESolution, N} probs = map(solutions) do s DataDrivenProblem(s; kwargs...) end @@ -54,7 +58,7 @@ function DirectDataset(s::NamedTuple; name = gensym(:DDSet), kwargs...) DataDrivenProblem(si[:X]; probtype = DDProbType(1), _kwargs...) end - DataDrivenDataset(probs...; name = name) + return DataDrivenDataset(probs...; name = name) end """ @@ -68,7 +72,7 @@ function DiscreteDataset(s::NamedTuple; name = gensym(:DDSet), kwargs...) _kwargs = collect_problem_kwargs(si; kwargs...) DataDrivenProblem(si[:X]; probtype = DDProbType(2), _kwargs...) end - DataDrivenDataset(probs...; name = name) + return DataDrivenDataset(probs...; name = name) end """ @@ -79,8 +83,10 @@ $(SIGNATURES) Automatically constructs derivatives via an additional collocation method, which can be either a collocation or an interpolation from `DataInterpolations.jl` wrapped by an `InterpolationMethod` provided by the `collocation` keyword argument. """ -function ContinuousDataset(s::NamedTuple; name = gensym(:DDSet), - collocation = InterpolationMethod(), kwargs...) +function ContinuousDataset( + s::NamedTuple; name = gensym(:DDSet), + collocation = InterpolationMethod(), kwargs... + ) probs = map(keys(s)) do k si = s[k] # Check for differential states @@ -94,7 +100,7 @@ function ContinuousDataset(s::NamedTuple; name = gensym(:DDSet), throw(ArgumentError("A continuous problem $(k) needs to have either derivative or time information specified!")) end end - DataDrivenDataset(probs...; name = name) + return DataDrivenDataset(probs...; name = name) end collect_problem_kwargs(s; kwargs...) = begin @@ -131,25 +137,27 @@ function Base.print(io::IO, x::DataDrivenDataset{N, C, P}) where {N, C, P} end function is_valid(x::DataDrivenDataset) - all(map(is_valid, x.probs)) + return all(map(is_valid, x.probs)) end function get_implicit_data(x::DataDrivenDataset) - reduce(hcat, map(get_implicit_data, x.probs)) + return reduce(hcat, map(get_implicit_data, x.probs)) end # We assume common parameters (for now) function ModelingToolkit.parameters(x::DataDrivenDataset, i = :) - parameters(first(x.probs), i) + return parameters(first(x.probs), i) end -function remake_problem(d::DataDrivenDataset{<:Any, <:Any, probType}; +function remake_problem( + d::DataDrivenDataset{<:Any, <:Any, probType}; p = parameters(d), - kwargs...) where {probType} + kwargs... + ) where {probType} probs = map(d.probs) do prob remake_problem(prob, p = p) end - DataDrivenDataset(probs...) + return DataDrivenDataset(probs...) end function get_oop_args(x::DataDrivenDataset) @@ -163,5 +171,5 @@ function get_oop_args(x::DataDrivenDataset) t = vcat(t, t̂) U = hcat(U, Û) end - X, p, t, U + return X, p, t, U end diff --git a/src/problem/type.jl b/src/problem/type.jl index 2927e54e8..6fc86af80 100644 --- a/src/problem/type.jl +++ b/src/problem/type.jl @@ -39,7 +39,7 @@ $(SIGNATURES) Check if the problem is parameterized. """ function is_parametrized(x::AbstractDataDrivenProblem{N, U, C}) where {N, U, C} - hasfield(typeof(x), :p) && !isempty(x.p) + return hasfield(typeof(x), :p) && !isempty(x.p) end """ @@ -90,7 +90,7 @@ prob = DiscreteDataDrivenProblem(X, t, input_signal) ``` """ struct DataDrivenProblem{dType, cType, probType} <: - AbstractDataDrivenProblem{dType, cType, probType} + AbstractDataDrivenProblem{dType, cType, probType} # Data """State measurements""" @@ -114,8 +114,10 @@ end Base.eltype(::AbstractDataDrivenProblem{T}) where {T} = T -function DataDrivenProblem(probType, X, t, DX, Y, U, p; name = gensym(:DDProblem), - kwargs...) +function DataDrivenProblem( + probType, X, t, DX, Y, U, p; name = gensym(:DDProblem), + kwargs... + ) dType = Base.promote_eltype(X, t, DX, Y, U, p) cType = !isempty(U) name = isa(name, Symbol) ? name : Symbol(name) @@ -132,11 +134,13 @@ function DataDrivenProblem(probType, X, t, DX, Y, U, p; name = gensym(:DDProblem return DataDrivenProblem{dType, cType, probType}(_promote(X, t, DX, Y, U, p)..., name) end -function remake_problem(d::DataDrivenProblem{<:Any, <:Any, probType}; +function remake_problem( + d::DataDrivenProblem{<:Any, <:Any, probType}; X = getfield(d, :X), t = getfield(d, :t), DX = getfield(d, :DX), Y = getfield(d, :Y), U = getfield(d, :U), p = getfield(d, :p), - kwargs...) where {probType} - DataDrivenProblem(probType, X, t, DX, Y, U, p; kwargs...) + kwargs... + ) where {probType} + return DataDrivenProblem(probType, X, t, DX, Y, U, p; kwargs...) end function DataDrivenProblem(probtype, X, t, DX, Y, U::F, p; kwargs...) where {F <: Function} @@ -149,14 +153,16 @@ function DataDrivenProblem(probtype, X, t, DX, Y, U::F, p; kwargs...) where {F < return DataDrivenProblem(probtype, _promote(X, t, DX, Y, u_, p)...; kwargs...) end -function DataDrivenProblem(X::AbstractMatrix; +function DataDrivenProblem( + X::AbstractMatrix; t::AbstractVector = collect(one(eltype(X)):size(X, 2)), DX::AbstractMatrix = Array{eltype(X)}(undef, 0, 0), Y::AbstractMatrix = Array{eltype(X)}(undef, 0, 0), U::F = Array{eltype(X)}(undef, 0, 0), p::Union{AbstractVector, MTKParameters} = Array{eltype(X)}(undef, 0), probtype = nothing, - kwargs...) where {F <: Union{AbstractMatrix, Function}} + kwargs... + ) where {F <: Union{AbstractMatrix, Function}} if SS.isscimlstructure(p) _p, _, _ = SS.canonicalize(SS.Tunable(), p) else @@ -184,7 +190,7 @@ function Base.print(io::IO, x::AbstractDataDrivenProblem{N, C, P}) where {N, C, n = size(x.Y, 1) isempty(x.Y) ? nothing : println("$n observed variable(s)") n = size(x.U, 1) - isempty(x.U) ? nothing : println(io, "$n control(s)") + return isempty(x.U) ? nothing : println(io, "$n control(s)") end Base.show(io::IO, x::DataDrivenProblem{N, C, P}) where {N, C, P} = summary(io, x) @@ -196,20 +202,24 @@ A time discrete `DataDrivenProblem` useable for problems of the form `f(x[i],p,t $(SIGNATURES) """ function DiscreteDataDrivenProblem(X::AbstractMatrix; kwargs...) - DataDrivenProblem(X; probtype = DDProbType(2), kwargs...) + return DataDrivenProblem(X; probtype = DDProbType(2), kwargs...) end function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector; kwargs...) - DataDrivenProblem(X; t = t, probtype = DDProbType(2), kwargs...) + return DataDrivenProblem(X; t = t, probtype = DDProbType(2), kwargs...) end -function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, U::AbstractMatrix; - kwargs...) +function DiscreteDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, U::AbstractMatrix; + kwargs... + ) return DataDrivenProblem(X; t = t, U = U, probtype = DDProbType(2), kwargs...) end -function DiscreteDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, U::Function; - kwargs...) +function DiscreteDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, U::Function; + kwargs... + ) return DataDrivenProblem(X; t = t, U = U, probtype = DDProbType(2), kwargs...) end @@ -226,30 +236,40 @@ function ContinuousDataDrivenProblem(X::AbstractMatrix, DX::AbstractMatrix; kwar return DataDrivenProblem(X; DX = DX, probtype = DDProbType(3), kwargs...) end -function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, - DX::AbstractMatrix; kwargs...) +function ContinuousDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, + DX::AbstractMatrix; kwargs... + ) return DataDrivenProblem(X; t = t, DX = DX, probtype = DDProbType(3), kwargs...) end -function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, - DX::AbstractMatrix, U::AbstractMatrix; kwargs...) +function ContinuousDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, + DX::AbstractMatrix, U::AbstractMatrix; kwargs... + ) return DataDrivenProblem(X; t = t, DX = DX, U = U, probtype = DDProbType(3), kwargs...) end -function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, +function ContinuousDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, DX::AbstractMatrix, U::F; - kwargs...) where {F <: Function} + kwargs... + ) where {F <: Function} return DataDrivenProblem(X; t = t, DX = DX, U = U, probtype = DDProbType(3), kwargs...) end -function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, - collocation = InterpolationMethod(); kwargs...) +function ContinuousDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, + collocation = InterpolationMethod(); kwargs... + ) dx, x, t = collocate_data(X, t, collocation; kwargs...) return DataDrivenProblem(x; t = t, DX = dx, probtype = DDProbType(3), kwargs...) end -function ContinuousDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, - U::AbstractMatrix, collocation; kwargs...) +function ContinuousDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, + U::AbstractMatrix, collocation; kwargs... + ) dx, x, t = collocate_data(X, t, collocation; kwargs...) return DataDrivenProblem(x; t = t, DX = dx, U = U, probtype = DDProbType(3), kwargs...) end @@ -264,14 +284,17 @@ function DirectDataDrivenProblem(X::AbstractMatrix, Y::AbstractMatrix; kwargs... return DataDrivenProblem(X; Y = Y, probtype = DDProbType(1), kwargs...) end -function DirectDataDrivenProblem(X::AbstractMatrix, t::AbstractVector, Y::AbstractMatrix; - kwargs...) +function DirectDataDrivenProblem( + X::AbstractMatrix, t::AbstractVector, Y::AbstractMatrix; + kwargs... + ) return DataDrivenProblem(X; t = t, Y = Y, probtype = DDProbType(1), kwargs...) end function DirectDataDrivenProblem( X::AbstractMatrix, t::AbstractVector, Y::AbstractMatrix, U; - kwargs...) + kwargs... + ) return DataDrivenProblem(X; t = t, Y = Y, U = U, probtype = DDProbType(1), kwargs...) end @@ -296,67 +319,71 @@ get_name(p::AbstractDataDrivenProblem) = getfield(p, :name) function states(p::AbstractDataDrivenProblem, i = :, j = :) x = getfield(p, :X) - isempty(x) ? x : getindex(x, i, j) + return isempty(x) ? x : getindex(x, i, j) end function ModelingToolkit.parameters(p::AbstractDataDrivenProblem, i = :) x = getfield(p, :p) - isempty(x) ? x : getindex(x, i) + return isempty(x) ? x : getindex(x, i) end function ModelingToolkit.independent_variable(p::AbstractDataDrivenProblem, i = :) x = getfield(p, :t) - isempty(x) ? x : getindex(x, i) + return isempty(x) ? x : getindex(x, i) end function ModelingToolkit.get_dvs(p::ABSTRACT_CONT_PROB, i = :, j = :) x = getfield(p, :DX) - isempty(x) ? x : getindex(x, i, j) + return isempty(x) ? x : getindex(x, i, j) end function ModelingToolkit.get_dvs(p::ABSTRACT_DIRECT_PROB, i = :, j = :) x = getfield(p, :Y) - isempty(x) ? x : getindex(x, i, j) + return isempty(x) ? x : getindex(x, i, j) end function ModelingToolkit.get_dvs(p::ABSTRACT_DISCRETE_PROB, i = :, j = :) - states(p, i, j) + return states(p, i, j) end function ModelingToolkit.observed(p::AbstractDataDrivenProblem, i = :, j = :) x = getfield(p, :Y) - isempty(x) ? x : getindex(x, i, j) + return isempty(x) ? x : getindex(x, i, j) end function controls(p::AbstractDataDrivenProblem, i = :, j = :) x = getfield(p, :U) - isempty(x) ? x : getindex(x, i, j) + return isempty(x) ? x : getindex(x, i, j) end function Base.getindex(p::AbstractDataDrivenProblem, i = :, j = :) - return (states(p, i, j), + return ( + states(p, i, j), ModelingToolkit.parameters(p), ModelingToolkit.independent_variable(p, j), - ModelingToolkit.controls(p, i, j)) + ModelingToolkit.controls(p, i, j), + ) end function (b::Basis{<:Any, <:Any})(p::AbstractDataDrivenProblem{<:Any, <:Any, <:Any}) f = get_f(b) - _apply_vec_function(f, get_implicit_data(p), get_oop_args(p)...) + return _apply_vec_function(f, get_implicit_data(p), get_oop_args(p)...) end -function (b::Basis{<:Any, <:Any})(res::AbstractMatrix, - p::AbstractDataDrivenProblem{<:Any, <:Any, <:Any}) +function (b::Basis{<:Any, <:Any})( + res::AbstractMatrix, + p::AbstractDataDrivenProblem{<:Any, <:Any, <:Any} + ) f = get_f(b) - _apply_vec_function!(f, res, get_implicit_data(p), get_oop_args(p)...) + return _apply_vec_function!(f, res, get_implicit_data(p), get_oop_args(p)...) end # Check for nans, inf etc function check_domain(x) - @assert all(.~isnan.(x))&&all(.~isinf.(x)) ("One or more measurements contain `NaN` or `Inf`.") + return @assert all(.~isnan.(x))&&all(.~isinf.(x)) ("One or more measurements contain `NaN` or `Inf`.") end function check_lengths(args...) - @assert all(map(x -> size(x)[end], args) .== size(first(args))[end]) "One or more measurements are not sized equally." + return @assert all(map(x -> size(x)[end], args) .== size(first(args))[end]) "One or more measurements are not sized equally." end # Return the target variables @@ -365,14 +392,16 @@ get_implicit_data(x::ABSTRACT_DISCRETE_PROB{N, C}) where {N, C} = x.X[:, 2:end] get_implicit_data(x::ABSTRACT_CONT_PROB{N, C}) where {N, C} = x.DX function get_oop_args(x::DataDrivenProblem) - map(f -> getfield(x, f), (:X, :p, :t, :U)) + return map(f -> getfield(x, f), (:X, :p, :t, :U)) end function get_oop_args(x::DataDrivenProblem{N, C, DDProbType(2)}) where {N, C} - return (x.X[:, 1:(end - 1)], + return ( + x.X[:, 1:(end - 1)], x.p, x.t[1:(end - 1)], - x.U[:, 1:(end - 1)]) + x.U[:, 1:(end - 1)], + ) end """ @@ -448,12 +477,12 @@ end macro is_applicable(problem, basis) return quote if isa($(esc(problem)), ABSTRACT_DIRECT_PROB) - @assert length(states($(esc(basis))))==size(observed($(esc(problem))), 1) "Problem and basis need to have same observed size" + @assert length(states($(esc(basis)))) == size(observed($(esc(problem))), 1) "Problem and basis need to have same observed size" else - @assert length(states($(esc(basis))))==size(states($(esc(problem))), 1) "Problem and basis need to have same state size" + @assert length(states($(esc(basis)))) == size(states($(esc(problem))), 1) "Problem and basis need to have same state size" end - @assert length(controls($(esc(basis))))==size(controls($(esc(problem))), 1) "Problem and basis need to have same control size" - @assert length(parameters($(esc(basis))))<=length(parameters($(esc(problem)))) "Problem and basis need to have consistent parameter size" + @assert length(controls($(esc(basis)))) == size(controls($(esc(problem))), 1) "Problem and basis need to have same control size" + @assert length(parameters($(esc(basis)))) <= length(parameters($(esc(problem)))) "Problem and basis need to have consistent parameter size" end end @@ -473,8 +502,10 @@ end ## DESolution dispatch -function DataDrivenProblem(sol::T; use_interpolation = false, - kwargs...) where {T <: DiffEqBase.DESolution} +function DataDrivenProblem( + sol::T; use_interpolation = false, + kwargs... + ) where {T <: DiffEqBase.DESolution} if sol.retcode != :Success throw(AssertionError("The solution is not successful. Abort.")) return diff --git a/src/solution.jl b/src/solution.jl index 92d7ea070..cf528ac44 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -27,11 +27,13 @@ struct DataDrivenSolution{T} <: AbstractDataDrivenSolution internal_problem::InternalDataDrivenProblem end -function DataDrivenSolution(b::AbstractBasis, p::AbstractDataDrivenProblem, +function DataDrivenSolution( + b::AbstractBasis, p::AbstractDataDrivenProblem, alg::AbstractDataDrivenAlgorithm, result::Vector{<:AbstractDataDrivenResult}, internal_problem::InternalDataDrivenProblem, - retcode = DDReturnCode(2)) + retcode = DDReturnCode(2) + ) ps = get_parameter_values(b) prob = remake_problem(p, p = ps) @@ -53,14 +55,16 @@ function DataDrivenSolution(b::AbstractBasis, p::AbstractDataDrivenProblem, rss = sum(abs2, residuals) end - return DataDrivenSolution{eltype(p)}(b, + return DataDrivenSolution{eltype(p)}( + b, retcode, alg, result, prob, rss, length(parameters(b)), - internal_problem) + internal_problem + ) end (r::DataDrivenSolution)(args...) = r.basis(args...) @@ -111,7 +115,7 @@ $(SIGNATURES) Returns the log-likelihood of the `DataDrivenSolution` assuming a normal distributed error. """ function StatsBase.loglikelihood(sol::DataDrivenSolution) - begin + return begin -nobs(sol) / 2 * log.(rss(sol) / nobs(sol)) end end @@ -122,7 +126,7 @@ $(SIGNATURES) Returns the number of observations of the `DataDrivenSolution`. """ function StatsBase.nobs(sol::DataDrivenSolution) - begin + return begin prod(size(get_implicit_data(getfield(sol, :prob)))) end end diff --git a/src/utils/basis_generators.jl b/src/utils/basis_generators.jl index f8538f99f..32a06adc8 100644 --- a/src/utils/basis_generators.jl +++ b/src/utils/basis_generators.jl @@ -1,11 +1,11 @@ for f in [ - :chebyshev_basis, - :sin_basis, - :cos_basis, - :fourier_basis, - :polynomial_basis, - :monomial_basis -] + :chebyshev_basis, + :sin_basis, + :cos_basis, + :fourier_basis, + :polynomial_basis, + :monomial_basis, + ] @eval $f(x, c) = $f(scalarize(x), c) end @@ -28,7 +28,7 @@ function chebyshev_basis(x::Array, coefficients::AbstractVector) eqs = Array{Num}(undef, size(x, 1) * length(coefficients)) f(x, t) = cos.(t .* acos.(x)) _generateBasis!(eqs, f, x, coefficients) - eqs + return eqs end chebyshev_basis(x::Array, terms::Int) = chebyshev_basis(x, 1:terms) @@ -43,7 +43,7 @@ function sin_basis(x::Array, coefficients::AbstractVector) eqs = Array{Num}(undef, size(x, 1) * length(coefficients)) f(x, t) = sin.(t .* x) _generateBasis!(eqs, f, x, coefficients) - eqs + return eqs end sin_basis(x::Array, terms::Int) = sin_basis(x, 1:terms) @@ -58,7 +58,7 @@ function cos_basis(x::Array, coefficients::AbstractVector) eqs = Array{Num}(undef, size(x, 1) * length(coefficients)) f(x, t) = cos.(t .* x) _generateBasis!(eqs, f, x, coefficients) - eqs + return eqs end cos_basis(x::Array, terms::Int) = cos_basis(x, 1:terms) @@ -73,7 +73,7 @@ function fourier_basis(x::Array, coefficients::AbstractVector{Int}) eqs = Array{Num}(undef, size(x, 1) * length(coefficients)) f(x, t) = iseven(t) ? cos.(t .* x ./ 2) : sin.(t .* x ./ 2) _generateBasis!(eqs, f, x, coefficients) - eqs + return eqs end fourier_basis(x::Array, terms::Int) = fourier_basis(x, 1:terms) @@ -103,7 +103,7 @@ function polynomial_basis(x::Array, degree::Int = 1) end end end - eqs + return eqs end """ @@ -126,5 +126,5 @@ function monomial_basis(x::AbstractArray, degree::Int = 1) idx = (i - 1) * n_e + j + 1 eqs[idx] = x[i]^exponents[j] end - eqs + return eqs end diff --git a/src/utils/build_basis.jl b/src/utils/build_basis.jl index 036938ed6..b139255eb 100644 --- a/src/utils/build_basis.jl +++ b/src/utils/build_basis.jl @@ -37,16 +37,16 @@ end function _generate_variables(sym::Symbol, n::Int, offset::Int = 0) xs = [Symbolics.variable(sym, i) for i in (offset + 1):(offset + n)] - Num.(map(ModelingToolkit.tovar, xs)) + return Num.(map(ModelingToolkit.tovar, xs)) end function _generate_parameters(sym::Symbol, n::Int, offset::Int = 0) xs = [Symbolics.variable(sym, i) for i in (offset + 1):(offset + n)] - Num.(map(ModelingToolkit.toparam, xs)) + return Num.(map(ModelingToolkit.toparam, xs)) end function _set_default_val(x::Num, val::T) where {T <: Number} - Num(Symbolics.setdefaultval(Symbolics.unwrap(x), val)) + return Num(Symbolics.setdefaultval(Symbolics.unwrap(x), val)) end function __build_eqs(coeff_mat, basis, prob) @@ -77,7 +77,7 @@ function __build_eqs(coeff_mat, basis, prob) end return is_implicit(basis) ? _implicit_build_eqs(basis, eqs, p, prob) : - _explicit_build_eqs(basis, eqs, p, prob) + _explicit_build_eqs(basis, eqs, p, prob) end function _explicit_build_eqs(basis, eqs, p, prob) @@ -154,20 +154,22 @@ function __construct_basis(X, b, prob, options) ps = parameters(b) eqs, ps, - implicits = is_implicit(b) ? _implicit_build_eqs(b, eqs, ps, prob) : - _explicit_build_eqs(b, eqs, ps, prob) + implicits = is_implicit(b) ? _implicit_build_eqs(b, eqs, ps, prob) : + _explicit_build_eqs(b, eqs, ps, prob) p_new = map(eachindex(p)) do i _set_default_val(Num(ps[i]), p[i]) end end - Basis(eqs, states(b), + return Basis( + eqs, states(b), parameters = p_new, iv = get_iv(b), controls = controls(b), observed = observed(b), implicits = implicits, name = gensym(:Basis), - eval_expression = eval_expresssion) + eval_expression = eval_expresssion + ) end function unit_basis(prob::DataDrivenProblem) @@ -181,5 +183,5 @@ function unit_basis(prob::DataDrivenProblem) p = _generate_parameters(:p, n_p) u = _generate_variables(:u, n_u) - Basis([x; u], x, controls = u, independent_variable = t, parameters = p) + return Basis([x; u], x, controls = u, independent_variable = t, parameters = p) end diff --git a/src/utils/collocation.jl b/src/utils/collocation.jl index 7c25a4151..652d5a75f 100644 --- a/src/utils/collocation.jl +++ b/src/utils/collocation.jl @@ -101,7 +101,7 @@ function calckernel(::TricubeKernel, t) end function calckernel(::GaussianKernel, t) - exp(-0.5 * t^2) / (sqrt(2 * π)) + return exp(-0.5 * t^2) / (sqrt(2 * π)) end function calckernel(::CosineKernel, t) @@ -113,28 +113,28 @@ function calckernel(::CosineKernel, t) end function calckernel(::LogisticKernel, t) - 1 / (exp(t) + 2 + exp(-t)) + return 1 / (exp(t) + 2 + exp(-t)) end function calckernel(::SigmoidKernel, t) - 2 / (π * (exp(t) + exp(-t))) + return 2 / (π * (exp(t) + exp(-t))) end function calckernel(::SilvermanKernel, t) - sin(abs(t) / 2 + π / 4) * 0.5 * exp(-abs(t) / sqrt(2)) + return sin(abs(t) / 2 + π / 4) * 0.5 * exp(-abs(t) / sqrt(2)) end function construct_t1(t, tpoints) - hcat(ones(eltype(tpoints), length(tpoints)), tpoints .- t) + return hcat(ones(eltype(tpoints), length(tpoints)), tpoints .- t) end function construct_t2(t, tpoints) - hcat(ones(eltype(tpoints), length(tpoints)), tpoints .- t, (tpoints .- t) .^ 2) + return hcat(ones(eltype(tpoints), length(tpoints)), tpoints .- t, (tpoints .- t) .^ 2) end function construct_w(t, tpoints, h, kernel) W = @. calckernel((kernel,), (tpoints - t) / h) / h - Diagonal(W) + return Diagonal(W) end """ @@ -211,27 +211,33 @@ function collocate_data(data, tpoints, kernel = TriangularKernel(); crop = false crop && return estimated_derivative[:, 2:(end - 1)], estimated_derivative[:, 2:(end - 1)], tpoints[2:(end - 1)] - estimated_derivative, estimated_solution, tpoints + return estimated_derivative, estimated_solution, tpoints end # Adapted to dispatch on InterpolationMethod function collocate_data(data, tpoints, interp::InterpolationMethod; kwargs...) - collocate_data(data, tpoints, tpoints, interp; kwargs...) + return collocate_data(data, tpoints, tpoints, interp; kwargs...) end -function collocate_data(data::AbstractVector, tpoints::AbstractVector, +function collocate_data( + data::AbstractVector, tpoints::AbstractVector, tpoints_sample::AbstractVector, - interp::InterpolationMethod; kwargs...) + interp::InterpolationMethod; kwargs... + ) u, du, - tpoints = collocate_data(reshape(data, 1, :), tpoints, tpoints_sample, interp; - kwargs...) + tpoints = collocate_data( + reshape(data, 1, :), tpoints, tpoints_sample, interp; + kwargs... + ) return du[1, :], u[1, :], tpoints end # Adapted to dispatch on InterpolationMethod -function collocate_data(data::AbstractMatrix{T}, tpoints::AbstractVector{T}, +function collocate_data( + data::AbstractMatrix{T}, tpoints::AbstractVector{T}, tpoints_sample::AbstractVector{T}, interp::InterpolationMethod; - crop = false, kwargs...) where {T} + crop = false, kwargs... + ) where {T} u = zeros(T, size(data, 1), length(tpoints_sample)) du = zeros(T, size(data, 1), length(tpoints_sample)) for d1 in 1:size(data, 1) diff --git a/src/utils/common_options.jl b/src/utils/common_options.jl index 98d3b971c..cc35bf646 100644 --- a/src/utils/common_options.jl +++ b/src/utils/common_options.jl @@ -1,4 +1,3 @@ - """ $(TYPEDEF) diff --git a/src/utils/data_processing.jl b/src/utils/data_processing.jl index c5843c9a5..1e6ef3a12 100644 --- a/src/utils/data_processing.jl +++ b/src/utils/data_processing.jl @@ -34,8 +34,8 @@ function (d::DataProcessing)(data::Tuple) xtrain, xtest = splitobs(data, at = split, shuffle = false) batchsize = batchsize <= 0 ? size(first(xtrain), 2) : batchsize - xtest, - DataLoader(xtrain, batchsize = batchsize, partial = partial, shuffle = true, rng = rng) + return xtest, + DataLoader(xtrain, batchsize = batchsize, partial = partial, shuffle = true, rng = rng) end (d::DataProcessing)(X, Y) = d((X, Y)) @@ -62,19 +62,19 @@ DataNormalization() = DataNormalization{Nothing}() DataNormalization(method::Type{T}) where {T} = DataNormalization{T}() function StatsBase.fit(::DataNormalization{Nothing}, data) - StatsBase.fit(ZScoreTransform, data, dims = 2, scale = false, center = false) + return StatsBase.fit(ZScoreTransform, data, dims = 2, scale = false, center = false) end function StatsBase.fit(::DataNormalization{UnitRangeTransform}, data) tf = StatsBase.fit(UnitRangeTransform, data, dims = 2) # Adapt for constants here tf.scale .= [isinf(s) ? one(eltype(s)) : s for s in tf.scale] - tf + return tf end function StatsBase.fit(::DataNormalization{ZScoreTransform}, data) tf = StatsBase.fit(ZScoreTransform, data, dims = 2, center = false) # Adapt for constants here tf.scale .= [iszero(s) ? one(eltype(s)) : s for s in tf.scale] - tf + return tf end diff --git a/src/utils/plot_recipes.jl b/src/utils/plot_recipes.jl index 9459f6e8e..30897d077 100644 --- a/src/utils/plot_recipes.jl +++ b/src/utils/plot_recipes.jl @@ -33,14 +33,18 @@ end function gather_plot_information(x::DataDrivenDataset{N, C, D}) where {N, C, D} data = map([:X, :Y, :DX, :U]) do s - reduce(hcat, map(x.probs) do prob - getproperty(prob, s) - end) + reduce( + hcat, map(x.probs) do prob + getproperty(prob, s) + end + ) end - t = reduce(vcat, map(x.probs) do prob - getproperty(prob, :t) - end) + t = reduce( + vcat, map(x.probs) do prob + getproperty(prob, :t) + end + ) ylab = is_direct(x) ? "Sample ID" : "t" @@ -102,7 +106,7 @@ end # TODO Dispatch for SINDy to show a separate plot of the active # terms in the candidates -> Need LatexStrings for that (most likely) # Get the symbols of the states for plotting -# Add extreme values for error ? +# Add extreme values for error ? @recipe function resplot(x::AbstractDataDrivenSolution) plotins, xlab = gather_plot_information(x) diff --git a/src/utils/utils.jl b/src/utils/utils.jl index a28b6132a..10a684407 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -20,7 +20,7 @@ function optimal_svht(m::Int64, n::Int64; known_noise::Bool = false) end function marcenko_pastur_density(t, lower, upper, beta) - sqrt((upper - t) .* (t - lower)) ./ (2π * beta * t) + return sqrt((upper - t) .* (t - lower)) ./ (2π * beta * t) end function incremental_marcenko_pastur(x, beta, gamma) @@ -29,7 +29,7 @@ function incremental_marcenko_pastur(x, beta, gamma) lower = (1 - sqrt(beta))^2 @inline function marcenko_pastur(x) - begin + return begin if (upper - x) * (x - lower) > 0 return marcenko_pastur_density(x, lower, upper, beta) else @@ -54,7 +54,7 @@ function median_marcenko_pastur(beta) change = true x = ones(eltype(upper), 5) y = similar(x) - while change && (upper - lower > 1e-5) + while change && (upper - lower > 1.0e-5) x = range(lower, upper, length = 5) for (i, xi) in enumerate(x) y[i] = one(eltype(x)) - incremental_marcenko_pastur(xi, beta, 0) diff --git a/test/basis/basis.jl b/test/basis/basis.jl index 3a5547265..d706e2e90 100644 --- a/test/basis/basis.jl +++ b/test/basis/basis.jl @@ -31,10 +31,14 @@ using ModelingToolkit @test isequal(b(x0), DataDrivenDiffEq.get_f(b)(x0, p, t, u)) @test isequal(b(x0, p), DataDrivenDiffEq.get_f(b)(x0, p, t, u)) @test isequal(b(x0, p, t), DataDrivenDiffEq.get_f(b)(x0, p, t, u)) - @test isequal(b(x0, p0, t, zeros(2)), - DataDrivenDiffEq.get_f(b)(x0, p0, t, zeros(2))) - @test isequal(b(x0, p0, t0, zeros(2)), - DataDrivenDiffEq.get_f(b)(x0, p0, t0, zeros(2))) + @test isequal( + b(x0, p0, t, zeros(2)), + DataDrivenDiffEq.get_f(b)(x0, p0, t, zeros(2)) + ) + @test isequal( + b(x0, p0, t0, zeros(2)), + DataDrivenDiffEq.get_f(b)(x0, p0, t0, zeros(2)) + ) @test isequal(b(x0, p0, t0, u0), DataDrivenDiffEq.get_f(b)(x0, p0, t0, u0)) end @@ -71,10 +75,14 @@ using ModelingToolkit @test isequal(b(x0), DataDrivenDiffEq.get_f(b)(x0, p, t, u)) @test isequal(b(x0, p), DataDrivenDiffEq.get_f(b)(x0, p, t, u)) @test isequal(b(x0, p, t), DataDrivenDiffEq.get_f(b)(x0, p, t, u)) - @test isequal(b(x0, p0, t, zeros(2)), - DataDrivenDiffEq.get_f(b)(x0, p0, t, zeros(2))) - @test isequal(b(x0, p0, t0, zeros(2)), - DataDrivenDiffEq.get_f(b)(x0, p0, t0, zeros(2))) + @test isequal( + b(x0, p0, t, zeros(2)), + DataDrivenDiffEq.get_f(b)(x0, p0, t, zeros(2)) + ) + @test isequal( + b(x0, p0, t0, zeros(2)), + DataDrivenDiffEq.get_f(b)(x0, p0, t0, zeros(2)) + ) @test isequal(b(x0, p0, t0, u0), DataDrivenDiffEq.get_f(b)(x0, p0, t0, u0)) end @@ -106,8 +114,10 @@ end @test isequal(ModelingToolkit.get_iv(basis), t) @test isequal(controls(basis), []) @test !DataDrivenDiffEq.is_implicit(basis) - @test DataDrivenDiffEq.count_operation((1 + cos(u[2]) * sin(u[1]))^3, - [+, cos, ^, *]) == 4 + @test DataDrivenDiffEq.count_operation( + (1 + cos(u[2]) * sin(u[1]))^3, + [+, cos, ^, *] + ) == 4 # Check array functionalities unique!(basis) diff --git a/test/basis/generators.jl b/test/basis/generators.jl index 3ee8431e6..2c511a416 100644 --- a/test/basis/generators.jl +++ b/test/basis/generators.jl @@ -9,9 +9,13 @@ u = collect(u) @test isequal(chebyshev_basis(u, 1), cos.(1 .* acos.(u))) @test isequal(fourier_basis(u, 1), sin.(1 .* u ./ 2)) @test isequal(monomial_basis(u, 1), [1; u .^ 1]) -@test isequal(polynomial_basis(u, 2), - [1; u[1]^1; u[1]^2; u[2]^1; u[1]^1 * u[2]^1; u[2]^2; u[3]^1; u[1]^1 * u[3]^1; - u[2]^1 * u[3]^1; u[3]^2]) +@test isequal( + polynomial_basis(u, 2), + [ + 1; u[1]^1; u[1]^2; u[2]^1; u[1]^1 * u[2]^1; u[2]^2; u[3]^1; u[1]^1 * u[3]^1; + u[2]^1 * u[3]^1; u[3]^2 + ] +) @test isequal(sin_basis(u, 1:2), vcat([sin.(i .* u) for i in 1:2]...)) @test isequal(cos_basis(u, 1:5), vcat([cos.(i .* u) for i in 1:5]...)) diff --git a/test/basis/implicit_basis.jl b/test/basis/implicit_basis.jl index 96cab9728..a1b88e173 100644 --- a/test/basis/implicit_basis.jl +++ b/test/basis/implicit_basis.jl @@ -56,12 +56,18 @@ end d = Difference(get_iv(basis), dt = 1.0) ∂ = Differential(get_iv(basis)) - direct_res = DataDrivenDiffEq.__construct_basis(Ξ, basis, direct_prob, - DataDrivenDiffEq.DataDrivenCommonOptions()) - discrete_res = DataDrivenDiffEq.__construct_basis(Ξ, basis, discrete_prob, - DataDrivenDiffEq.DataDrivenCommonOptions()) - cont_res = DataDrivenDiffEq.__construct_basis(Ξ, basis, cont_prob, - DataDrivenDiffEq.DataDrivenCommonOptions()) + direct_res = DataDrivenDiffEq.__construct_basis( + Ξ, basis, direct_prob, + DataDrivenDiffEq.DataDrivenCommonOptions() + ) + discrete_res = DataDrivenDiffEq.__construct_basis( + Ξ, basis, discrete_prob, + DataDrivenDiffEq.DataDrivenCommonOptions() + ) + cont_res = DataDrivenDiffEq.__construct_basis( + Ξ, basis, cont_prob, + DataDrivenDiffEq.DataDrivenCommonOptions() + ) for r in [direct_res, discrete_res, cont_res] lhs = Num.(map(eq -> eq.rhs, equations(direct_res))) @@ -74,7 +80,9 @@ end # Note: This is purely testing for functionality! basis = Basis(du .+ u, u, implicits = du) K = Float32[0 3 0; 2 0 1; 0 0 0.5] - imp_basis = DataDrivenDiffEq.__construct_basis(K, basis, discrete_prob, - DataDrivenDiffEq.DataDrivenCommonOptions()) + imp_basis = DataDrivenDiffEq.__construct_basis( + K, basis, discrete_prob, + DataDrivenDiffEq.DataDrivenCommonOptions() + ) @test all(isequal.(equations(imp_basis), collect(du .~ -u))) end diff --git a/test/commonsolve/commonsolve.jl b/test/commonsolve/commonsolve.jl index c96ce704a..9d2c3dc13 100644 --- a/test/commonsolve/commonsolve.jl +++ b/test/commonsolve/commonsolve.jl @@ -18,7 +18,7 @@ using StatsBase normalizer = DataNormalization(ZScoreTransform) transformation = StatsBase.fit(normalizer, xs) y = StatsBase.transform(transformation, xs) - @test var(y)≈one(eltype(y)) atol=1e-1 + @test var(y) ≈ one(eltype(y)) atol = 1.0e-1 end @testset "UnitRange" begin @@ -104,9 +104,11 @@ end internal::IP end - function CommonSolve.solve!(p::DataDrivenDiffEq.InternalDataDrivenProblem{ - DummyDataDrivenAlgorithm - }) + function CommonSolve.solve!( + p::DataDrivenDiffEq.InternalDataDrivenProblem{ + DummyDataDrivenAlgorithm, + } + ) return DummyDataDrivenResult(p) end @@ -161,9 +163,12 @@ end us = Symbolics.collect(us) ys = Symbolics.collect(ys) basis = Basis( - [xs .* us[3] .* exp(-ys[2]); xs[1] * sin(us[2]); xs[3] * us[1]; - sum(ys)], xs, - controls = us, implicits = ys) + [ + xs .* us[3] .* exp(-ys[2]); xs[1] * sin(us[2]); xs[3] * us[1]; + sum(ys) + ], xs, + controls = us, implicits = ys + ) prob = DirectDataDrivenProblem(x, y, U = u) alg = DummyDataDrivenAlgorithm() diff --git a/test/problem/problem.jl b/test/problem/problem.jl index 909ebbf38..b9e3d0da6 100644 --- a/test/problem/problem.jl +++ b/test/problem/problem.jl @@ -68,8 +68,10 @@ end @testset "Problem Basis Interaction" begin @variables x y z t α β u - b1 = Basis([α * x; β * y; z * t^2 + u], [x; y; z], iv = t, parameters = [α; β], - controls = [u]) + b1 = Basis( + [α * x; β * y; z * t^2 + u], [x; y; z], iv = t, parameters = [α; β], + controls = [u] + ) b2 = Basis([α * x; β * y; z * t; α], [x; y; z], iv = t, parameters = [α; β]) sample_size = 100 X1 = randn(3, sample_size) @@ -125,8 +127,10 @@ end p2 = ContinuousDataDrivenProblem(X, t, DX = DX) p3 = ContinuousDataDrivenProblem(X, t, DX = DX) - data = (prob1 = (X = X, t = t, Y = Y), - prob2 = (X = X, t = t, Y = Y, DX = DX)) + data = ( + prob1 = (X = X, t = t, Y = Y), + prob2 = (X = X, t = t, Y = Y, DX = DX), + ) s1 = DataDrivenDataset(p1, p2) s2 = ContinuousDataset(data) @@ -144,11 +148,11 @@ end # Sizes for s in sets @test size(s) == - (first(size(p1)), is_discrete(s) ? 2 * size(X, 2) - 2 : 2 * size(X, 2)) + (first(size(p1)), is_discrete(s) ? 2 * size(X, 2) - 2 : 2 * size(X, 2)) @test DataDrivenDiffEq.is_valid(s) end - # Basis handling + # Basis handling @variables x[1:size(X, 1)] b = Basis(x, x) @test b(s1) == hcat(b(p1), b(p2)) @@ -158,8 +162,10 @@ end @test hcat(X, X, X) == b(s5) # Check if misspecified data is detected - wrong_data = (prob1 = (X = X, Y = Y), - prob2 = (X = X, t = t, Y = Y)) + wrong_data = ( + prob1 = (X = X, Y = Y), + prob2 = (X = X, t = t, Y = Y), + ) @test_throws ArgumentError ContinuousDataset(wrong_data) end @@ -169,13 +175,13 @@ end # Define autoregulation system without @mtkmodel macro # (avoids macro import issues with SafeTestsets) - @parameters α=1.0 β=1.3 γ=2.0 δ=0.5 - @variables (x(time))[1:2]=[20.0, 12.0] + @parameters α = 1.0 β = 1.3 γ = 2.0 δ = 0.5 + @variables (x(time))[1:2] = [20.0, 12.0] x = collect(x) eqs = [ D(x[1]) ~ α / (1 + x[2]) - β * x[1], - D(x[2]) ~ γ / (1 + x[1]) - δ * x[2] + D(x[2]) ~ γ / (1 + x[1]) - δ * x[2], ] @named sys = System(eqs, time) diff --git a/test/problem/samplers.jl b/test/problem/samplers.jl index 5b2d0a1ed..3e91f4139 100644 --- a/test/problem/samplers.jl +++ b/test/problem/samplers.jl @@ -1,7 +1,7 @@ using DataDrivenDiffEq using LinearAlgebra -# Random Test Problem of size 100 +# Random Test Problem of size 100 X = randn(2, 100) Y = randn(1, 100) @@ -19,7 +19,7 @@ end @testset "Mini Batching" begin for n in 2:5, repeat_ in [true, false], - shuffle_ in [true, false] + shuffle_ in [true, false] batch = Batcher(n = n, repeated = repeat_, shuffle = shuffle_) train, test = batch(prob) @test length(train) == n @@ -31,8 +31,10 @@ end @testset "DataSampler" begin for repeat_ in [true, false], shuffle_ in [true, false] - ds = DataSampler(Batcher(n = 2, repeated = repeat_, shuffle = shuffle_), - Split(ratio = 0.8)) + ds = DataSampler( + Batcher(n = 2, repeated = repeat_, shuffle = shuffle_), + Split(ratio = 0.8) + ) train, test = ds(prob) @test length(test) == 20 @test sum(length, train) == 80 diff --git a/test/runtests.jl b/test/runtests.jl index a47f6ee57..afa45bbf9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,14 +6,14 @@ const GROUP = get(ENV, "GROUP", "All") function dev_subpkg(subpkg) subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) - Pkg.develop(PackageSpec(path = subpkg_path)) + return Pkg.develop(PackageSpec(path = subpkg_path)) end function activate_subpkg_env(subpkg) subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) Pkg.activate(subpkg_path) Pkg.develop(PackageSpec(path = subpkg_path)) - Pkg.instantiate() + return Pkg.instantiate() end @time begin diff --git a/test/solution/solution.jl b/test/solution/solution.jl index 6a2c59d8a..0648d94da 100644 --- a/test/solution/solution.jl +++ b/test/solution/solution.jl @@ -9,7 +9,7 @@ Random.seed!(0) p0 = randn(rng, 2) @variables x[1:3] -@parameters p₁=p0[1] p₂=p0[2] +@parameters p₁ = p0[1] p₂ = p0[2] p = [p₁; p₂] eqs = [sin(p[1] * x[1]) + x[2]; p[2] * x[3]] @@ -25,18 +25,22 @@ struct DummyDataDrivenResult{IP} <: DataDrivenDiffEq.AbstractDataDrivenResult internal::IP end -function CommonSolve.solve!(p::DataDrivenDiffEq.InternalDataDrivenProblem{ - DummyDataDrivenAlgorithm -}) +function CommonSolve.solve!( + p::DataDrivenDiffEq.InternalDataDrivenProblem{ + DummyDataDrivenAlgorithm, + } + ) return DummyDataDrivenResult(p) end prob = DirectDataDrivenProblem(x, y, p = p0) dummy_sol = solve(prob, b, DummyDataDrivenAlgorithm()) internal_problem = dummy_sol.internal -sol = DataDrivenSolution(b, prob, DataDrivenDiffEq.ZeroDataDrivenAlgorithm(), +sol = DataDrivenSolution( + b, prob, DataDrivenDiffEq.ZeroDataDrivenAlgorithm(), DataDrivenDiffEq.AbstractDataDrivenResult[dummy_sol], - internal_problem) + internal_problem +) @test dof(sol) == 2 @test rss(sol) == 0 @@ -58,14 +62,16 @@ ŷ = y .+ 0.01 * randn(rng, size(y)) prob = DirectDataDrivenProblem(x, ŷ, p = p0) dummy_sol = solve(prob, b, DummyDataDrivenAlgorithm()) internal_problem = dummy_sol.internal -sol_2 = DataDrivenSolution(b, prob, DataDrivenDiffEq.ZeroDataDrivenAlgorithm(), +sol_2 = DataDrivenSolution( + b, prob, DataDrivenDiffEq.ZeroDataDrivenAlgorithm(), DataDrivenDiffEq.AbstractDataDrivenResult[dummy_sol], - internal_problem) + internal_problem +) @test aic(sol_2) <= -1800.0 @test bic(sol_2) <= -1700.0 @test aicc(sol_2) <= -1800.0 @test 900.0 <= loglikelihood(sol_2) <= 1000.0 -@test r2(sol_2)≈0.99998 atol=1e-3 +@test r2(sol_2) ≈ 0.99998 atol = 1.0e-3 @test nobs(sol_2) == prod(size(y)) @test_nowarn summarystats(sol_2) diff --git a/test/utils.jl b/test/utils.jl index d9ea22966..650e4e913 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -21,35 +21,37 @@ end x = 0:0.1:10.0 y = permutedims(x) z = ones(1, length(x)) - # This list does not cover all kernels since some + # This list does not cover all kernels since some # are singular for m in [ - EpanechnikovKernel(), - UniformKernel(), - TriangularKernel(), - GaussianKernel(), - LogisticKernel(), - SigmoidKernel(), - SilvermanKernel() - ] + EpanechnikovKernel(), + UniformKernel(), + TriangularKernel(), + GaussianKernel(), + LogisticKernel(), + SigmoidKernel(), + SilvermanKernel(), + ] ẑ, ŷ, x̂ = collocate_data(y, x, m) - @test ẑ≈z atol=1e-1 rtol=1e-1 - @test ŷ≈y atol=1e-1 rtol=1e-1 - @test x̂≈x atol=1e-1 rtol=1e-1 + @test ẑ ≈ z atol = 1.0e-1 rtol = 1.0e-1 + @test ŷ ≈ y atol = 1.0e-1 rtol = 1.0e-1 + @test x̂ ≈ x atol = 1.0e-1 rtol = 1.0e-1 end x = 0:0.1:10.0 y = permutedims(sin.(x)) z = permutedims(cos.(x)) - for m in InterpolationMethod.([ - LinearInterpolation, - QuadraticInterpolation, - CubicSpline - ]) + for m in InterpolationMethod.( + [ + LinearInterpolation, + QuadraticInterpolation, + CubicSpline, + ] + ) ẑ, ŷ, x̂ = collocate_data(y, x, m) - @test ẑ≈z atol=1e-1 rtol=1e-1 - @test ŷ≈y atol=1e-1 rtol=1e-1 - @test x̂≈x atol=1e-1 rtol=1e-1 + @test ẑ ≈ z atol = 1.0e-1 rtol = 1.0e-1 + @test ŷ ≈ y atol = 1.0e-1 rtol = 1.0e-1 + @test x̂ ≈ x atol = 1.0e-1 rtol = 1.0e-1 end end