Changes from prev commit.

nmodes is now forced to be an integer.
Silly transpose error fix
Update relative energy for truncated SVD algos

TODO: Modify type of POD, similar to the type of snapshots

Author: Rahul Manavalan

src/DataReduction/POD.jl

@@ -1,65 +1,65 @@
 function matricize(VoV::Vector{Vector{FT}}) where {FT}
-    Matrix(reduce(hcat,VoV)')
+    Matrix(reduce(hcat,VoV))
 end
 
-mutable struct POD{FT,IT} <: AbstractDRProblem
+mutable struct POD{FT} <: AbstractDRProblem
     snapshots::Union{Vector{Vector{FT}},Matrix{FT}}
-    nmodes::IT
+    nmodes::Int
     rbasis::Matrix{FT}
-    energy::FT
+    renergy::FT
 
-    function POD(snaps::Vector{Vector{FT}},nmodes::IT) where {FT,IT}
+    function POD(snaps::Vector{Vector{FT}},nmodes::Int) where {FT}
         errorhandle(matricize(snaps),nmodes)
-        new{eltype(snaps[1]),typeof(nmodes)}(snaps,nmodes,Array{FT,2}(undef,size(snaps,1),nmodes),FT(0))
+        new{eltype(snaps[1])}(snaps,nmodes,Array{FT,2}(undef,size(snaps,1),nmodes),FT(0))
     end
 
-    function POD(snaps::Matrix{FT},nmodes::IT) where {FT,IT}
+    function POD(snaps::Matrix{FT},nmodes::Int) where {FT}
         errorhandle(snaps,nmodes)
-        new{eltype(snaps),typeof(nmodes)}(snaps,nmodes,Array{FT,2}(undef,size(snaps,1),nmodes),FT(0))
+        new{eltype(snaps)}(snaps,nmodes,Array{FT,2}(undef,size(snaps,1),nmodes),FT(0))
     end
 
-    function POD(snaps::Adjoint{FT,Matrix{FT}},nmodes::IT) where {FT,IT}
+    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))
     end
 end
 
-function reduce!(pod::POD{FT,IT},::SVD) where {FT,IT}
+function reduce!(pod::POD{FT},::SVD) where {FT}
     op_matrix = pod.snapshots
     if typeof(pod.snapshots) == Vector{Vector{FT}}
         op_matrix = matricize(pod.snapshots)
     end
     u,s,v = svd(op_matrix)
     pod.rbasis .= u[:,1:pod.nmodes]
     sr = s[1:pod.nmodes]
-    pod.energy = sum(sr)/sum(s)
+    pod.renergy = sum(s)/(sum(s) + (size(op_matrix,1)-pod.nmodes)*s[end])
     nothing
 end
 
-function reduce!(pod::POD{FT,IT},::TSVD) where {FT,IT}
+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.energy = NaN
+    pod.renergy = sum(s)/(sum(s) + (size(op_matrix,1)-pod.nmodes)*s[end])
     pod.rbasis .= u
     nothing
 end
 
-function reduce!(pod::POD{FT,IT},::RSVD) where {FT,IT}
+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.energy = NaN
+    pod.renergy = sum(s)/(sum(s) + (size(op_matrix,1)-pod.nmodes)*s[end])
     pod.rbasis .= u
     nothing
 end
 
 function Base.show(io::IO,pod::POD)
     print(io,"POD \n")
     print(io,"Reduction Order = ",pod.nmodes,"\n")
-    print(io,"Snapshot size = (", length(pod.snapshots),",",length(pod.snapshots[1]),")\n")
-    print(io,"Energy = ", pod.energy,"\n")
+    print(io,"Snapshot size = (", size(pod.snapshots,1),",",size(pod.snapshots[1],2),")\n")
+    print(io,"Relative Energy = ", pod.renergy,"\n")
 end

test/DataReduction.jl

@@ -1,33 +1,43 @@
 #--------- Data Reduction -----------------#
 
-@testset "POD - Attractor Test" begin
-    function lorenz_prob()
-        function lorenz!(du,u,p,t)
-            du[1] = p[1]*(u[2]-u[1])
-            du[2] = u[1]*(p[2]-u[3]) - u[2]
-            du[3] = u[1]*u[2] - p[3]*u[3]
-        end
-
-        u0 = [1,0,0]
-        p = [10,28,8/3]
-        tspan = (0,100)
-        prob = ODEProblem(lorenz!,u0,tspan,p)
-        sol = solve(prob,Tsit5())
-        sol
+function lorenz_prob()
+    function lorenz!(du,u,p,t)
+        du[1] = p[1]*(u[2]-u[1])
+        du[2] = u[1]*(p[2]-u[3]) - u[2]
+        du[3] = u[1]*u[2] - p[3]*u[3]
     end
 
+    u0 = [1,0,0]
+    p = [10,28,8/3]
+    tspan = (0,100)
+    prob = ODEProblem(lorenz!,u0,tspan,p)
+    sol = solve(prob,Tsit5())
+    sol
+end
+
+@testset "POD-Utils" begin
+    solution = lorenz_prob()
+    VoV = solution.u
+    M = matricize(VoV)
+    @test size(M,1) == length(VoV[1]) # Parameters
+    @test size(M,2) == length(VoV) # Time
+end
+
+@testset "POD - Attractor Test" begin
     sol = lorenz_prob()
-    solution = Matrix(reduce(hcat,sol.u)')
+    solution = Array(sol)
 
     order = 2
     solver = SVD()
     reducer = POD(solution,order)
     reduce!(reducer,solver)
 
+    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.energy > 0.9
+    @test reducer.renergy > 0.9
 
     order = 2
     solver = TSVD()
@@ -43,6 +53,7 @@
     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)