From 01d87f5206e982cab0a4d19e8e442463d999444a Mon Sep 17 00:00:00 2001 From: smallpondtom Date: Tue, 4 Jun 2024 15:23:18 -0400 Subject: [PATCH 1/5] Implement qdeim and odeim --- src/deim.jl | 69 ++++++++++++++++++++++++++++++++++++++++++++++------ test/deim.jl | 31 ++++++++++++++++++++--- 2 files changed, 89 insertions(+), 11 deletions(-) diff --git a/src/deim.jl b/src/deim.jl index 95a19a7..bf3b07f 100644 --- a/src/deim.jl +++ b/src/deim.jl @@ -26,6 +26,37 @@ function deim_interpolation_indices(basis::AbstractMatrix)::Vector{Int} return indices end +function qdeim_interpolation_indices(basis::AbstractMatrix)::Vector{Int} + dim = size(basis, 2) + return qr(basis', ColumnNorm()).p[1:dim] +end + +function odeim_interpolation_indices(basis::AbstractMatrix, sampling_dim::Int)::Vector{Int} + dim = size(basis, 2) + @assert sampling_dim >= dim && sampling_dim <= size(basis,1) "Invalid sampling dimension" + + # Compute the first dim points with QDEIM + p = qdeim_interpolation_indices(basis) + + # select points n, ..., m + for _ in (length(p) + 1):m + F = svd(p' * basis) + S = F.S + W = F.V + gap = S[end - 1]^2 - S[end]^2 # eigengap + proj_basis = transpose(W) * basis + r = gap .+ sum(proj_basis.^2, dims=1) + r .= r .- sqrt.((gap + sum(proj_basis.^2, dims=1)).^2 .- 4 * gap * proj_basis[end, :].^2) + indices = sortperm(r, rev=true) + e = 1 + while any(indices[e] .== p) + e += 1 + end + push!(p, indices[e]) + end + return p +end + """ $(SIGNATURES) @@ -51,6 +82,10 @@ where ``P=[\\mathbf e_{\\rho_1},\\dots,\\mathbf e_{\\rho_m}]\\in\\mathbb R^{n\\t algorithm, and ``\\mathbf e_{\\rho_i}=[0,\\ldots,0,1,0,\\ldots,0]^T\\in\\mathbb R^n`` is the ``\\rho_i``-th column of the identity matrix ``I_n\\in\\mathbb R^{n\\times n}``. +Besides the standard DEIM algorithm for interpolation, this method also supports the QDEIM +and the ODEIM algorithms. The ODEIM algorithm requires an additional parameter `odeim_dim` +to specify the number of the oversampled interpolation points. + # Arguments - `full_vars::AbstractVector`: the dependent variables ``\\underset{n\\times 1}{\\mathbf y}`` in FOM. - `linear_coeffs::AbstractMatrix`: the coefficient matrix ``\\underset{n\\times n}A`` of linear terms in FOM. @@ -59,15 +94,22 @@ the ``\\rho_i``-th column of the identity matrix ``I_n\\in\\mathbb R^{n\\times n - `reduced_vars::AbstractVector`: the dependent variables ``\\underset{k\\times 1}{\\hat{\\mathbf y}}`` in the reduced-order model. - `linear_projection_matrix::AbstractMatrix`: the projection matrix ``\\underset{n\\times k}V`` for the dependent variables ``\\mathbf y``. - `nonlinear_projection_matrix::AbstractMatrix`: the projection matrix ``\\underset{n\\times m}U`` for the nonlinear functions ``\\mathbf F``. +- `interpolation_algo::Symbol`: the interpolation algorithm, which can be `:deim`, `:qdeim`, or `:odeim`. # Return - `reduced_rhss`: the right-hand side of ROM. - `linear_projection_eqs`: the linear projection mapping ``\\mathbf y=V\\hat{\\mathbf y}``. + +# References +- [DEIM](https://epubs.siam.org/doi/abs/10.1137/110822724): Chaturantabut and Sorensen, 2012. +- [QDEIM](http://epubs.siam.org/doi/10.1137/15M1019271): Drmac and Gugercin, 2016. +- [ODEIM](https://epubs.siam.org/doi/10.1137/19M1307391): Peherstorfer, Drmac, and Gugercin, 2020. """ function deim(full_vars::AbstractVector, linear_coeffs::AbstractMatrix, constant_part::AbstractVector, nonlinear_part::AbstractVector, reduced_vars::AbstractVector, linear_projection_matrix::AbstractMatrix, - nonlinear_projection_matrix::AbstractMatrix; kwargs...) + nonlinear_projection_matrix::AbstractMatrix, + interpolation_algo::Symbol; odeim_dim::Integer, kwargs...) # rename variables for convenience y = full_vars A = linear_coeffs @@ -81,7 +123,13 @@ function deim(full_vars::AbstractVector, linear_coeffs::AbstractMatrix, linear_projection_eqs = Symbolics.scalarize(y .~ V * ŷ) linear_projection_dict = Dict(eq.lhs => eq.rhs for eq in linear_projection_eqs) - indices = deim_interpolation_indices(U) # DEIM interpolation indices + if interpolation_algo == :deim + indices = deim_interpolation_indices(U) # DEIM interpolation indices + elseif interpolation_algo == :qdeim + indices = qdeim_interpolation_indices(U) # QDEIM interpolation indices + elseif interpolation_algo == :odeim + indices = odeim_interpolation_indices(U, odeim_dim) # ODEIM interpolation indices + end # the DEIM projector (not DEIM basis) satisfies # F(original_vars) ≈ projector * F(pod_basis * reduced_vars)[indices] projector = ((@view U[indices, :])' \ (U' * V))' @@ -91,8 +139,9 @@ function deim(full_vars::AbstractVector, linear_coeffs::AbstractMatrix, Â = V' * A * V ĝ = V' * g reduced_rhss = Â * ŷ + ĝ + F̂ - reduced_rhss, linear_projection_eqs + return reduced_rhss, linear_projection_eqs end + """ $(FUNCTIONNAME)( sys::ModelingToolkit.ODESystem, @@ -100,6 +149,7 @@ end pod_dim::Integer; deim_dim::Integer = pod_dim, name::Symbol = Symbol(nameof(sys), :_deim), + interpolation_algo::Symbol = :deim, kwargs... ) -> ModelingToolkit.ODESystem @@ -116,11 +166,16 @@ The LHS of equations in `sys` are all assumed to be 1st order derivatives. Use The POD basis used for DEIM interpolation is obtained from the snapshot matrix of the nonlinear terms, which is computed by executing the runtime-generated function for -nonlinear expressions. +nonlinear expressions. + +Additional to the DEIM algorithm, this function also supports the QDEIM and ODEIM. For ODEIM, +the `odeim_dim` parameter specifies the number of oversampled interpolation points. """ function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer; - deim_dim::Integer = pod_dim, name::Symbol = Symbol(nameof(sys), :_deim), - kwargs...)::ODESystem + deim_dim::Integer = pod_dim, odeim_dim::Integer = pod_dim, + name::Symbol = Symbol(nameof(sys), :_deim), + interpolation_algo::Symbol = :deim, kwargs...)::ODESystem + @assert interpolation_algo ∈ (:deim, :qdeim, :odeim) "Invalid interpolation algorithm" sys = deepcopy(sys) @set! sys.name = name @@ -158,7 +213,7 @@ function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer; reduce!(deim_reducer, TSVD()) U = deim_reducer.rbasis # DEIM projection basis - reduced_rhss, linear_projection_eqs = deim(dvs, A, g, F, ŷ, V, U; kwargs...) + reduced_rhss, linear_projection_eqs = deim(dvs, A, g, F, ŷ, V, U, interpolation_algo; odeim_dim, kwargs...) reduced_deqs = D.(ŷ) ~ reduced_rhss @set! sys.eqs = [Symbolics.scalarize(reduced_deqs); eqs] diff --git a/test/deim.jl b/test/deim.jl index 05d61d7..5310afb 100644 --- a/test/deim.jl +++ b/test/deim.jl @@ -37,18 +37,41 @@ sol = solve(ode_prob, Tsit5(), saveat = 1.0) snapshot_simpsys = Array(sol.original_sol) pod_dim = 3 + +# test DEIM deim_sys = @test_nowarn deim(simp_sys, snapshot_simpsys, pod_dim) +# test QDEIM +qdeim_sys = @test_nowarn deim(simp_sys, snapshot_simpsys, pod_dim; interpolation_algo=:qdeim) +# test ODEIM +odeim_sys = @test_nowarn deim(simp_sys, snapshot_simpsys, pod_dim; interpolation_algo=:odeim) -# check the number of dependent variables in the new system +# DEIM: check the number of dependent variables in the new system @test length(ModelingToolkit.get_states(deim_sys)) == pod_dim - deim_prob = ODEProblem(deim_sys, nothing, tspan) - deim_sol = solve(deim_prob, Tsit5(), saveat = 1.0) nₓ = length(sol[x]) nₜ = length(sol[t]) -# test solution retrieva +# test solution retrieval @test size(deim_sol[v(x, t)]) == (nₓ, nₜ) @test size(deim_sol[w(x, t)]) == (nₓ, nₜ) + +# QDEIM: check the number of dependent variables in the new system +@test length(ModelingToolkit.get_states(qdeim_sys)) == pod_dim +deim_prob = ODEProblem(qdeim_sys, nothing, tspan) +deim_sol = solve(qdeim_prob, Tsit5(), saveat = 1.0) + +# test solution retrieval +@test size(deim_sol[v(x, t)]) == (nₓ, nₜ) +@test size(deim_sol[w(x, t)]) == (nₓ, nₜ) + + +# ODEIM: check the number of dependent variables in the new system +@test length(ModelingToolkit.get_states(qdeim_sys)) == pod_dim +deim_prob = ODEProblem(qdeim_sys, nothing, tspan) +deim_sol = solve(qdeim_prob, Tsit5(), saveat = 1.0) + +# test solution retrieval +@test size(deim_sol[v(x, t)]) == (nₓ, nₜ) +@test size(deim_sol[w(x, t)]) == (nₓ, nₜ) \ No newline at end of file From d38d13f2d4fda805a8283ff9cf1f92969eaec301 Mon Sep 17 00:00:00 2001 From: smallpondtom Date: Tue, 4 Jun 2024 15:23:30 -0400 Subject: [PATCH 2/5] add script folder to ignore --- .gitignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 5a4b8af..3fbbad5 100644 --- a/.gitignore +++ b/.gitignore @@ -23,4 +23,5 @@ docs/site/ # environment. Manifest.toml -.vscode \ No newline at end of file +.vscode +scripts/ \ No newline at end of file From 2e0af6bba1692e36b09940456dbe12b36bc055c1 Mon Sep 17 00:00:00 2001 From: smallpondtom Date: Tue, 4 Jun 2024 15:24:06 -0400 Subject: [PATCH 3/5] scripts to script in ignore --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 3fbbad5..679106e 100644 --- a/.gitignore +++ b/.gitignore @@ -24,4 +24,4 @@ docs/site/ Manifest.toml .vscode -scripts/ \ No newline at end of file +script/ \ No newline at end of file From bf49b66b5d33fab391c52a81d75b018e73985ff1 Mon Sep 17 00:00:00 2001 From: smallpondtom Date: Tue, 4 Jun 2024 17:38:15 -0400 Subject: [PATCH 4/5] implemented and tested qdeim and odeim --- src/deim.jl | 32 ++++++++++++++++++++------------ test/deim.jl | 9 ++++----- test/runtests.jl | 6 +++--- 3 files changed, 27 insertions(+), 20 deletions(-) diff --git a/src/deim.jl b/src/deim.jl index bf3b07f..7ecf101 100644 --- a/src/deim.jl +++ b/src/deim.jl @@ -26,28 +26,36 @@ function deim_interpolation_indices(basis::AbstractMatrix)::Vector{Int} return indices end +""" +$(TYPEDSIGNATURES) + +Compute the QDEIM interpolation indices for the given projection basis. +""" function qdeim_interpolation_indices(basis::AbstractMatrix)::Vector{Int} dim = size(basis, 2) return qr(basis', ColumnNorm()).p[1:dim] end -function odeim_interpolation_indices(basis::AbstractMatrix, sampling_dim::Int)::Vector{Int} +""" +$(TYPEDSIGNATURES) + +Compute the ODEIM interpolation indices for the given projection basis. +""" +function odeim_interpolation_indices(basis::AbstractMatrix, m::Int)::Vector{Int} dim = size(basis, 2) - @assert sampling_dim >= dim && sampling_dim <= size(basis,1) "Invalid sampling dimension" + @assert m >= dim && m <= size(basis,1) "Invalid sampling dimension" # Compute the first dim points with QDEIM p = qdeim_interpolation_indices(basis) - # select points n, ..., m + # select points n+1, ..., m for _ in (length(p) + 1):m - F = svd(p' * basis) - S = F.S - W = F.V + _, S, W = svd(basis[p, :]) gap = S[end - 1]^2 - S[end]^2 # eigengap - proj_basis = transpose(W) * basis + proj_basis = W' * basis' r = gap .+ sum(proj_basis.^2, dims=1) - r .= r .- sqrt.((gap + sum(proj_basis.^2, dims=1)).^2 .- 4 * gap * proj_basis[end, :].^2) - indices = sortperm(r, rev=true) + r -= sqrt.((gap .+ sum(proj_basis.^2, dims=1)).^2 - 4 * gap * (proj_basis[end, :].^2)') + indices = sortperm(vec(r), rev=true) e = 1 while any(indices[e] .== p) e += 1 @@ -109,7 +117,7 @@ function deim(full_vars::AbstractVector, linear_coeffs::AbstractMatrix, constant_part::AbstractVector, nonlinear_part::AbstractVector, reduced_vars::AbstractVector, linear_projection_matrix::AbstractMatrix, nonlinear_projection_matrix::AbstractMatrix, - interpolation_algo::Symbol; odeim_dim::Integer, kwargs...) + interpolation_algo::Symbol, odeim_dim::Integer; kwargs...) # rename variables for convenience y = full_vars A = linear_coeffs @@ -172,7 +180,7 @@ Additional to the DEIM algorithm, this function also supports the QDEIM and ODEI the `odeim_dim` parameter specifies the number of oversampled interpolation points. """ function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer; - deim_dim::Integer = pod_dim, odeim_dim::Integer = pod_dim, + deim_dim::Integer = pod_dim, odeim_dim::Integer = 2*pod_dim, name::Symbol = Symbol(nameof(sys), :_deim), interpolation_algo::Symbol = :deim, kwargs...)::ODESystem @assert interpolation_algo ∈ (:deim, :qdeim, :odeim) "Invalid interpolation algorithm" @@ -213,7 +221,7 @@ function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer; reduce!(deim_reducer, TSVD()) U = deim_reducer.rbasis # DEIM projection basis - reduced_rhss, linear_projection_eqs = deim(dvs, A, g, F, ŷ, V, U, interpolation_algo; odeim_dim, kwargs...) + reduced_rhss, linear_projection_eqs = deim(dvs, A, g, F, ŷ, V, U, interpolation_algo, odeim_dim; kwargs...) reduced_deqs = D.(ŷ) ~ reduced_rhss @set! sys.eqs = [Symbolics.scalarize(reduced_deqs); eqs] diff --git a/test/deim.jl b/test/deim.jl index 5310afb..afb1f3d 100644 --- a/test/deim.jl +++ b/test/deim.jl @@ -60,17 +60,16 @@ nₜ = length(sol[t]) # QDEIM: check the number of dependent variables in the new system @test length(ModelingToolkit.get_states(qdeim_sys)) == pod_dim deim_prob = ODEProblem(qdeim_sys, nothing, tspan) -deim_sol = solve(qdeim_prob, Tsit5(), saveat = 1.0) +deim_sol = solve(deim_prob, Tsit5(), saveat = 1.0) # test solution retrieval @test size(deim_sol[v(x, t)]) == (nₓ, nₜ) @test size(deim_sol[w(x, t)]) == (nₓ, nₜ) - # ODEIM: check the number of dependent variables in the new system -@test length(ModelingToolkit.get_states(qdeim_sys)) == pod_dim -deim_prob = ODEProblem(qdeim_sys, nothing, tspan) -deim_sol = solve(qdeim_prob, Tsit5(), saveat = 1.0) +@test length(ModelingToolkit.get_states(odeim_sys)) == pod_dim +deim_prob = ODEProblem(odeim_sys, nothing, tspan) +deim_sol = solve(deim_prob, Tsit5(), saveat = 1.0) # test solution retrieval @test size(deim_sol[v(x, t)]) == (nₓ, nₜ) diff --git a/test/runtests.jl b/test/runtests.jl index 821f7d6..f8849e6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ using SafeTestsets -@safetestset "Quality Assurance" begin include("qa.jl") end -@safetestset "POD" begin include("DataReduction.jl") end -@safetestset "utils" begin include("utils.jl") end +# @safetestset "Quality Assurance" begin include("qa.jl") end +# @safetestset "POD" begin include("DataReduction.jl") end +# @safetestset "utils" begin include("utils.jl") end @safetestset "DEIM" begin include("deim.jl") end From 0e2df65ec00d61ca0c879efbfea5c077cf2ad4c4 Mon Sep 17 00:00:00 2001 From: smallpondtom Date: Tue, 4 Jun 2024 17:52:14 -0400 Subject: [PATCH 5/5] add back all the tests --- test/runtests.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index f8849e6..821f7d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ using SafeTestsets -# @safetestset "Quality Assurance" begin include("qa.jl") end -# @safetestset "POD" begin include("DataReduction.jl") end -# @safetestset "utils" begin include("utils.jl") end +@safetestset "Quality Assurance" begin include("qa.jl") end +@safetestset "POD" begin include("DataReduction.jl") end +@safetestset "utils" begin include("utils.jl") end @safetestset "DEIM" begin include("deim.jl") end