Reduced Bases VoV->Matrix

Considering that we need the reduced bases in Matrix form down the line,
for POD based MOR methods, a change in the DS is made.

Author: Rahul Manavalan

Project.toml

@@ -5,7 +5,10 @@ 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

@@ -0,0 +1,40 @@
+using ModelOrderReduction
+using Plots
+using OrdinaryDiffEq
+
+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
+
+sol = lorenz_prob()
+solution = Matrix(reduce(hcat,sol.u)')
+plot(solution[:,1],solution[:,2],solution[:,3])
+# savefig("lorenz_attractor.png")
+
+## Two way POD
+reducer = POD(solution,2)
+reduce!(reducer,SVD())
+bases = reducer.rbasis
+plot(bases[:,1],bases[:,2],label="POD2")
+# savefig("pod2.png")
+
+
+## One way POD
+reducer = POD(solution,1)
+reduce!(reducer,SVD())
+bases = reducer.rbasis
+plot(bases[:,1],label="POD1")
+# savefig("pod1.png")
+plot(solution[:,3],label='z')
+# savefig("z.png")

src/DataReduction/POD.jl

@@ -1,23 +1,25 @@
-# Leverage PCA from MultivariateStats.jl ?
-
 function matricize(VoV::Vector{Vector{FT}}) where {FT}
     Matrix(reduce(hcat,VoV)')
 end
 
 mutable struct POD{FT,IT} <: AbstractDRProblem
     snapshots::Union{Vector{Vector{FT}},Matrix{FT}}
     nmodes::IT
-    rbasis::Vector{Vector{FT}}
+    rbasis::Matrix{FT}
     energy::FT
 
     function POD(snaps::Vector{Vector{FT}},nmodes::IT) where {FT,IT}
         errorhandle(matricize(snaps),nmodes)
-        new{eltype(snaps[1]),typeof(nmodes)}(snaps,nmodes)
+        new{eltype(snaps[1]),typeof(nmodes)}(snaps,nmodes,Array{FT,2}(undef,size(snaps,1),nmodes),FT(0))
     end
 
     function POD(snaps::Matrix{FT},nmodes::IT) where {FT,IT}
         errorhandle(snaps,nmodes)
-        new{eltype(snaps),typeof(nmodes)}(snaps,nmodes)
+        new{eltype(snaps),typeof(nmodes)}(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}
+        POD(Matrix(snaps),nmodes,Array{FT,2}(undef,size(snaps,1),nmodes),FT(0))
     end
 end
 
@@ -27,10 +29,9 @@ function reduce!(pod::POD{FT,IT},::SVD) where {FT,IT}
         op_matrix = matricize(pod.snapshots)
     end
     u,s,v = svd(op_matrix)
-    vr = v[:,1:pod.nmodes]
+    pod.rbasis .= u[:,1:pod.nmodes]
     sr = s[1:pod.nmodes]
     pod.energy = sum(sr)/sum(s)
-    pod.rbasis = [vr[:,i] for i=1:pod.nmodes]
     nothing
 end
 
@@ -41,7 +42,7 @@ function reduce!(pod::POD{FT,IT},::TSVD) where {FT,IT}
     end
     u,s,v = tsvd(op_matrix,pod.nmodes)
     pod.energy = NaN
-    pod.rbasis = [v[:,i] for i=1:pod.nmodes]
+    pod.rbasis .= u
     nothing
 end
 
@@ -52,7 +53,7 @@ function reduce!(pod::POD{FT,IT},::RSVD) where {FT,IT}
     end
     u,s,v = rsvd(op_matrix,pod.nmodes)
     pod.energy = NaN
-    pod.rbasis = [v[:,i] for i=1:pod.nmodes]
+    pod.rbasis .= u
     nothing
 end
 

src/ModelOrderReduction.jl

@@ -2,18 +2,19 @@ module ModelOrderReduction
 #========================Data Reduction=========================================#
     include("Types.jl")
     include("ErrorHandle.jl")
-    
+
     using LinearAlgebra
     using TSVD
     using RandomizedLinAlg
 
     include("DataReduction/POD.jl")
-    include("DataReduction/DifussionMaps.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/DataReduction/DifussionMaps.jl src/ModelReduction/LiftAndLearn.jlsrc/DataReduction/DifussionMaps.jl renamed to src/ModelReduction/LiftAndLearn.jl

@@ -1,5 +1,6 @@
 #--------- Data Reduction -----------------#
-@testset "POD" begin
+
+@testset "POD - Attractor Test" begin
     function lorenz_prob()
         function lorenz!(du,u,p,t)
             du[1] = p[1]*(u[2]-u[1])
@@ -24,8 +25,8 @@
     reduce!(reducer,solver)
 
     # Ad-hoc tests. To be checked with Chris.
-    @test length(reducer.rbasis) == reducer.nmodes
-    @test length(reducer.rbasis[1]) == size(solution,2)
+    @test size(reducer.rbasis,2) == reducer.nmodes
+    @test size(reducer.rbasis,1) == size(solution,1)
     @test reducer.energy > 0.9
 
     order = 2
@@ -34,15 +35,15 @@
     reduce!(reducer,solver)
 
     # Ad-hoc tests. To be checked with Chris.
-    @test length(reducer.rbasis) == reducer.nmodes
-    @test length(reducer.rbasis[1]) == size(solution,2)
+    @test size(reducer.rbasis,2) == reducer.nmodes
+    @test size(reducer.rbasis,1) == size(solution,1)
 
     order = 2
     solver = RSVD()
     reducer = POD(solution,order)
     reduce!(reducer,solver)
 
     # Ad-hoc tests. To be checked with Chris.
-    @test length(reducer.rbasis) == reducer.nmodes
-    @test length(reducer.rbasis[1]) == size(solution,2)
+    @test size(reducer.rbasis,2) == reducer.nmodes
+    @test size(reducer.rbasis,1) == size(solution,1)
 end