Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion docs/pages.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
pages = [
"Home" => "index.md",
"Functions" => "functions.md",
"Tutorials" => Any["tutorials/deim_FitzHugh-Nagumo.md"],
"Functions" => "functions.md",
"Developer Interface" => "interfaces.md",
]
54 changes: 54 additions & 0 deletions docs/src/interfaces.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
# Developer Interface

The supported user workflow is deliberately small: construct a [`POD`](@ref) from
state snapshots, then call [`reduce!`](@ref) with one of the built-in backends,
[`SVD`](@ref), [`TSVD`](@ref), or [`RSVD`](@ref). The public API does not require users
to interact with the abstract types described below.

## POD reduction contract

`POD` accepts either an `AbstractMatrix{T}` whose columns are state snapshots or an
`AbstractVector{<:AbstractVector{T}}` containing state snapshots, with `T <: AbstractFloat`.
All snapshots must describe the same state dimension. The positional `nmodes` form fixes
the retained rank. The keyword form selects the rank from `min_renergy`, `min_nmodes`, and
`max_nmodes`.

The reduction contract is:

1. Construct `pod = POD(snapshots, nmodes)` or `POD(snapshots; kwargs...)`.
2. Call `reduce!(pod, SVD())`, `reduce!(pod, TSVD())`, or `reduce!(pod, RSVD())`.
3. Read `pod.rbasis`, `pod.spectrum`, `pod.nmodes`, and `pod.renergy` after the call.

`reduce!` mutates `pod` and returns `nothing`. The three concrete backends are the
supported choices; their keyword arguments are forwarded to their respective SVD
implementations as described in the API reference.

```jldoctest
julia> using ModelOrderReduction

julia> snapshots = [1.0 0.0; 0.0 1.0];

julia> pod = POD(snapshots, 1);

julia> reduce!(pod, SVD()); pod.nmodes
1
```

## Internal type hierarchy

The following abstract types are implementation details. They are documented so that
contributors can understand the source hierarchy, but they are not exported and are not
stable extension points:

```@docs
ModelOrderReduction.AbstractReductionProblem
ModelOrderReduction.AbstractMORProblem
ModelOrderReduction.AbstractDRProblem
ModelOrderReduction.AbstractSVD
```

Do not subtype these types, dispatch on them from downstream packages, or add new
`reduce!` methods for custom algorithm types. Such methods would depend on internal
dispatch and field invariants that are not covered by the public API. If a supported
extension point is needed, it should first be designed and documented as a separate
interface with generic tests.
75 changes: 63 additions & 12 deletions src/DataReduction/POD.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,30 +27,44 @@ end
_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)
POD(snapshots; min_renergy = 1.0, min_nmodes = 1,
max_nmodes = length(snapshots[1])) -> POD
POD(snapshots, nmodes::Int) -> POD

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.
- `snapshots::AbstractMatrix{T}`: a state-by-snapshot matrix, where each column is a
state snapshot.
- `snapshots::AbstractVector{<:AbstractVector{T}}`: a vector of state snapshots.
- `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.
# Keywords
- `min_renergy::T = 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.
- `snapshots::S`: the input snapshot data.
- `min_renergy::T`: the minimum relative spectral energy used for truncation.
- `min_nmodes::Int`: the lower bound on the number of retained modes.
- `max_nmodes::Int`: the upper bound on the number of retained modes.
- `nmodes::Int`: the selected number of retained modes.
- `rbasis::Union{Missing, Matrix{T}}`: the reduced basis after [`reduce!`](@ref), or
`missing` before reduction.
- `renergy::T`: the captured relative spectral energy.
- `spectrum::Union{Missing, Vector{T}}`: the singular-value spectrum after
[`reduce!`](@ref), or `missing` before reduction.

# Throws
- `AssertionError`: if the snapshots are not vector-valued or the truncation bounds are
invalid.

# Returns
- `POD`: an initialized reduction problem. The basis and spectrum are populated by
[`reduce!`](@ref).

# Examples
```jldoctest
Expand Down Expand Up @@ -126,6 +140,9 @@ Compute the reduced basis and full singular-value spectrum for `pod` using dense
- `pod::POD`: reduction problem to update in place.
- `alg::SVD`: dense singular value decomposition backend.

# Returns
- `nothing`: `pod` is updated in place.

# Examples
```jldoctest
julia> using ModelOrderReduction
Expand Down Expand Up @@ -157,6 +174,23 @@ truncated SVD.
# Arguments
- `pod::POD`: reduction problem to update in place.
- `alg::TSVD`: truncated singular value decomposition backend.

# Returns
- `nothing`: `pod` is updated in place.

# Throws
- An exception propagated by `TSVD.tsvd` if the snapshot matrix is not accepted by the
truncated SVD backend.

# Examples
```jldoctest
julia> using ModelOrderReduction

julia> pod = POD([3.0 0.0; 0.0 1.0], 1);

julia> reduce!(pod, TSVD()); pod.nmodes
1
```
"""
function reduce!(pod::POD{S, T}, alg::TSVD)::Nothing where {S, T}
u, s, v = _tsvd(pod.snapshots, pod.nmodes; alg.kwargs...)
Expand All @@ -176,6 +210,23 @@ randomized SVD.
# Arguments
- `pod::POD`: reduction problem to update in place.
- `alg::RSVD`: randomized singular value decomposition backend.

# Returns
- `nothing`: `pod` is updated in place.

# Throws
- An exception propagated by `RandomizedLinAlg.rsvd` if the snapshot matrix is not
accepted by the randomized SVD backend.

# Examples
```jldoctest
julia> using ModelOrderReduction

julia> pod = POD([3.0 0.0; 0.0 1.0], 1);

julia> reduce!(pod, RSVD()); pod.nmodes
1
```
"""
function reduce!(pod::POD{S, T}, alg::RSVD)::Nothing where {S, T}
u, s, v = _rsvd(pod.snapshots, pod.nmodes, alg.p)
Expand Down
50 changes: 44 additions & 6 deletions src/Types.jl
Original file line number Diff line number Diff line change
@@ -1,15 +1,50 @@
"""
AbstractReductionProblem

Internal common supertype for model-reduction problem types.

This type is not exported and is not a supported extension point. User code should
construct one of the documented concrete problem types instead.
"""
abstract type AbstractReductionProblem end

"""
AbstractMORProblem <: AbstractReductionProblem

Internal marker type for model-order-reduction problems.

This type is not exported and is not a supported extension point. Its hierarchy may
change without a public API guarantee.
"""
abstract type AbstractMORProblem <: AbstractReductionProblem end

"""
AbstractDRProblem <: AbstractReductionProblem

Internal marker type for data-reduction problems such as [`POD`](@ref).

This type is not exported and is not a supported extension point. Use the documented
[`POD`](@ref) constructors and [`reduce!`](@ref) methods instead.
"""
abstract type AbstractDRProblem <: AbstractReductionProblem end

"""
AbstractSVD

Internal supertype for the built-in singular-value-decomposition backends.

This type is not exported and is not a supported extension point. The supported
backends are [`SVD`](@ref), [`TSVD`](@ref), and [`RSVD`](@ref); custom subtypes and
additional `reduce!` dispatches are not part of the public API.
"""
abstract type AbstractSVD end

"""
SVD(; kwargs...)
SVD(; kwargs...) -> SVD

Dense singular value decomposition backend for [`reduce!`](@ref).

# Keyword Arguments
# Keywords
- `kwargs...`: keyword arguments forwarded to `LinearAlgebra.svd`.

# Fields
Expand All @@ -32,12 +67,12 @@ struct SVD{K <: NamedTuple} <: AbstractSVD
end

"""
TSVD(; kwargs...)
TSVD(; kwargs...) -> TSVD

Truncated singular value decomposition backend for [`reduce!`](@ref). Use this backend
when only the requested reduced modes should be computed.

# Keyword Arguments
# Keywords
- `kwargs...`: keyword arguments forwarded to `TSVD.tsvd`.

# Fields
Expand All @@ -60,16 +95,19 @@ struct TSVD{K <: NamedTuple} <: AbstractSVD
end

"""
RSVD([p = 0])
RSVD([p::Int = 0]) -> RSVD

Randomized singular value decomposition backend for [`reduce!`](@ref).

# Arguments
- `p::Integer = 0`: number of oversampling vectors used by `RandomizedLinAlg.rsvd`.
- `p::Int = 0`: number of oversampling vectors used by `RandomizedLinAlg.rsvd`.

# Fields
- `p::Int`: the oversampling parameter.

# Throws
- `MethodError`: if `p` is not an `Int`.

# Examples
```jldoctest
julia> using ModelOrderReduction
Expand Down
45 changes: 32 additions & 13 deletions src/deim.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@ $(TYPEDSIGNATURES)

Compute the DEIM interpolation indices for the given projection basis.

The orthonormal `basis` should not be a sparse matrix.
The orthonormal `basis` should not be a sparse matrix. This helper is internal and is
not part of the exported API.
"""
function deim_interpolation_indices(basis::AbstractMatrix)::Vector{Int}
dim = size(basis, 2)
Expand Down Expand Up @@ -88,22 +89,34 @@ algorithm, and ``\\mathbf e_{\\rho_i}=[0,\\ldots,0,1,0,\\ldots,0]^T\\in\\mathbb
the ``\\rho_i``-th column of the identity matrix ``I_n\\in\\mathbb R^{n\\times n}``.

# 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.
- `constant_part::AbstractVector`: the constant terms ``\\underset{n\\times 1}{\\mathbf g}`` in FOM.
- `nonlinear_part::AbstractVector`: the nonlinear functions ``\\underset{n\\times 1}{\\mathbf F}`` in FOM.
- `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``.

# Keyword 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.
- `constant_part::AbstractVector`: the constant terms
``\\underset{n\\times 1}{\\mathbf g}`` in FOM.
- `nonlinear_part::AbstractVector`: the nonlinear functions
``\\underset{n\\times 1}{\\mathbf F}`` in FOM.
- `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``.

# Keywords
- `kwargs...`: keyword arguments forwarded to `Symbolics.substitute` when constructing
the nonlinear reduced model.

# Return
- `reduced_rhss`: the right-hand side of ROM.
# Returns
- `Tuple{AbstractVector, Vector{Equation}}`: the reduced right-hand side and
linear projection equations.
- `linear_projection_eqs`: the linear projection mapping ``\\mathbf y=V\\hat{\\mathbf y}``.

# Throws
- An exception from the matrix operations if the projection matrices have incompatible
dimensions.

# Examples
```julia
reduced_rhs, projection_equations = deim(
Expand Down Expand Up @@ -174,7 +187,7 @@ nonlinear expressions.
- `snapshot::AbstractMatrix`: state-by-time snapshot matrix for `sys`.
- `pod_dim::Integer`: number of POD state modes to retain.

# Keyword Arguments
# Keywords
- `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
Expand All @@ -184,6 +197,12 @@ nonlinear expressions.
```julia
reduced_system = deim(compiled_system, snapshots, 4; deim_dim = 6)
```

# Returns
- `ModelingToolkit.ODESystem`: the reduced and completed system.

# Throws
- `ArgumentError`: if the observed equations contain a dependency cycle.
"""
function deim(
sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
Expand Down
17 changes: 17 additions & 0 deletions test/interface.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
using Test, ModelOrderReduction

@testset "Public POD interface" begin
matrix_snapshots = [3.0 0.0; 0.0 1.0]
vector_snapshots = [[3.0, 0.0], [0.0, 1.0]]

for snapshots in (matrix_snapshots, vector_snapshots)
for alg in (SVD(), TSVD(), RSVD())
pod = POD(snapshots, 1)
@test reduce!(pod, alg) === nothing
@test size(pod.rbasis) == (2, 1)
@test length(pod.spectrum) >= pod.nmodes
@test pod.nmodes == 1
@test 0.0 < pod.renergy <= 1.0
end
end
end
Loading