Skip to content

Commit dfb39ba

Browse files
BoundaryValueProblemsBoundaryValueProblems
committed
Revamped HGLET.jl by adding HGLET_BasisBasis and HGLET_Analysis for users, and modified other files.
1 parent 84e7123 commit dfb39ba

File tree

3 files changed

+137
-38
lines changed

3 files changed

+137
-38
lines changed

‎src/HGLET.jl

Lines changed: 130 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -7,26 +7,26 @@ using ..GraphSignal, ..GraphPartition, ..BasisSpecification, ..GHWT, SparseArray
77
include("common.jl")
88

99
export HGLET_Synthesis, HGLET_jkl, HGLET_Analysis, HGLET_Analysis_All,
10-
HGLET_GHWT_BestBasis, HGLET_GHWT_Synthesis, HGLET_GHWT_BestBasis_minrelerror
10+
HGLET_BestBasis, HGLET_GHWT_BestBasis, HGLET_GHWT_Synthesis, HGLET_GHWT_BestBasis_minrelerror
1111

1212

1313
"""
1414
function HGLET_Synthesis(dvec::Vector{Float64}, GP::GraphPart, BS::BasisSpec,
15-
G::GraphSig; method::Symbol = :L)
15+
G::GraphSig; gltype::Symbol = :L)
1616
1717
### Input Arguments
1818
* `dvec`: the expansion coefficients corresponding to the chosen basis
1919
* `GP`: a GraphPart object
2020
* `BS`: a BasisSpec object
2121
* `G`: a GraphSig object
22-
* `method`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
22+
* `gltype`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
2323
2424
### Output Argument
2525
* `f`: the reconstructed signal
2626
* `GS`: the reconstructed GraphSig object
2727
"""
2828
function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
29-
G::GraphSig; method::Symbol = :L)
29+
G::GraphSig; gltype::Symbol = :L)
3030
# Preliminaries
3131
W = G.W
3232
jmax = size(GP.rs,2)
@@ -63,12 +63,12 @@ function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
6363
normalizep = true
6464
end
6565
# compute the GL eigenvectors via svd.
66-
if ( method == :Lrw || method == :Lsym ) && normalizep
66+
if ( gltype == :Lrw || gltype == :Lsym ) && normalizep
6767
v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
6868
else
6969
v,_,_ = svd(D - Matrix(W_temp))
70-
if method != :L
71-
@warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(method))
70+
if gltype != :L
71+
@warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(gltype))
7272
end
7373
end
7474
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,
9090
end
9191

9292
# reconstruct the signal
93-
if method == :Lrw && normalizep
93+
if gltype == :Lrw && normalizep
9494
f[rs1:rs3-1,:] = D2*v*dmatrix[rs1:rs3-1,j,:]
9595
else
9696
f[rs1:rs3-1,:] = v*dmatrix[rs1:rs3-1,j,:]
@@ -138,7 +138,7 @@ function HGLET_jkl(GP::GraphPart, drow::Int, dcol::Int)
138138
end
139139

140140
"""
141-
function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
141+
function HGLET_Analysis(G::GraphSig, GP::GraphPart, gltype::Symbol = :L)
142142
143143
For a GraphSig object 'G', generate the matrix of HGLET expansion
144144
coefficients corresponding to the eigenvectors of specified version of
@@ -147,12 +147,12 @@ function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
147147
### Input Arguments
148148
* `G`: a GraphSig object
149149
* `GP`: a GraphPart object
150-
* `method`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
150+
* `gltype`: :L, :Lrw, or :Lsym, indicating which eigenvectors are used
151151
152152
### Output Argument
153153
* `dmatrix`: the matrix of expansion coefficients using the specific GL matrix
154154
"""
155-
function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
155+
function HGLET_Analysis(G::GraphSig, GP::GraphPart, gltype::Symbol = :L)
156156
# Preliminaries
157157
W = G.W
158158
ind = GP.ind
@@ -189,12 +189,12 @@ function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
189189
end
190190

191191
# compute the GL eigenvectors via svd
192-
if ( method == :Lrw || method == :Lsym ) && normalizep
192+
if ( gltype == :Lrw || gltype == :Lsym ) && normalizep
193193
v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
194194
else
195195
v,_,_ = svd(D - Matrix(W_temp))
196-
if method != :L
197-
@warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(method))
196+
if gltype != :L
197+
@warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(gltype))
198198
end
199199
end
200200
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)
216216
end
217217

218218
# obtain the expansion coeffcients
219-
if method == :Lrw && normalizep
219+
if gltype == :Lrw && normalizep
220220
dmatrix[rs1:rs3-1,j,:] = v'*(D.^0.5)*G.f[indrs,:]
221221
else
222222
dmatrix[rs1:rs3-1,j,:] = v'*G.f[indrs,:]
@@ -247,6 +247,8 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
247247
W = G.W
248248
ind = GP.ind
249249
rs = GP.rs
250+
# This is the type of graph Laplacian matrix used for graph partitioning;
251+
# Not the eigenvector dictionary computed for the HGLET.
250252
method = GP.method
251253
N = size(G.W,1)
252254
jmax = size(rs,2)
@@ -377,13 +379,95 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
377379
end
378380

379381

382+
"""
383+
function HGLET_BestBasis(GP::GraphPart;
384+
dmatrix::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
385+
gltype::Symbol = :L, cfspec::Any = 0.1, flatten::Any = 1)
386+
387+
Select the HGLET best basis from the input matrix of expansion coefficients
388+
389+
### Input Arguments
390+
* `GP`: a GraphPart object
391+
* `dmatrix`: the matrix of HGLET expansion coefficients
392+
* `gltype`: the type of graph Laplacian matrix used for HGLET dictionary
393+
(default = :L)
394+
* `cfspec`: the cost functional specification to be used: see utils.jl
395+
for the detail. If it's a number, say, p, then the l^p norm is used.
396+
If it's a function, then that function is used. Default is 0.1, i.e., l^0.1
397+
* `flatten`: the method for flattening vector-valued data to scalar-valued data;
398+
If it's a number, say, p, then the sum of the l^p norm is computed.
399+
There are many other options, such as :ash, :entropy, etc. See utils.jl
400+
for the detail. Default is 1, i.e, the sum of the l^1 norm of the coefs.
401+
402+
### Output Argument
403+
* `dvec`: the vector of expansion coefficients corresponding to the bestbasis
404+
* `BS`: a BasisSpec object which specifies the best basis
405+
"""
406+
function HGLET_BestBasis(GP::GraphPart;
407+
dmatrix::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
408+
gltype::Symbol = :L, cfspec::Any = 0.1, flatten::Any = 1)
409+
410+
# the cost functional to be used
411+
costfun = cost_functional(cfspec)
412+
413+
# constants and dmatrix cleanup
414+
N, jmax, fcols = size(dmatrix)
415+
dmatrix[abs.(dmatrix) .< 10^2*eps()] .= 0
416+
417+
# flatten dmatrix for more than one input signals
418+
if fcols > 1
419+
dmatrix0 = deepcopy(dmatrix)
420+
dmatrix = dropdims(dmatrix_flatten(dmatrix, flatten), dims = 3)
421+
end
422+
423+
#
424+
# Find the HGLET best basis now!
425+
#
426+
427+
# allocate/initialize ==> order matters here
428+
dvec = dmatrix[:,jmax]
429+
levlist = jmax*ones(Int,N)
430+
431+
# set the tolerance
432+
tol = 10^4*eps()
433+
434+
# perform the basis search
435+
for j = jmax:-1:1
436+
regioncount = count(!iszero, GP.rs[:,j]) - 1
437+
for r = 1:regioncount
438+
indr = GP.rs[r,j]:(GP.rs[r+1,j]-1)
439+
### compute the cost of the current best basis
440+
costBB = costfun(dvec[indr])
441+
442+
### compute the cost of the HGLET coefficients
443+
costNew = costfun(dmatrix[indr,j])
444+
# change the best basis if the new cost is less expensive
445+
if costBB >= costNew - tol
446+
costBB, dvec[indr], levlist[indr] = BBchange(costNew, dmatrix[indr,j], j)
447+
end
448+
end
449+
end
450+
451+
levlist = collect(enumerate(levlist))
452+
453+
BS = BasisSpec(levlist, c2f = true, description = "HGLET Best Basis")
454+
#levlist2levlengths!(GP, BS)
455+
456+
# if we flattened dmatrix, then "unflatten" the expansion coefficients
457+
if fcols > 1
458+
dvec = dmatrix2dvec(dmatrix0, GP, BS)
459+
end
460+
461+
return dvec, BS
462+
end
463+
380464
"""
381465
function HGLET_GHWT_BestBasis(GP::GraphPart;
382466
dmatrixH::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
383467
dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
384468
dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
385469
dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
386-
costfun::Any = 0.1,flatten::Any = 1)
470+
cfspec::Any = 0.1,flatten::Any = 1)
387471
388472
Select the best basis from several matrices of expansion coefficients
389473
@@ -393,12 +477,17 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
393477
* `dmatrixHsym`: the matrix of HGLET expansion coefficients ==> eigenvectors of Lsym
394478
* `dmatrixG`: the matrix of GHWT expansion coefficients
395479
* `GP`: a GraphPart object
396-
* `costfun`: the cost functional to be used
397-
* `flatten`: the method for flattening vector-valued data to scalar-valued data
480+
* `cfspec`: the cost functional specification to be used: see utils.jl
481+
for the detail. If it's a number, say, p, then the l^p norm is used.
482+
If it's a function, then that function is used. Default is 0.1, i.e., l^0.1
483+
* `flatten`: the method for flattening vector-valued data to scalar-valued data;
484+
If it's a number, say, p, then the sum of the l^p norm is computed.
485+
There are many other options, such as :ash, :entropy, etc. See utils.jl
486+
for the detail. Default is 1, i.e, the sum of the l^1 norm of the coefs.
398487
399488
### Output Argument
400489
* `dvec`: the vector of expansion coefficients corresponding to the bestbasis
401-
* `BS`: a BasisSpec object which specifies the best-basis
490+
* `BS`: a BasisSpec object which specifies the best basis
402491
* `trans`: specifies which transform was used for that portion of the signal:
403492
00 = HGLET with L
404493
01 = HGLET with Lrw
@@ -410,15 +499,15 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
410499
dmatrixHrw::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
411500
dmatrixHsym::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
412501
dmatrixG::Array{Float64,3} = Array{Float64,3}(undef,0,0,0),
413-
costfun::Any = 0.1,flatten::Any = 1)
502+
cfspec::Any = 0.1,flatten::Any = 1)
414503
# specify transform codes
415504
transHsym = [true false]
416505
transG = [true true]
417506
transHrw = [false true]
418507
transH = [false false]
419508

420509
# the cost functional to be used
421-
costfun = cost_functional(costfun)
510+
costfun = cost_functional(cfspec)
422511

423512
# constants and dmatrix cleanup
424513
if !isempty(dmatrixHsym)
@@ -442,23 +531,23 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
442531
if fcols > 1
443532
if !isempty(dmatrixHsym)
444533
dmatrix0Hsym = deepcopy(dmatrixHsym)
445-
dmatrixdHsym = dropdims(dmatrix_flatten(dmatrixHsym,flatten),dims = 1)
534+
dmatrixHsym = dropdims(dmatrix_flatten(dmatrixHsym,flatten),dims = 3)
446535
end
447536
if !isempty(dmatrixG)
448537
dmatrix0G = deepcopy(dmatrixG)
449-
dmatrixG = dropdims(dmatrix_flatten(dmatrixG,flatten),dims = 1)
538+
dmatrixG = dropdims(dmatrix_flatten(dmatrixG,flatten),dims = 3)
450539
end
451540
if !isempty(dmatrixHrw)
452541
dmatrix0Hrw = deepcopy(dmatrixHrw)
453-
dmatrixHrw = dropdims(dmatrix_flatten(dmatrixHrw,flatten),dims = 1)
542+
dmatrixHrw = dropdims(dmatrix_flatten(dmatrixHrw,flatten),dims = 3)
454543
end
455544
if !isempty(dmatrixH)
456545
dmatrix0H = deepcopy(dmatrixH)
457-
dmatrixH = dropdims(dmatrix_flatten(dmatrixH,flatten),dims = 1)
546+
dmatrixH = dropdims(dmatrix_flatten(dmatrixH,flatten),dims = 3)
458547
end
459548
end
460549

461-
# Find the HGLET/GHWT best-basis
550+
# Find the HGLET/GHWT best basis
462551

463552
# allocate/initialize ==> order matters here
464553
if !isempty(dmatrixHsym)
@@ -508,7 +597,7 @@ function HGLET_GHWT_BestBasis(GP::GraphPart;
508597
end
509598
end
510599

511-
### compute the cost of the GHWT coefficients
600+
### compute the cost of the HGLET-Lrw coefficients
512601
if !isempty(dmatrixHrw)
513602
costNew = costfun(dmatrixHrw[indr,j])
514603
# change the best basis if the new cost is less expensive
@@ -569,6 +658,16 @@ return dvec, BS, transfull
569658

570659
end
571660

661+
function BBchange(costNew::Float64, dvec::Vector{Float64}, j::Int)
662+
#change to the new best basis
663+
costBB = costNew
664+
665+
n = length(dvec)
666+
667+
levlist = fill(j, n)
668+
669+
return costBB, dvec, levlist
670+
end
572671

573672
function BBchange(costNew::Float64, dvec::Vector{Float64}, j::Int, trans::Array{Bool,2})
574673
#change to the new best basis
@@ -619,8 +718,8 @@ function HGLET_GHWT_Synthesis(dvec::Matrix{Float64},GP::GraphPart,BS::BasisSpec,
619718
# Synthesize using the transforms separately
620719

621720
fH,_ = HGLET_Synthesis(dvecH, GP, BS, G)
622-
fHrw,_ = HGLET_Synthesis(dvecHrw, GP, BS, G, method = :Lrw)
623-
fHsym,_ = HGLET_Synthesis(dvecHsym, GP, BS, G, method = :Lsym)
721+
fHrw,_ = HGLET_Synthesis(dvecHrw, GP, BS, G, gltype = :Lrw)
722+
fHsym,_ = HGLET_Synthesis(dvecHsym, GP, BS, G, gltype = :Lsym)
624723
fG = ghwt_synthesis(dvecG, GP, BS)
625724

626725
f = fH + fHrw + fHsym + fG
@@ -656,7 +755,7 @@ function HGLET_GHWT_BestBasis_minrelerror(GP::GraphPart, G::GraphSig;
656755
657756
### Output argument:
658757
* `dvec`: the vector of expansion coefficients corresponding to the bestbasis
659-
* `BS`: a BasisSpec object which specifies the best-basis
758+
* `BS`: a BasisSpec object which specifies the best basis
660759
* `trans`: specifies which transform was used for that portion of the signal
661760
00 = HGLET with L
662761
01 = HGLET with Lrw
@@ -684,7 +783,7 @@ function HGLET_GHWT_BestBasis_minrelerror(GP::GraphPart, G::GraphSig;
684783
dvec_temp, BS_temp, trans_temp = HGLET_GHWT_BestBasis(GP,
685784
dmatrixH = dmatrixH, dmatrixG = dmatrixG,
686785
dmatrixHrw = dmatrixHrw, dmatrixHsym = dmatrixHsym,
687-
costfun = tau_temp)
786+
cfspec = tau_temp)
688787

689788
# check whether any HGLET Lrw basis vectors are in the best basis
690789
orthbasis = true

‎src/LP-HGLET.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -182,26 +182,26 @@ function standardize_eigenvector_signs!(vec)
182182
end
183183

184184
"""
185-
HGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L)
185+
HGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L)
186186
187187
assemble the whole HGLET dictionary
188188
189189
### Input Arguments
190190
* `GP`: a GraphPart object
191191
* `G`: a GraphSig object
192-
* `method`: `:L` or `:Lsym`
192+
* `gltype`: `:L` or `:Lsym`
193193
194194
### Output Argument
195195
* `dictionary`: the HGLET dictionary
196196
197197
"""
198-
function HGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L)
198+
function HGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L)
199199
N = size(G.W, 1)
200200
jmax = size(GP.rs, 2)
201201
dictionary = zeros(N, jmax, N)
202202
for j = 1:jmax
203203
BS = BasisSpec(collect(enumerate(j * ones(Int, N))))
204-
dictionary[:, j, :] = HGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; method = method)[1]'
204+
dictionary[:, j, :] = HGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; gltype = method)[1]'
205205
end
206206
return dictionary
207207
end

‎src/helpers.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -306,7 +306,7 @@ function getall_expansioncoeffs(G_Sig::GraphSig, GP_star::GraphPart, VM_NGWP::Ar
306306
############# plain Laplacian eigenvectors coefficients
307307
GP = partition_tree_fiedler(G_Sig)
308308
dmatrixH = HGLET_Analysis_All(G_Sig, GP)[1]
309-
dvec_hglet, BS_hglet, _ = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, costfun = 1)
309+
dvec_hglet, BS_hglet, _ = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, cfspec = 1)
310310
dmatrix = ghwt_analysis!(G_Sig, GP = GP)
311311
############# Haar
312312
BS_haar = bs_haar(GP)
@@ -628,10 +628,10 @@ function getall_expansioncoeffs2(G_Sig::GraphSig, GP_star::GraphPart,
628628
############# HGLET
629629
GP = partition_tree_fiedler(G_Sig)
630630
dmatrixH, _, dmatrixHsym = HGLET_Analysis_All(G_Sig, GP)
631-
dvec_hglet, BS_hglet, trans_hglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, dmatrixHsym = dmatrixHsym, costfun = 1)
631+
dvec_hglet, BS_hglet, trans_hglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH, dmatrixHsym = dmatrixHsym, cfspec = 1)
632632
############# LP-HGLET
633633
dmatrixsH, dmatrixsHsym = LPHGLET_Analysis_All(G_Sig, GP; ϵ = 0.3)
634-
dvec_lphglet, BS_lphglet, trans_lphglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixsH, dmatrixHsym = dmatrixsHsym, costfun = 1)
634+
dvec_lphglet, BS_lphglet, trans_lphglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixsH, dmatrixHsym = dmatrixsHsym, cfspec = 1)
635635
############# GHWT dictionaries
636636
dmatrix = ghwt_analysis!(G_Sig, GP = GP)
637637
############# Haar

0 commit comments

Comments
 (0)