Merge pull request #1 from dynamic-queries/Main

Main

Author: dynamic-queries

Project.toml

@@ -5,10 +5,7 @@ version = "0.1.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
-MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
 RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46"
 TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 

examples/POD/lorenz.jl

@@ -18,8 +18,8 @@ function lorenz_prob()
 end
 
 sol = lorenz_prob()
-solution = Matrix(reduce(hcat,sol.u)')
-plot(solution[:,1],solution[:,2],solution[:,3])
+solution = Array(sol)
+plot(solution[1,:],solution[2,:],solution[3,:])
 # savefig("lorenz_attractor.png")
 
 ## Two way POD

src/DataReduction/DiffusionMaps.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

src/DataReduction/POD.jl

@@ -1,5 +1,5 @@
 function errorhandle(data::Matrix{FT},modes::IT) where {FT,IT}
     @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 [1,$(s)]."
+    @assert (modes>0)&(modes<s) "Number of modes should be in {1,2,...,$(s-1)}."
 end

src/DataReduction/VAE.jl

@@ -8,13 +8,10 @@ module ModelOrderReduction
     using RandomizedLinAlg
 
     include("DataReduction/POD.jl")
-    include("DataReduction/DiffusionMaps.jl")
-    include("DataReduction/VAE.jl")
 
     export SVD, TSVD, RSVD
     export POD, reduce!, matricize
 #========================Model Reduction========================================#
-    include("ModelReduction/LiftAndLearn.jl")
 
 #===============================================================================#
 end

src/ErrorHandle.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)