Adopt strict SciMLTesting 2.4 QA

docs/make.jl

@@ -1,3 +1,7 @@
+using Pkg
+
+Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
+
 using Documenter, ModelOrderReduction
 
 cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
@@ -9,8 +13,7 @@ makedocs(
     sitename = "ModelOrderReduction.jl",
     authors = "Bowen S. Zhu",
     modules = [ModelOrderReduction],
-    clean = true, doctest = false, linkcheck = true,
-    warnonly = [:missing_docs, :example_block],
+    clean = true, doctest = true, checkdocs = :exports, linkcheck = true,
     format = Documenter.HTML(
         assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/ModelOrderReduction/stable/"

docs/src/index.md

@@ -12,7 +12,7 @@ Pkg.add("ModelOrderReduction")
 ## Contributing
 
   - Please refer to the
-    [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md)
+    [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://docs.sciml.ai/ColPrac/stable/)
     for guidance on PRs, issues, and other matters relating to contributing to SciML.
 
   - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions.

src/DataReduction/POD.jl

@@ -1,4 +1,4 @@
-using TSVD: tsvd
+import TSVD as TruncatedSVD
 using RandomizedLinAlg: rsvd
 
 function matricize(VoV::Vector{Vector{T}})::Matrix{T} where {T}
@@ -17,7 +17,7 @@ function _tsvd(data::Vector{Vector{T}}, n::Int = 1; kwargs...) where {T}
     return _tsvd(mat_data, n; kwargs...)
 end
 
-_tsvd(data, n::Int = 1; kwargs...) = tsvd(data, n; kwargs...)
+_tsvd(data, n::Int = 1; kwargs...) = TruncatedSVD.tsvd(data, n; kwargs...)
 
 function _rsvd(data::Vector{Vector{T}}, n::Int, p::Int) where {T}
     mat_data = matricize(data)
@@ -30,7 +30,37 @@ _rsvd(data, n::Int, p::Int) = rsvd(data, n, p)
     POD(snapshots; min_renergy = 1.0, min_nmodes = 1, max_nmodes = length(snapshots[1]))
     POD(snapshots, nmodes)
 
-Proper orthogonal decomposition reduction problem built from state snapshots.
+Proper orthogonal decomposition reduction problem built from state snapshots. Call
+[`reduce!`](@ref) with an SVD backend to compute its basis and spectrum.
+
+# Arguments
+- `snapshots`: a state-by-snapshot matrix or a vector of state vectors.
+- `nmodes::Int`: the fixed number of retained modes. This positional form disables
+  energy-based truncation.
+
+# Keyword Arguments
+- `min_renergy = 1.0`: minimum captured relative spectral energy when selecting modes.
+- `min_nmodes::Int = 1`: lower bound on the number of retained modes.
+- `max_nmodes::Int = length(snapshots[1])`: upper bound on the number of retained modes.
+
+# Fields
+- `snapshots`: the input snapshot data.
+- `min_renergy`, `min_nmodes`, `max_nmodes`: the truncation policy.
+- `nmodes`: the selected number of retained modes.
+- `rbasis`: the reduced basis after [`reduce!`](@ref), or `missing` before reduction.
+- `renergy`: the captured relative spectral energy.
+- `spectrum`: the singular-value spectrum after [`reduce!`](@ref), or `missing` before
+  reduction.
+
+# Examples
+```jldoctest
+julia> using ModelOrderReduction
+
+julia> pod = POD([3.0 0.0; 0.0 1.0], 1);
+
+julia> pod.nmodes
+1
+```
 """
 mutable struct POD{S, T <: AbstractFloat} <: AbstractDRProblem
     # specified
@@ -88,9 +118,23 @@ function determine_truncation(
 end
 
 """
-    reduce!(pod, alg)
+    reduce!(pod::POD, alg::SVD) -> nothing
+
+Compute the reduced basis and full singular-value spectrum for `pod` using dense SVD.
+
+# Arguments
+- `pod::POD`: reduction problem to update in place.
+- `alg::SVD`: dense singular value decomposition backend.
+
+# Examples
+```jldoctest
+julia> using ModelOrderReduction
 
-Compute the reduced basis and spectrum for `pod` using the SVD backend `alg`.
+julia> pod = POD([3.0 0.0; 0.0 1.0], 1);
+
+julia> reduce!(pod, SVD()); size(pod.rbasis)
+(2, 1)
+```
 """
 function reduce!(pod::POD{S, T}, alg::SVD)::Nothing where {S, T}
     u, s, v = _svd(pod.snapshots; alg.kwargs...)
@@ -104,6 +148,16 @@ function reduce!(pod::POD{S, T}, alg::SVD)::Nothing where {S, T}
     return nothing
 end
 
+"""
+    reduce!(pod::POD, alg::TSVD) -> nothing
+
+Compute the reduced basis and truncated singular-value spectrum for `pod` using a
+truncated SVD.
+
+# Arguments
+- `pod::POD`: reduction problem to update in place.
+- `alg::TSVD`: truncated singular value decomposition backend.
+"""
 function reduce!(pod::POD{S, T}, alg::TSVD)::Nothing where {S, T}
     u, s, v = _tsvd(pod.snapshots, pod.nmodes; alg.kwargs...)
     n_max = min(size(u, 1), size(v, 1))
@@ -113,6 +167,16 @@ function reduce!(pod::POD{S, T}, alg::TSVD)::Nothing where {S, T}
     return nothing
 end
 
+"""
+    reduce!(pod::POD, alg::RSVD) -> nothing
+
+Compute the reduced basis and approximate singular-value spectrum for `pod` using a
+randomized SVD.
+
+# Arguments
+- `pod::POD`: reduction problem to update in place.
+- `alg::RSVD`: randomized singular value decomposition backend.
+"""
 function reduce!(pod::POD{S, T}, alg::RSVD)::Nothing where {S, T}
     u, s, v = _rsvd(pod.snapshots, pod.nmodes, alg.p)
     n_max = min(size(u, 1), size(v, 1))

src/Types.jl

@@ -7,7 +7,21 @@ abstract type AbstractSVD end
 """
     SVD(; kwargs...)
 
-Dense singular value decomposition backend for projection basis construction.
+Dense singular value decomposition backend for [`reduce!`](@ref).
+
+# Keyword Arguments
+- `kwargs...`: keyword arguments forwarded to `LinearAlgebra.svd`.
+
+# Fields
+- `kwargs`: the named tuple of keyword arguments forwarded to the decomposition.
+
+# Examples
+```jldoctest
+julia> using ModelOrderReduction
+
+julia> SVD() isa SVD
+true
+```
 """
 struct SVD{K <: NamedTuple} <: AbstractSVD
     kwargs::K
@@ -20,7 +34,22 @@ end
 """
     TSVD(; kwargs...)
 
-Truncated singular value decomposition backend for projection basis construction.
+Truncated singular value decomposition backend for [`reduce!`](@ref). Use this backend
+when only the requested reduced modes should be computed.
+
+# Keyword Arguments
+- `kwargs...`: keyword arguments forwarded to `TSVD.tsvd`.
+
+# Fields
+- `kwargs`: the named tuple of keyword arguments forwarded to the decomposition.
+
+# Examples
+```jldoctest
+julia> using ModelOrderReduction
+
+julia> TSVD() isa TSVD
+true
+```
 """
 struct TSVD{K <: NamedTuple} <: AbstractSVD
     kwargs::K
@@ -31,9 +60,23 @@ struct TSVD{K <: NamedTuple} <: AbstractSVD
 end
 
 """
-    RSVD([p])
+    RSVD([p = 0])
+
+Randomized singular value decomposition backend for [`reduce!`](@ref).
+
+# Arguments
+- `p::Integer = 0`: number of oversampling vectors used by `RandomizedLinAlg.rsvd`.
+
+# Fields
+- `p::Int`: the oversampling parameter.
+
+# Examples
+```jldoctest
+julia> using ModelOrderReduction
 
-Randomized singular value decomposition backend with oversampling parameter `p`.
+julia> RSVD(2).p
+2
+```
 """
 struct RSVD <: AbstractSVD
     p::Int

src/deim.jl

@@ -26,6 +26,42 @@ function deim_interpolation_indices(basis::AbstractMatrix)::Vector{Int}
     return indices
 end
 
+function _sort_observed_equations(equations::Vector{Equation})::Vector{Equation}
+    assignments = Dict{Any, Int}()
+    for (index, equation) in enumerate(equations)
+        assignments[Symbolics.unwrap(equation.lhs)] = index
+    end
+
+    dependents = [Int[] for _ in equations]
+    degrees = zeros(Int, length(equations))
+    for (index, equation) in enumerate(equations)
+        dependencies = Set(Symbolics.unwrap.(Symbolics.get_variables(Symbolics.unwrap(equation.rhs))))
+        for variable in dependencies
+            dependency = get(assignments, variable, nothing)
+            if !isnothing(dependency) && dependency != index
+                push!(dependents[dependency], index)
+                degrees[index] += 1
+            end
+        end
+    end
+
+    available = findall(iszero, degrees)
+    ordered = Equation[]
+    sizehint!(ordered, length(equations))
+    while !isempty(available)
+        index = popfirst!(available)
+        push!(ordered, equations[index])
+        for dependent in dependents[index]
+            degrees[dependent] -= 1
+            degrees[dependent] == 0 && push!(available, dependent)
+        end
+    end
+
+    length(ordered) == length(equations) ||
+        throw(ArgumentError("observed equations contain a dependency cycle"))
+    return ordered
+end
+
 """
 $(SIGNATURES)
 
@@ -60,9 +96,21 @@ the ``\\rho_i``-th column of the identity matrix ``I_n\\in\\mathbb R^{n\\times n
 - `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``.
 
+# Keyword Arguments
+- `kwargs...`: keyword arguments forwarded to `Symbolics.substitute` when constructing
+  the nonlinear reduced model.
+
 # Return
 - `reduced_rhss`: the right-hand side of ROM.
 - `linear_projection_eqs`: the linear projection mapping ``\\mathbf y=V\\hat{\\mathbf y}``.
+
+# Examples
+```julia
+reduced_rhs, projection_equations = deim(
+    full_variables, linear_coefficients, constant_terms, nonlinear_terms,
+    reduced_variables, state_basis, nonlinear_basis,
+)
+```
 """
 function deim(
         full_vars::AbstractVector, linear_coeffs::AbstractMatrix,
@@ -119,6 +167,23 @@ 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.
+
+# Arguments
+- `sys::ModelingToolkit.ODESystem`: compiled first-order ODE system without internal
+  subsystems.
+- `snapshot::AbstractMatrix`: state-by-time snapshot matrix for `sys`.
+- `pod_dim::Integer`: number of POD state modes to retain.
+
+# Keyword Arguments
+- `deim_dim::Integer = pod_dim`: number of DEIM modes for nonlinear terms.
+- `name::Symbol = Symbol(nameof(sys), :_deim)`: name assigned to the reduced system.
+- `kwargs...`: keyword arguments forwarded to ModelingToolkit transformations and
+  generated nonlinear functions.
+
+# Examples
+```julia
+reduced_system = deim(compiled_system, snapshots, 4; deim_dim = 6)
+```
 """
 function deim(
         sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
@@ -168,13 +233,8 @@ function deim(
     @set! sys.eqs = [Symbolics.scalarize(reduced_deqs); eqs]
 
     old_observed = ModelingToolkit.get_observed(sys)
-    fullstates = [map(eq -> eq.lhs, old_observed); dvs; ModelingToolkit.get_unknowns(sys)]
     new_observed = [old_observed; linear_projection_eqs]
-    new_sorted_observed = ModelingToolkit.topsort_equations(
-        new_observed, fullstates;
-        kwargs...
-    )
-    @set! sys.observed = new_sorted_observed
+    @set! sys.observed = _sort_observed_equations(new_observed)
 
     # Numeric initial conditions for the reduced unknowns from the snapshot's first column.
     # The snapshot is assumed to start at t = tspan[1], matching the FOM initial state.

test/Project.toml

@@ -13,4 +13,4 @@ MethodOfLines = "0.11"
 ModelingToolkit = "11"
 OrdinaryDiffEq = "7"
 SafeTestsets = "0.1"
-SciMLTesting = "2.1"
+SciMLTesting = "2.4"

test/qa/Project.toml

@@ -12,5 +12,5 @@ Aqua = "0.8"
 JET = "0.9, 0.10, 0.11"
 LinearAlgebra = "1"
 SafeTestsets = "0.0.1, 0.1"
-SciMLTesting = "2.1"
+SciMLTesting = "2.4"
 Test = "1"

test/qa/qa.jl

@@ -1,28 +1,3 @@
 using SciMLTesting, ModelOrderReduction, Test
 
-run_qa(
-    ModelOrderReduction;
-    explicit_imports = true,
-    api_docs_kwargs = (; rendered = true),
-    # Whole-package JET (`report_package`/`test_package`) hits a toplevel
-    # `invalid redefinition of constant ModelOrderReduction.TSVD` error: the package
-    # exports a `TSVD` struct whose name collides with the `TSVD` dependency package
-    # under JET's toplevel virtualizer. JET coverage is provided by the targeted
-    # `report_call` checks in jet_tests.jl, which sidestep the toplevel re-eval.
-    jet = false,
-    ei_kwargs = (;
-        # `topsort_equations` is owned by ModelingToolkitBase, re-exported but not
-        # marked public by ModelingToolkit, and there is no make-public plan for it.
-        # `via_owners` flags the re-export; `are_public` flags the non-public status.
-        all_qualified_accesses_via_owners = (;
-            ignore = (
-                :topsort_equations,
-            ),
-        ),
-        all_qualified_accesses_are_public = (;
-            ignore = (
-                :topsort_equations,
-            ),
-        ),
-    )
-)
+run_qa(ModelOrderReduction)

test/utils.jl

@@ -41,3 +41,11 @@ using ModelingToolkit
         @test_throws ArgumentError ModelOrderReduction.separate_terms(exprs, vars, t)
     end
 end
+
+@testset "observed equation ordering" begin
+    @variables a b c
+    equations = [c ~ b + 1, b ~ a + 1]
+    ordered = ModelOrderReduction._sort_observed_equations(equations)
+
+    @test ordered == reverse(equations)
+end