|
| 1 | +""" |
| 2 | + function LPHGLET_Synthesis(dvec::Vector{Float64}, GP::GraphPart, BS::BasisSpec, G::GraphSig; method::Symbol = :L, ϵ::Float64 = 0.3) |
| 3 | +
|
| 4 | +Perform Lapped-HGLET Synthesis transform |
| 5 | +
|
| 6 | +### Input Arguments |
| 7 | +* `dvec`: the expansion coefficients corresponding to the chosen basis |
| 8 | +* `GP`: a GraphPart object |
| 9 | +* `BS`: a BasisSpec object |
| 10 | +* `G`: a GraphSig object |
| 11 | +* `method`: :L or :Lsym, indicating which eigenvectors are used |
| 12 | +* `ϵ`: relative action bandwidth (default: 0.3) |
| 13 | +
|
| 14 | +### Output Argument |
| 15 | +* `f`: the reconstructed signal |
| 16 | +* `GS`: the reconstructed GraphSig object |
| 17 | +""" |
| 18 | +function LPHGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec, G::GraphSig; method::Symbol = :L, ϵ::Float64 = 0.3) |
| 19 | + # Preliminaries |
| 20 | + W = G.W |
| 21 | + inds = GP.inds |
| 22 | + rs = GP.rs |
| 23 | + N = size(W, 1) |
| 24 | + jmax = size(rs, 2) |
| 25 | + Uf = Matrix{Float64}(I, N, N) |
| 26 | + used_node = Set() |
| 27 | + |
| 28 | + # fill in the appropriate entries of dmatrix |
| 29 | + dmatrix = dvec2dmatrix(dvec, GP, BS) |
| 30 | + |
| 31 | + f = zeros(size(dmatrix[:, jmax, :])) |
| 32 | + |
| 33 | + |
| 34 | + # Perform the synthesis transform |
| 35 | + for j = 1:jmax |
| 36 | + regioncount = count(!iszero, rs[:,j]) - 1 |
| 37 | + # assemble orthogonal folding operator at level j - 1 |
| 38 | + keep_folding!(Uf, used_node, W, GP; ϵ = ϵ, j = j - 1) |
| 39 | + for r = 1:regioncount |
| 40 | + # indices of current region |
| 41 | + indr = rs[r, j]:(rs[r + 1, j] - 1) |
| 42 | + # indices of current region's nodes |
| 43 | + indrs = inds[indr, j] |
| 44 | + # number of nodes in current region |
| 45 | + n = length(indrs) |
| 46 | + |
| 47 | + # only proceed forward if coefficients do not exist |
| 48 | + if (j == jmax || count(!iszero, dmatrix[indr, j + 1, :]) == 0) && count(!iszero, dmatrix[indr, j, :]) > 0 |
| 49 | + # compute the eigenvectors |
| 50 | + W_temp = W[indrs,indrs] |
| 51 | + D_temp = sparse(Diagonal(dropdims(sum(W_temp, dims = 1), dims = 1))) |
| 52 | + if method == :L |
| 53 | + # compute the eigenvectors of L ==> svd(L) |
| 54 | + vec = svd(Matrix(D_temp - W_temp)).U |
| 55 | + elseif method == :Lsym |
| 56 | + # check if one can assemble the Lsym |
| 57 | + if minimum(sum(W[indrs, indrs], dims = 1)) > 10^3 * eps() |
| 58 | + ### eigenvectors of L_sym ==> svd(L_sym) |
| 59 | + D_temp_p = sparse(Diagonal(dropdims(sum(W_temp, dims = 1), dims = 1).^(-1/2))) |
| 60 | + vec = svd(Matrix(D_temp_p * (D_temp - W_temp) * D_temp_p)).U |
| 61 | + else |
| 62 | + ### eigenvectors of L ==> svd(L) |
| 63 | + vec = svd(Matrix(D_temp - W_temp)).U |
| 64 | + end |
| 65 | + end |
| 66 | + vec = vec[:, end:-1:1] |
| 67 | + |
| 68 | + |
| 69 | + # standardize the eigenvector signs |
| 70 | + standardize_eigenvector_signs!(vec) |
| 71 | + |
| 72 | + # construct unfolder operator custom to current region |
| 73 | + P = Uf[indrs, :]' |
| 74 | + |
| 75 | + # reconstruct the signal |
| 76 | + f += (P * vec) * dmatrix[indr, j, :] |
| 77 | + |
| 78 | + end |
| 79 | + end |
| 80 | + end |
| 81 | + |
| 82 | + # creat a GraphSig object with the reconstructed data |
| 83 | + GS = deepcopy(G) |
| 84 | + replace_data!(GS, f) |
| 85 | + |
| 86 | + return f, GS |
| 87 | +end |
| 88 | + |
| 89 | + |
| 90 | + |
| 91 | +""" |
| 92 | + function LPHGLET_Analysis_All(G::GraphSig, GP::GraphPart; ϵ::Float64 = 0.3) |
| 93 | +
|
| 94 | +For a GraphSig object 'G', generate the 2 matrices of Lapped-HGLET expansion coefficients |
| 95 | +corresponding to the eigenvectors of L and Lsym |
| 96 | +
|
| 97 | +### Input Arguments |
| 98 | +* `G`: a GraphSig object |
| 99 | +* `GP`: a GraphPart object |
| 100 | +* `ϵ`: relative action bandwidth (default: 0.3) |
| 101 | +
|
| 102 | +### Output Argument |
| 103 | +* `dmatrixlH`: the matrix of expansion coefficients for L |
| 104 | +* `dmatrixlHsym`: the matrix of expansion coefficients for Lsym |
| 105 | +* `GP`: a GraphPart object |
| 106 | +""" |
| 107 | +function LPHGLET_Analysis_All(G::GraphSig, GP::GraphPart; ϵ::Float64 = 0.3) |
| 108 | + # Preliminaries |
| 109 | + W = G.W |
| 110 | + inds = GP.inds |
| 111 | + rs = GP.rs |
| 112 | + N = size(W, 1) |
| 113 | + jmax = size(rs, 2) |
| 114 | + fcols = size(G.f, 2) |
| 115 | + Uf = Matrix{Float64}(I, N, N) |
| 116 | + used_node = Set() |
| 117 | + dmatrixlH = zeros(N, jmax, fcols) |
| 118 | + dmatrixlHsym = deepcopy(dmatrixlH) |
| 119 | + |
| 120 | + for j = 1:jmax |
| 121 | + regioncount = count(!iszero, rs[:,j]) - 1 |
| 122 | + # assemble orthogonal folding operator at level j - 1 |
| 123 | + keep_folding!(Uf, used_node, W, GP; ϵ = ϵ, j = j - 1) |
| 124 | + for r = 1:regioncount |
| 125 | + # indices of current region |
| 126 | + indr = rs[r, j]:(rs[r + 1, j] - 1) |
| 127 | + # indices of current region's nodes |
| 128 | + indrs = inds[indr, j] |
| 129 | + # number of nodes in current region |
| 130 | + n = length(indrs) |
| 131 | + |
| 132 | + # compute the eigenvectors |
| 133 | + W_temp = W[indrs,indrs] |
| 134 | + D_temp = sparse(Diagonal(dropdims(sum(W_temp, dims = 1), dims = 1))) |
| 135 | + ## eigenvectors of L ==> svd(L) |
| 136 | + vec = svd(Matrix(D_temp - W_temp)).U |
| 137 | + ## eigenvectors of L_sym ==> svd(L_sym) |
| 138 | + if minimum(sum(W[indrs, indrs], dims = 1)) > 10^3 * eps() |
| 139 | + ### eigenvectors of L_sym ==> svd(L_sym) |
| 140 | + D_temp_p = sparse(Diagonal(dropdims(sum(W_temp, dims = 1), dims = 1).^(-1/2))) |
| 141 | + vec_sym = svd(Matrix(D_temp_p * (D_temp - W_temp) * D_temp_p)).U |
| 142 | + else |
| 143 | + ### eigenvectors of L ==> svd(L) |
| 144 | + vec_sym = deepcopy(vec) |
| 145 | + end |
| 146 | + |
| 147 | + # standardize the eigenvector signs |
| 148 | + vec = vec[:, end:-1:1] |
| 149 | + standardize_eigenvector_signs!(vec) |
| 150 | + vec_sym = vec_sym[:, end:-1:1] |
| 151 | + standardize_eigenvector_signs!(vec_sym) |
| 152 | + |
| 153 | + # construct unfolding operator custom to current region |
| 154 | + P = Uf[indrs, :]' |
| 155 | + # obtain the expansion coefficients |
| 156 | + dmatrixlH[indr, j, :] = (P * vec)' * G.f |
| 157 | + dmatrixlHsym[indr, j, :] = (P * vec_sym)' * G.f |
| 158 | + end |
| 159 | + end |
| 160 | + |
| 161 | + return dmatrixlH, dmatrixlHsym |
| 162 | + |
| 163 | +end |
| 164 | + |
| 165 | + |
| 166 | + |
| 167 | +function standardize_eigenvector_signs!(vec) |
| 168 | + # standardize the eigenvector signs for HGLET (different with NGWPs) |
| 169 | + for col = 1:size(vec, 2) |
| 170 | + row = 1 |
| 171 | + standardized = false |
| 172 | + while !standardized |
| 173 | + if vec[row, col] > 10^3 * eps() |
| 174 | + standardized = true |
| 175 | + elseif vec[row,col] < -10^3 * eps() |
| 176 | + vec[:, col] = -vec[:, col] |
| 177 | + else |
| 178 | + row += 1 |
| 179 | + end |
| 180 | + end |
| 181 | + end |
| 182 | +end |
| 183 | + |
| 184 | +""" |
| 185 | + HGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L) |
| 186 | +
|
| 187 | +assemble the whole HGLET dictionary |
| 188 | +
|
| 189 | +### Input Arguments |
| 190 | +* `GP`: a GraphPart object |
| 191 | +* `G`: a GraphSig object |
| 192 | +* `method`: `:L` or `:Lsym` |
| 193 | +
|
| 194 | +### Output Argument |
| 195 | +* `dictionary`: the HGLET dictionary |
| 196 | +
|
| 197 | +""" |
| 198 | +function HGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L) |
| 199 | + N = size(G.W, 1) |
| 200 | + jmax = size(GP.rs, 2) |
| 201 | + dictionary = zeros(N, jmax, N) |
| 202 | + for j = 1:jmax |
| 203 | + BS = BasisSpec(collect(enumerate(j * ones(Int, N)))) |
| 204 | + dictionary[:, j, :] = HGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; method = method)[1]' |
| 205 | + end |
| 206 | + return dictionary |
| 207 | +end |
| 208 | + |
| 209 | +""" |
| 210 | + LPHGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L, ϵ::Float64 = 0.3) |
| 211 | +
|
| 212 | +assemble the whole LP-HGLET dictionary |
| 213 | +
|
| 214 | +### Input Arguments |
| 215 | +* `GP`: a GraphPart object |
| 216 | +* `G`: a GraphSig object |
| 217 | +* `method`: `:L` or `:Lsym` |
| 218 | +* `ϵ`: relative action bandwidth (default: 0.3) |
| 219 | +
|
| 220 | +### Output Argument |
| 221 | +* `dictionary`: the LP-HGLET dictionary |
| 222 | +
|
| 223 | +""" |
| 224 | +function LPHGLET_dictionary(GP::GraphPart, G::GraphSig; method::Symbol = :L, ϵ::Float64 = 0.3) |
| 225 | + N = size(G.W, 1) |
| 226 | + jmax = size(GP.rs, 2) |
| 227 | + dictionary = zeros(N, jmax, N) |
| 228 | + for j = 1:jmax |
| 229 | + BS = BasisSpec(collect(enumerate(j * ones(Int, N)))) |
| 230 | + dictionary[:, j, :] = LPHGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; method = method, ϵ = ϵ)[1]' |
| 231 | + end |
| 232 | + return dictionary |
| 233 | +end |
0 commit comments