Revamped HGLET.jl by adding HGLET_BasisBasis and HGLET_Analysis for users, and modified other files.

Author: BoundaryValueProblems

src/HGLET.jl

@@ -7,26 +7,26 @@ using ..GraphSignal, ..GraphPartition, ..BasisSpecification, ..GHWT, SparseArray
 include("common.jl")
 
 export HGLET_Synthesis, HGLET_jkl, HGLET_Analysis, HGLET_Analysis_All,
-HGLET_GHWT_BestBasis, HGLET_GHWT_Synthesis, HGLET_GHWT_BestBasis_minrelerror
+HGLET_BestBasis, HGLET_GHWT_BestBasis, HGLET_GHWT_Synthesis, HGLET_GHWT_BestBasis_minrelerror
 
 
 """
 function HGLET_Synthesis(dvec::Vector{Float64}, GP::GraphPart, BS::BasisSpec,
-                         G::GraphSig; method::Symbol = :L)
+                         G::GraphSig; gltype::Symbol = :L)
 
 ### Input Arguments
  * `dvec`: the expansion coefficients corresponding to the chosen basis
  * `GP`: a GraphPart object
  * `BS`: a BasisSpec object
  * `G`: a GraphSig object
- * `method`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
+ * `gltype`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
 
 ### Output Argument
  * `f`: the reconstructed signal
  * `GS`: the reconstructed GraphSig object
 """
 function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
-                         G::GraphSig; method::Symbol = :L)
+                         G::GraphSig; gltype::Symbol = :L)
     # Preliminaries
     W = G.W
     jmax = size(GP.rs,2)
@@ -63,12 +63,12 @@ function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
                         normalizep = true
                     end
                     # compute the GL eigenvectors via svd.
-                    if ( method == :Lrw || method == :Lsym ) && normalizep
+                    if ( gltype == :Lrw || gltype == :Lsym ) && normalizep
                         v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
                     else
                         v,_,_ = svd(D - Matrix(W_temp))
-                        if method != :L
-                            @warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(method))
+                        if gltype != :L
+                            @warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(gltype))
                         end
                     end
                     v = v[:,end:-1:1] # reorder the ev's in the decreasing ew's
@@ -90,7 +90,7 @@ function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
                     end
 
                     # reconstruct the signal
-                    if method == :Lrw && normalizep
+                    if gltype == :Lrw && normalizep
                         f[rs1:rs3-1,:] = D2*v*dmatrix[rs1:rs3-1,j,:]
                     else
                         f[rs1:rs3-1,:] = v*dmatrix[rs1:rs3-1,j,:]
@@ -138,7 +138,7 @@ function HGLET_jkl(GP::GraphPart, drow::Int, dcol::Int)
 end
 
 """
-function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
+function HGLET_Analysis(G::GraphSig, GP::GraphPart, gltype::Symbol = :L)
 
     For a GraphSig object 'G', generate the matrix of HGLET expansion
     coefficients corresponding to the eigenvectors of specified version of
@@ -147,12 +147,12 @@ function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
 ### Input Arguments
  * `G`:  a GraphSig object
  * `GP`: a GraphPart object
- * `method`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
+ * `gltype`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
 
 ### Output Argument
  * `dmatrix`: the matrix of expansion coefficients using the specific GL matrix
 """
-function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
+function HGLET_Analysis(G::GraphSig, GP::GraphPart, gltype::Symbol = :L)
     # Preliminaries
     W = G.W
     ind = GP.ind
@@ -189,12 +189,12 @@ function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
                 end
 
                 # compute the GL eigenvectors via svd
-                if ( method == :Lrw || method == :Lsym ) && normalizep
+                if ( gltype == :Lrw || gltype == :Lsym ) && normalizep
                     v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
                 else
                     v,_,_ = svd(D - Matrix(W_temp))
-                    if method != :L
-                        @warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(method))
+                    if gltype != :L
+                        @warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(gltype))
                     end
                 end
                 v = v[:,end:-1:1] # reorder the ev's in the decreasing ew's
@@ -216,7 +216,7 @@ function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
                 end
 
                 # obtain the expansion coeffcients
-                if method == :Lrw && normalizep
+                if gltype == :Lrw && normalizep
                     dmatrix[rs1:rs3-1,j,:] = v'*(D.^0.5)*G.f[indrs,:]
                 else
                     dmatrix[rs1:rs3-1,j,:] = v'*G.f[indrs,:]
@@ -247,6 +247,8 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
     W = G.W
     ind = GP.ind
     rs = GP.rs
+    # This is the type of graph Laplacian matrix used for graph partitioning;
+    # Not the eigenvector dictionary computed for the HGLET.
     method = GP.method
     N = size(G.W,1)
     jmax = size(rs,2)
@@ -377,13 +379,95 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
 end
 
 
+"""
+function HGLET_BestBasis(GP::GraphPart; 
+    dmatrix::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
+    gltype::Symbol = :L, cfspec::Any = 0.1, flatten::Any = 1)
+
+    Select the HGLET best basis from the input matrix of expansion coefficients
+
+### Input Arguments
+* `GP`:      a GraphPart object
+* `dmatrix`: the matrix of HGLET expansion coefficients 
+* `gltype`:  the type of graph Laplacian matrix used for HGLET dictionary
+    (default = :L)
+* `cfspec`: the cost functional specification to be used: see utils.jl
+    for the detail. If it's a number, say, p, then the l^p norm is used.
+    If it's a function, then that function is used. Default is 0.1, i.e., l^0.1
+* `flatten`: the method for flattening vector-valued data to scalar-valued data;
+    If it's a number, say, p, then the sum of the l^p norm is computed.
+    There are many other options, such as :ash, :entropy, etc. See utils.jl
+    for the detail. Default is 1, i.e, the sum of the l^1 norm of the coefs.
+
+### Output Argument
+* `dvec`: the vector of expansion coefficients corresponding to the bestbasis
+* `BS`:   a BasisSpec object which specifies the best basis
+"""
+function HGLET_BestBasis(GP::GraphPart; 
+    dmatrix::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
+    gltype::Symbol = :L, cfspec::Any = 0.1, flatten::Any = 1)
+
+    # the cost functional to be used
+    costfun = cost_functional(cfspec)
+
+    # constants and dmatrix cleanup
+    N, jmax, fcols = size(dmatrix)
+    dmatrix[abs.(dmatrix) .< 10^2*eps()] .= 0
+
+    # flatten dmatrix for more than one input signals
+    if fcols > 1
+        dmatrix0 = deepcopy(dmatrix)
+        dmatrix = dropdims(dmatrix_flatten(dmatrix, flatten), dims = 3)
+    end
+
+    #
+    # Find the HGLET best basis now!
+    #
+
+    # allocate/initialize ==> order matters here
+    dvec = dmatrix[:,jmax]
+    levlist = jmax*ones(Int,N)
+
+    # set the tolerance
+    tol = 10^4*eps()
+
+    # perform the basis search
+    for j = jmax:-1:1
+        regioncount = count(!iszero, GP.rs[:,j]) - 1
+        for r = 1:regioncount
+            indr = GP.rs[r,j]:(GP.rs[r+1,j]-1)
+            ### compute the cost of the current best basis
+            costBB = costfun(dvec[indr])
+
+            ### compute the cost of the HGLET coefficients
+            costNew = costfun(dmatrix[indr,j])
+            # change the best basis if the new cost is less expensive
+            if costBB >= costNew - tol
+                costBB, dvec[indr], levlist[indr] = BBchange(costNew, dmatrix[indr,j], j)
+            end
+        end
+    end
+
+    levlist = collect(enumerate(levlist))
+
+    BS = BasisSpec(levlist, c2f = true, description = "HGLET Best Basis")
+    #levlist2levlengths!(GP, BS)
+
+    # if we flattened dmatrix, then "unflatten" the expansion coefficients
+    if fcols > 1
+        dvec = dmatrix2dvec(dmatrix0, GP, BS)
+    end
+
+    return dvec, BS
+end
+
 """
 function HGLET_GHWT_BestBasis(GP::GraphPart; 
     dmatrixH::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
     dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0), 
     dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
     dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0), 
-    costfun::Any = 0.1,flatten::Any = 1)
+    cfspec::Any = 0.1,flatten::Any = 1)
 
     Select the best basis from several matrices of expansion coefficients
 
@@ -393,12 +477,17 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
 * `dmatrixHsym`: the matrix of HGLET expansion coefficients ==> eigenvectors of Lsym
 * `dmatrixG`:    the matrix of GHWT expansion coefficients
 * `GP`:          a GraphPart object
-* `costfun`:     the cost functional to be used
-* `flatten`:     the method for flattening vector-valued data to scalar-valued data
+* `cfspec`: the cost functional specification to be used: see utils.jl
+    for the detail. If it's a number, say, p, then the l^p norm is used.
+    If it's a function, then that function is used. Default is 0.1, i.e., l^0.1
+* `flatten`: the method for flattening vector-valued data to scalar-valued data;
+    If it's a number, say, p, then the sum of the l^p norm is computed.
+    There are many other options, such as :ash, :entropy, etc. See utils.jl
+    for the detail. Default is 1, i.e, the sum of the l^1 norm of the coefs.
 
 ### Output Argument
 * `dvec`:     the vector of expansion coefficients corresponding to the bestbasis
-* `BS`:       a BasisSpec object which specifies the best-basis
+* `BS`:       a BasisSpec object which specifies the best basis
 * `trans`:    specifies which transform was used for that portion of the signal:
     00 = HGLET with L
     01 = HGLET with Lrw
@@ -410,15 +499,15 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
     dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
     dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
     dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0), 
-    costfun::Any = 0.1,flatten::Any = 1)
+    cfspec::Any = 0.1,flatten::Any = 1)
     # specify transform codes
     transHsym = [true false]
     transG = [true true]
     transHrw = [false true]
     transH = [false false]
 
     # the cost functional to be used
-    costfun = cost_functional(costfun)
+    costfun = cost_functional(cfspec)
 
     # constants and dmatrix cleanup
     if !isempty(dmatrixHsym)
@@ -442,23 +531,23 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
     if fcols > 1
         if !isempty(dmatrixHsym)
             dmatrix0Hsym = deepcopy(dmatrixHsym)
-            dmatrixdHsym = dropdims(dmatrix_flatten(dmatrixHsym,flatten),dims = 1)
+            dmatrixHsym = dropdims(dmatrix_flatten(dmatrixHsym,flatten),dims = 3)
         end
         if !isempty(dmatrixG)
             dmatrix0G = deepcopy(dmatrixG)
-            dmatrixG = dropdims(dmatrix_flatten(dmatrixG,flatten),dims = 1)
+            dmatrixG = dropdims(dmatrix_flatten(dmatrixG,flatten),dims = 3)
         end
         if !isempty(dmatrixHrw)
             dmatrix0Hrw = deepcopy(dmatrixHrw)
-            dmatrixHrw = dropdims(dmatrix_flatten(dmatrixHrw,flatten),dims = 1)
+            dmatrixHrw = dropdims(dmatrix_flatten(dmatrixHrw,flatten),dims = 3)
         end
         if !isempty(dmatrixH)
             dmatrix0H = deepcopy(dmatrixH)
-            dmatrixH = dropdims(dmatrix_flatten(dmatrixH,flatten),dims = 1)
+            dmatrixH = dropdims(dmatrix_flatten(dmatrixH,flatten),dims = 3)
         end
     end
 
-    # Find the HGLET/GHWT best-basis
+    # Find the HGLET/GHWT best basis
 
     # allocate/initialize ==> order matters here
     if !isempty(dmatrixHsym)
@@ -508,7 +597,7 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
                 end
             end
 
-            ### compute the cost of the GHWT coefficients
+            ### compute the cost of the HGLET-Lrw coefficients
             if !isempty(dmatrixHrw)
                 costNew = costfun(dmatrixHrw[indr,j])
                 # change the best basis if the new cost is less expensive
@@ -569,6 +658,16 @@ return dvec, BS, transfull
 
 end
 
+function BBchange(costNew::Float64, dvec::Vector{Float64}, j::Int)
+    #change to the new best basis
+    costBB = costNew
+
+    n = length(dvec)
+
+    levlist = fill(j, n)
+
+    return costBB, dvec, levlist
+end
 
 function BBchange(costNew::Float64, dvec::Vector{Float64}, j::Int, trans::Array{Bool,2})
     #change to the new best basis
@@ -619,8 +718,8 @@ function HGLET_GHWT_Synthesis(dvec::Matrix{Float64},GP::GraphPart,BS::BasisSpec,
     # Synthesize using the transforms separately
 
     fH,_ = HGLET_Synthesis(dvecH, GP, BS, G)
-    fHrw,_ = HGLET_Synthesis(dvecHrw, GP, BS, G, method = :Lrw)
-    fHsym,_ = HGLET_Synthesis(dvecHsym, GP, BS, G, method = :Lsym)
+    fHrw,_ = HGLET_Synthesis(dvecHrw, GP, BS, G, gltype = :Lrw)
+    fHsym,_ = HGLET_Synthesis(dvecHsym, GP, BS, G, gltype = :Lsym)
     fG = ghwt_synthesis(dvecG, GP, BS)
 
     f = fH + fHrw + fHsym + fG
@@ -656,7 +755,7 @@ function HGLET_GHWT_BestBasis_minrelerror(GP::GraphPart, G::GraphSig;
 
 ### Output argument:
 * `dvec`: the vector of expansion coefficients corresponding to the bestbasis
-* `BS`: a BasisSpec object which specifies the best-basis
+* `BS`: a BasisSpec object which specifies the best basis
 * `trans`: specifies which transform was used for that portion of the signal
     00 = HGLET with L
     01 = HGLET with Lrw
@@ -684,7 +783,7 @@ function HGLET_GHWT_BestBasis_minrelerror(GP::GraphPart, G::GraphSig;
             dvec_temp, BS_temp, trans_temp = HGLET_GHWT_BestBasis(GP, 
                 dmatrixH = dmatrixH, dmatrixG = dmatrixG, 
                 dmatrixHrw = dmatrixHrw, dmatrixHsym = dmatrixHsym, 
-                costfun = tau_temp)
+                cfspec = tau_temp)
 
             # check whether any HGLET Lrw basis vectors are in the best basis
             orthbasis = true

src/LP-HGLET.jl

@@ -182,26 +182,26 @@ function standardize_eigenvector_signs!(vec)
 end
 
 """
-    HGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L)
+    HGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L)
 
 assemble the whole HGLET dictionary
 
 ### Input Arguments
 * `GP`: a GraphPart object
 * `G`:  a GraphSig object
-* `method`: `:L` or `:Lsym`
+* `gltype`: `:L` or `:Lsym`
 
 ### Output Argument
 * `dictionary`: the HGLET dictionary
 
 """
-function HGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L)
+function HGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L)
     N = size(G.W, 1)
     jmax = size(GP.rs, 2)
     dictionary = zeros(N, jmax, N)
     for j = 1:jmax
         BS = BasisSpec(collect(enumerate(j * ones(Int, N))))
-        dictionary[:, j, :] = HGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; method = method)[1]'
+        dictionary[:, j, :] = HGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; gltype = method)[1]'
     end
     return dictionary
 end

src/helpers.jl

@@ -306,7 +306,7 @@ function getall_expansioncoeffs(G_Sig::GraphSig, GP_star::GraphPart, VM_NGWP::Ar
     ############# plain Laplacian eigenvectors coefficients
     GP = partition_tree_fiedler(G_Sig)
     dmatrixH = HGLET_Analysis_All(G_Sig, GP)[1]
-    dvec_hglet, BS_hglet, _ = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, costfun = 1)
+    dvec_hglet, BS_hglet, _ = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, cfspec = 1)
     dmatrix = ghwt_analysis!(G_Sig, GP = GP)
     ############# Haar
     BS_haar = bs_haar(GP)
@@ -628,10 +628,10 @@ function getall_expansioncoeffs2(G_Sig::GraphSig, GP_star::GraphPart,
     ############# HGLET
     GP = partition_tree_fiedler(G_Sig)
     dmatrixH, _, dmatrixHsym = HGLET_Analysis_All(G_Sig, GP)
-    dvec_hglet, BS_hglet, trans_hglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, dmatrixHsym = dmatrixHsym, costfun = 1)
+    dvec_hglet, BS_hglet, trans_hglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, dmatrixHsym = dmatrixHsym, cfspec = 1)
     ############# LP-HGLET
     dmatrixsH, dmatrixsHsym = LPHGLET_Analysis_All(G_Sig, GP; ϵ = 0.3)
-    dvec_lphglet, BS_lphglet, trans_lphglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixsH, dmatrixHsym = dmatrixsHsym, costfun = 1)
+    dvec_lphglet, BS_lphglet, trans_lphglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixsH, dmatrixHsym = dmatrixsHsym, cfspec = 1)
     ############# GHWT dictionaries
     dmatrix = ghwt_analysis!(G_Sig, GP = GP)
     ############# Haar