Restructure POD interface

src/DataReduction/POD.jl

@@ -1,61 +1,89 @@
-function matricize(VoV::Vector{Vector{FT}}) where {FT}
-    Matrix(reduce(hcat, VoV))
+function matricize(VoV::Vector{Vector{T}}) where {T}
+    reduce(hcat, VoV)
 end
 
-mutable struct POD{FT} <: AbstractDRProblem
-    snapshots::Union{Vector{Vector{FT}}, Matrix{FT}}
-    nmodes::Int
-    rbasis::Matrix{FT}
-    renergy::FT
+function _svd(data::Vector{Vector{T}}; kwargs...) where {T}
+    mat_data = matricize(data)
+    return _svd(mat_data; kwargs...)
+end
 
-    function POD(snaps::Vector{Vector{FT}}, nmodes::Int) where {FT}
-        errorhandle(matricize(snaps), nmodes)
-        new{eltype(snaps[1])}(snaps, nmodes, Array{FT, 2}(undef, size(snaps, 1), nmodes),
-                              FT(0))
-    end
+_svd(data; kwargs...) = svd(data; kwargs...)
 
-    function POD(snaps::Matrix{FT}, nmodes::Int) where {FT}
-        errorhandle(snaps, nmodes)
-        new{eltype(snaps)}(snaps, nmodes, Array{FT, 2}(undef, size(snaps, 1), nmodes),
-                           FT(0))
-    end
+function _tsvd(data::Vector{Vector{T}}, n::Int = 1; kwargs...) where {T}
+    mat_data = matricize(data)
+    return _tsvd(mat_data, n; kwargs...)
+end
+
+_tsvd(data, n::Int = 1; kwargs...) = tsvd(data, n; kwargs...)
+
+function _rsvd(data::Vector{Vector{T}}, n::Int, p::Int) where {T}
+    mat_data = matricize(data)
+    return _rsvd(mat_data, n, p)
+end
 
-    function POD(snaps::Adjoint{FT, Matrix{FT}}, nmodes::Int) where {FT}
-        POD(Matrix(snaps), nmodes, Array{FT, 2}(undef, size(snaps, 1), nmodes), FT(0))
+_rsvd(data, n::Int, p::Int) = rsvd(data, n, p)
+
+mutable struct POD <: AbstractDRProblem
+    # specified 
+    snapshots::Any
+    min_renergy::Any
+    min_nmodes::Int
+    max_nmodes::Int
+    # computed
+    nmodes::Int
+    rbasis::Any
+    renergy::Any
+    spectrum::Any
+    # constructers
+    function POD(snaps;
+                 min_renergy = 1.0,
+                 min_nmodes::Int = 1,
+                 max_nmodes::Int = length(snaps[1]))
+        nmodes = min_nmodes
+        errorhandle(snaps, nmodes, min_renergy, min_nmodes, max_nmodes)
+        new(snaps, min_renergy, min_nmodes, max_nmodes, nmodes, missing, 1.0, missing)
+    end
+    function POD(snaps, nmodes::Int)
+        errorhandle(snaps, nmodes, 0.0, nmodes, nmodes)
+        new(snaps, 0.0, nmodes, nmodes, nmodes, missing, 1.0, missing)
     end
 end
 
-function reduce!(pod::POD{FT}, ::SVD) where {FT}
-    op_matrix = pod.snapshots
-    if typeof(pod.snapshots) == Vector{Vector{FT}}
-        op_matrix = matricize(pod.snapshots)
+function determine_truncation(s, min_nmodes, min_renergy, max_nmodes)
+    nmodes = min_nmodes
+    overall_energy = sum(s)
+    energy = sum(s[1:nmodes]) / overall_energy
+    while energy < min_renergy && nmodes < max_nmodes
+        nmodes += 1
+        energy += s[nmodes + 1] / overall_energy
     end
-    u, s, v = svd(op_matrix)
-    pod.rbasis .= u[:, 1:(pod.nmodes)]
-    sr = s[1:(pod.nmodes)]
-    pod.renergy = sum(s) / (sum(s) + (size(op_matrix, 1) - pod.nmodes) * s[end])
+    return nmodes, energy
+end
+
+function reduce!(pod::POD, alg::SVD)
+    u, s, v = _svd(pod.snapshots; alg.kwargs...)
+    pod.nmodes, pod.renergy = determine_truncation(s, pod.min_nmodes, pod.max_nmodes,
+                                                   pod.min_renergy)
+    pod.rbasis = u[:, 1:(pod.nmodes)]
+    pod.spectrum = s
     nothing
 end
 
-function reduce!(pod::POD{FT}, ::TSVD) where {FT}
-    op_matrix = pod.snapshots
-    if typeof(pod.snapshots) == Vector{Vector{FT}}
-        op_matrix = matricize(pod.snapshots)
-    end
-    u, s, v = tsvd(op_matrix, pod.nmodes)
-    pod.renergy = sum(s) / (sum(s) + (size(op_matrix, 1) - pod.nmodes) * s[end])
-    pod.rbasis .= u
+function reduce!(pod::POD, alg::TSVD)
+    u, s, v = _tsvd(pod.snapshots, pod.nmodes; alg.kwargs...)
+    n_max = min(size(u, 1), size(v, 1))
+    pod.renergy = sum(s) / (sum(s) + (n_max - pod.nmodes) * s[end])
+    pod.rbasis = u
+    pod.spectrum = s
     nothing
 end
 
-function reduce!(pod::POD{FT}, ::RSVD) where {FT}
-    op_matrix = pod.snapshots
-    if typeof(pod.snapshots) == Vector{Vector{FT}}
-        op_matrix = matricize(pod.snapshots)
-    end
-    u, s, v = rsvd(op_matrix, pod.nmodes)
-    pod.renergy = sum(s) / (sum(s) + (size(op_matrix, 1) - pod.nmodes) * s[end])
-    pod.rbasis .= u
+function reduce!(pod::POD, alg::RSVD)
+    u, s, v = _rsvd(pod.snapshots, pod.nmodes, alg.p)
+    n_max = min(size(u, 1), size(v, 1))
+    pod.renergy = sum(s) / (sum(s) + (n_max - pod.nmodes) * s[end])
+    pod.rbasis = u
+    pod.spectrum = s
     nothing
 end
 

src/ErrorHandle.jl

@@ -1,5 +1,16 @@
-function errorhandle(data::Matrix{FT}, modes::IT) where {FT, IT}
+function errorhandle(data::AbstractMatrix, nmodes::Int, max_energy, min_nmodes, max_nmodes)
     @assert size(data, 1)>1 "State vector is expected to be vector valued."
-    s = size(data, 2)
-    @assert (modes > 0)&(modes < s) "Number of modes should be in {1,2,...,$(s-1)}."
+    s = minimum(size(data))
+    @assert 0<nmodes<=s "Number of modes should be in {1,2,...,$s}."
+    @assert min_nmodes<=max_nmodes "Minimum number of modes must lie below maximum number of modes"
+    @assert 0.0<=max_energy<=1.0 "Maxmimum relative energy must be in [0,1]"
+end
+
+function errorhandle(data::AbstractVector{T}, nmodes::Int, max_energy, min_nmodes,
+                     max_nmodes) where {T <: AbstractVector}
+    @assert size(data[1], 1)>1 "State vector is expected to be vector valued."
+    s = min(size(data, 1), size(data[1], 1))
+    @assert 0<nmodes<=s "Number of modes should be in {1,2,...,$s}."
+    @assert min_nmodes<=max_nmodes "Minimum number of modes must lie below maximum number of modes"
+    @assert 0.0<=max_energy<=1.0 "Maxmimum relative energy must be in [0,1]"
 end

src/Types.jl

@@ -3,6 +3,24 @@ abstract type AbstractMORProblem <: AbstractReductionProblem end
 abstract type AbstractDRProblem <: AbstractReductionProblem end
 
 abstract type AbstractSVD end
-struct SVD <: AbstractSVD end
-struct TSVD <: AbstractSVD end
-struct RSVD <: AbstractSVD end
+
+struct SVD <: AbstractSVD
+    kwargs::Any
+    function SVD(; kwargs...)
+        new(kwargs)
+    end
+end
+
+struct TSVD <: AbstractSVD
+    kwargs::Any
+    function TSVD(; kwargs...)
+        new(kwargs)
+    end
+end
+
+struct RSVD <: AbstractSVD
+    p::Int
+    function RSVD(p::Int = 0)
+        new(p)
+    end
+end

test/DataReduction.jl

@@ -29,31 +29,38 @@ end
 
     order = 2
     solver = SVD()
-    reducer = POD(solution, order)
-    reduce!(reducer, solver)
+    matrix_reducer = POD(solution, order)
+    snapshot_reducer = POD(sol.u, order)
+    reduce!(matrix_reducer, solver)
+    reduce!(snapshot_reducer, solver)
 
-    reducer.renergy
+    @test all(matrix_reducer.rbasis .≈ snapshot_reducer.rbasis)
+    @test matrix_reducer.renergy ≈ snapshot_reducer.renergy
 
-    # Ad-hoc tests. To be checked with Chris.
-    @test size(reducer.rbasis, 2) == reducer.nmodes
-    @test size(reducer.rbasis, 1) == size(solution, 1)
-    @test reducer.renergy > 0.9
+    @test size(matrix_reducer.rbasis, 2) == matrix_reducer.nmodes
+    @test size(matrix_reducer.rbasis, 1) == size(solution, 1)
+    @test matrix_reducer.renergy > 0.9
+
+    reducer = POD(solution, min_nmodes = 1, max_nmodes = 2, min_renergy = 0.1)
+    reduce!(reducer, solver)
+    @test reducer.renergy > 0.1
+    @test reducer.nmodes == 1
 
     order = 2
     solver = TSVD()
     reducer = POD(solution, order)
     reduce!(reducer, solver)
 
-    # Ad-hoc tests. To be checked with Chris.
     @test size(reducer.rbasis, 2) == reducer.nmodes
     @test size(reducer.rbasis, 1) == size(solution, 1)
+    @test reducer.renergy > 0.7
 
     order = 2
     solver = RSVD()
     reducer = POD(solution, order)
     reduce!(reducer, solver)
 
-    # Ad-hoc tests. To be checked with Chris.
     @test size(reducer.rbasis, 2) == reducer.nmodes
     @test size(reducer.rbasis, 1) == size(solution, 1)
+    @test reducer.renergy > 0.7
 end