Added the figure scripts for the eGHWT journal paper

Author: BoundaryValueProblems

test/paperscripts/eGHWT2021/5block.png

@@ -0,0 +1,130 @@
+################################################################################
+####################### Approximation error plot ###############################
+################################################################################
+function approx_error2(DVEC::Array{Array{Float64,1},1}, 
+                       T::Vector{String}, L::Vector{Tuple{Symbol,Symbol}},
+                       frac::Float64 = 0.30)
+    # This version plots the relative L2 errors against
+    # the FRACTION of coefficients retained.
+    plot(xaxis = "Fraction of Coefficients Retained", 
+         yaxis = "Relative Approximation Error")
+    for i = 1:length(DVEC)
+        dvec = DVEC[i]
+        N = length(dvec)
+        dvec_norm = norm(dvec,2)
+        dvec_sort = sort(dvec.^2) # the smallest first
+        er = sqrt.(reverse(cumsum(dvec_sort)))/dvec_norm
+        er[er .== 0.0] .= minimum(er[er .!= 0.0]) # avoid blowup by taking log in the plot below
+        # er is the relative L^2 error of the whole thing: length(er)=N.
+        p = Int64(floor(frac*N)) + 1 # upper limit
+        plot!(frac*(0:(p-1))/(p-1), er[1:p], yaxis=:log, xlims = (0.,frac), 
+              label = T[i], line = L[i], linewidth = 2, grid = false)
+    end
+    display(current())
+end
+
+
+function approx_error3(DVEC::Array{Array{Float64,1},1},
+                       T::Vector{String}, L::Vector{Tuple{Symbol,Symbol}},
+                       N::Int64)
+    # This version plots the relative L2 errors against
+    # the NUMBER of coefficients retained.
+    plot(xaxis = "Number of Coefficients Retained", 
+         yaxis = "Relative Approximation Error")
+    for i = 1:length(DVEC)
+        dvec = DVEC[i]
+        dvec_norm = norm(dvec,2)
+        dvec_sort = sort(dvec.^2) # the smallest first
+        er = sqrt.(reverse(cumsum(dvec_sort)))/dvec_norm 
+        er[er .== 0.0] .= minimum(er[er .!= 0.0]) # avoid blowup by taking log in the plot below
+        # er is the relative L^2 error of the whole thing: length(er)=N.
+        plot!(0:N-1, er[1:N], yaxis=:log, label = T[i], line = L[i], 
+              linewidth = 2, grid = false)
+    end
+    display(current())
+end
+
+################################################################################
+############ Computing a weight matrix using the Gaussian affinity #############
+################################################################################
+function image_Gaussian_affinity(img::Matrix{Float64}, r::Int64, σ::Float64)
+# Get the neighbors for affinity matrix computation
+    l = 2*r + 1
+    temp_x, temp_y = fill(1,l^2), fill(1,l^2)
+    temp_xy = CartesianIndices((1:l,1:l))
+    for i = 1:l^2
+        temp_x[i], temp_y[i] = temp_xy[i][1], temp_xy[i][2]
+    end
+    # (r+1, r+1) is the index of the center location.
+    temp_ind = ((temp_x .- (r + 1)).^2 + (temp_y .- (r + 1)).^2) .<= r^2
+    # Now, temp_ind indicates those points withinin the circle of radius r 
+    # from the center while neighbor_x, neighbor_y represent relative positions
+    # of those points from the center.
+    neighbor_x = temp_x[temp_ind] .- (r + 1) 
+    neighbor_y = temp_y[temp_ind] .- (r + 1)
+   # So, for any index (x, y), points within (x ± neighbor_x, y ± neighbor_y) are 
+   # its neighbors for the purpose of calculating affinity
+
+# Create affinity matrix
+    m, n = size(img)
+    sig = img[:]
+    W = fill(0., (m*n,m*n))
+    for i = 1:m*n
+        cur = CartesianIndices((m,n))[i]
+        for j = 1:length(neighbor_x) # assuming dim(neighbor_x) == dim(neighbor_y) 
+            if 1 <= cur[1] + neighbor_x[j] <= m && 1 <= cur[2] + neighbor_y[j] <= n
+                tempd = LinearIndices((m,n))[cur[1] + neighbor_x[j], cur[2] + neighbor_y[j]]
+                W[i,tempd] = exp(-(sig[i] - sig[tempd])^2/σ)
+            end
+        end
+    end
+    return sparse(W)
+end # end of the image_Gaussian_affinity function
+
+################################################################################
+######## Display top 9 basis vectors of various bases for an image data ########
+################################################################################
+function top_vectors_plot2(dvec::Array{Float64, 2}, m::Int64, n::Int64, BS::BasisSpec, GP::GraphPart;
+                           clims::Tuple{Float64,Float64} = (-0.01, 0.01))
+
+    # Get the indices from max to min in terms of absolute value of the coefs
+    # Note that m*n == length(dvec); m and n are the image size.
+    sorted_ind = sortperm(abs.(dvec[:]), rev = true);
+
+    # Set up the layout as 3 x 3 subplots
+    plot(9, layout = (3,3), framestyle = :none, legend = false)
+
+    # Do the display!
+    for i=1:9
+        dvecT = fill(0., size(dvec))
+        dvecT[sorted_ind[i]] = 1
+        f = ghwt_synthesis(dvecT, GP, BS)
+        heatmap!(reshape(f, (m,n)), subplot=i, ratio=1, yaxis=:flip, 
+                 showaxis=false, ticks = false, color = :grays, clims = clims)
+    end
+    display(current())
+end
+################################################################################
+
+function top_vectors_plot3(dvec::Array{Float64, 2}, m::Int64, n::Int64, BS::BasisSpec,
+    GP::GraphPart; clims::Tuple{Float64,Float64} = (-0.01, 0.01), K::Int64 = 9)
+
+    # Get the indices from max to min in terms of absolute value of the coefs
+    # Note that m*n == length(dvec); m and n are the image size.
+    sorted_ind = sortperm(abs.(dvec[:]), rev = true);
+
+    # Set up the layout as 3 x 3 subplots
+    plot(K, layout = (Int(sqrt(K)), Int(sqrt(K))), framestyle = :none, legend = false)
+
+    # Do the display!
+    for i=1:K
+        dvecT = fill(0., size(dvec))
+        dvecT[sorted_ind[i]] = 1
+        f = ghwt_synthesis(dvecT, GP, BS)
+        heatmap!(reshape(f, (m,n)), subplot=i, ratio=1, yaxis=:flip, 
+                 showaxis=false, ticks = false, color = :grays, clims = clims)
+    end
+    display(current())
+end
+################################################################################
+

test/paperscripts/eGHWT2021/Toronto_new.jld2

@@ -0,0 +1,60 @@
+# Script to generate figures in Figure 2.
+using Plots, SparseArrays, JLD2, LinearAlgebra, MultiscaleGraphSignalTransforms
+
+# First define a utility function "Glevel", which generate a Graph Signal that 
+# are piecewise constants whose values are actual region indices "k".
+function Glevel(G::GraphSig, GP::GraphPart, j::Int64)
+    f = zeros(size(G.f))
+    for k in 1:size(GP.rs,1)
+        a = GP.rs[k,j+1]
+        b = GP.rs[k+1,j+1] - 1
+        if b == -1
+            break
+        end
+        f[a:b] .= k*1.0
+    end
+    Gsub = deepcopy(G)
+    Gsub.f[GP.ind] = f
+    return Gsub
+end
+
+# Set up the resolution and display size
+gr(dpi=200, size=(800,600))
+
+# Load the Toronto street network (vehicular volume counts as its graph signal)
+JLD2.@load "Toronto_new.jld2"
+tmp1 = toronto["G"]
+G=GraphSig(tmp1["W"],xy=tmp1["xy"],f=tmp1["f"],name =tmp1["name"],
+    plotspecs = tmp1["plotspecs"])
+
+# Assign correct edge weights via 1/Euclidean distance
+G = Adj2InvEuc(G)
+# Perform the hierarchical graph partitioning using the Fiedler vectors
+GP = partition_tree_fiedler(G) # Lrw by default
+# Compute the expansion coefficients of the full GHWT c2f dictionary
+dmatrix = ghwt_analysis!(G, GP=GP)
+N = length(G.f)
+
+# Fig. 2a: Level 1 first since Level 0 is not interesting
+j = 1
+GraphSig_Plot(Glevel(G,GP,j), linewidth = 1., markersize = 4., 
+    markercolor = :viridis, markerstrokealpha =0., notitle = true, nocolorbar = true)
+plot!(showaxis = false, ticks = false)
+savefig("toronto_j1.pdf")
+savefig("toronto_j1.png")
+
+# Fig. 2b: Level 2
+j = 2
+GraphSig_Plot(Glevel(G,GP,j), linewidth = 1., markersize = 4., 
+    markercolor = :viridis, markerstrokealpha =0., notitle = true, nocolorbar = true)
+plot!(showaxis = false, ticks = false)
+savefig("toronto_j2.pdf")
+savefig("toronto_j2.png")
+
+# Fig. 2c: Level 3
+j = 3
+GraphSig_Plot(Glevel(G,GP,j), linewidth = 1., markersize = 4., 
+    markercolor = :viridis, markerstrokealpha =0., notitle = true, nocolorbar = true)
+plot!(showaxis = false, ticks = false)
+savefig("toronto_j3.pdf")
+savefig("toronto_j3.png")

test/paperscripts/eGHWT2021/auxilaries.jl

@@ -0,0 +1,67 @@
+# Script to generate figures in Figure 4.
+using Plots, SparseArrays, JLD2, LinearAlgebra, MultiscaleGraphSignalTransforms
+
+# Set up the resolution and display size
+gr(dpi=200, size=(800,600))
+
+# Load the Toronto street network (vehicular volume counts as its graph signal)
+JLD2.@load "Toronto_new.jld2"
+tmp1 = toronto["G"]
+G=GraphSig(tmp1["W"],xy=tmp1["xy"],f=tmp1["f"],name =tmp1["name"],
+    plotspecs = tmp1["plotspecs"])
+
+# Assign correct edge weights via 1/Euclidean distance
+G = Adj2InvEuc(G)
+# Perform the hierarchical graph partitioning using the Fiedler vectors
+GP = partition_tree_fiedler(G) # Lrw by default
+# Compute the expansion coefficients of the full GHWT c2f dictionary
+dmatrix = ghwt_analysis!(G, GP=GP)
+N = length(G.f)
+
+# Fig. 4a: Sacling vector ψ^1_{0,0}, i.e., (j,k,l)=(1,0,0).
+j = 1
+loc = 1
+BS = bs_level(GP, j)
+dvec = zeros(N, 1)
+dvec[loc,1] = 1.0
+(f, GS) = ghwt_synthesis(dvec, GP, BS, G)
+GraphSig_Plot(GS, linewidth = 1., markersize = 4., markerstrokealpha = 0., 
+    markercolor = :viridis, notitle = true, nocolorbar = true, clim = (-0.01, 0.01))
+plot!(showaxis = false, ticks = false)
+k, l = BS.levlist[loc][1], BS.levlist[loc][2]
+print(j," ",(rs_to_region(GP.rs, GP.tag))[k,l]," ",GP.tag[k,l])
+savefig("toronto_psi100.pdf")
+savefig("toronto_psi100.png")
+
+# Fig. 4b: Haar vector ψ^2_{0,1}, i.e., (j,k,l)=(2,0,1).
+j = 2
+loc = 2
+BS = bs_level(GP, j)
+dvec = zeros(N, 1)
+dvec[loc,1] = 1.0
+(f, GS) = ghwt_synthesis(dvec, GP, BS, G)
+GraphSig_Plot(GS, linewidth = 1., markersize = 4., markerstrokealpha = 0., 
+    markercolor = :viridis, notitle = true, nocolorbar = true, clim = (-0.01, 0.01))
+plot!(showaxis = false, ticks = false)
+k, l = BS.levlist[loc][1], BS.levlist[loc][2]
+print(j," ",(rs_to_region(GP.rs, GP.tag))[k,l]," ",GP.tag[k,l])
+savefig("toronto_psi201.pdf")
+savefig("toronto_psi201.png")
+
+# Fig. 4c: Walsh vector ψ^3_{0,9}, i.e., (j,k,l)=(3,0,9).
+#j = 5
+#loc = 1000
+j = 3
+loc = 10
+BS = bs_level(GP, j)
+dvec = zeros(N, 1)
+dvec[loc,1] = 1
+(f, GS) = ghwt_synthesis(dvec, GP, BS, G)
+GraphSig_Plot(GS, linewidth = 1., markersize = 4., markerstrokealpha = 0., 
+    markercolor = :viridis, notitle = true, nocolorbar = true, clim = (-0.01, 0.01))
+plot!(showaxis = false, ticks = false)
+k, l = BS.levlist[loc][1], BS.levlist[loc][2]
+print(j," ",(rs_to_region(GP.rs, GP.tag))[k,l]," ",GP.tag[k,l])
+savefig("toronto_psi309.pdf")
+savefig("toronto_psi309.png")
+

test/paperscripts/eGHWT2021/barbara.jld2

@@ -0,0 +1,75 @@
+# Script to generate Figure 15.
+using Plots, SparseArrays, JLD2, LinearAlgebra, Wavelets, MultiscaleGraphSignalTransforms
+include("auxilaries.jl")
+include("../../../src/unbalanced_haar_image.jl")
+
+# Set up the resolution and display size
+gr(dpi=200, size=(800,600))
+
+# Load the Barbara image data
+JLD2.@load "barbara.jld2"
+
+# smaller face region of size 100x100
+row_zoom = 71:170
+col_zoom = 341:440
+display(heatmap(barbara[row_zoom, col_zoom],ratio=1, yaxis =:flip, 
+    showaxis = false, ticks = false, color = :grays, clims = (0,1)))
+
+# Extract that portion of the image
+matrix = deepcopy(barbara[row_zoom, col_zoom])
+
+#
+# Perform the image partition using the penalized total variation and then
+# compute the expansion coefficients
+#
+p = 3.0
+maxlev = 9 # need this to reach the single node leaves
+GPcols = PartitionTreeMatrix_unbalanced_haar_binarytree(matrix, maxlev, p)
+maxlev = 8 # need this to reach the single node leaves
+GProws = PartitionTreeMatrix_unbalanced_haar_binarytree(copy(matrix'), maxlev, p)
+# copy() is necessary to switch the Adjoint type to a regular Matrix{Float64}.
+dmatrix = ghwt_analysis_2d(matrix, GProws, GPcols)
+
+#
+# Compute the graph Haar coefficients from the previous PTV partition tree
+#
+BS_haar_rows = bs_haar(GProws)
+BS_haar_cols = bs_haar(GPcols)
+dvec_haar, loc_haar = BS2loc(dmatrix, GProws, GPcols, BS_haar_rows, BS_haar_cols)
+
+#
+# Compute the eGHWT best basis
+#
+dvec_eghwt, loc_eghwt = eghwt_bestbasis_2d(matrix, GProws, GPcols)
+
+#
+# Classical Haar transform via Wavelets.jl with direct input
+#
+dvec_classichaar = dwt(matrix, wavelet(WT.haar))
+
+#
+# Classical Haar transform via Wavelets.jl via putting in the dyadic image
+#
+matrix2 = zeros(128,128)
+matrix2[15:114,15:114] = deepcopy(matrix)
+dvec_classichaar2 = dwt(matrix2, wavelet(WT.haar))
+
+#
+# Classical Haar transform via Wavelets.jl via putting in the dyadic image
+#
+matrix3 = zeros(128,128)
+matrix3[1:100,1:100] = deepcopy(matrix)
+matrix3[101:end,1:100] = deepcopy(matrix[end:-1:end-27,1:100])
+matrix3[:,101:end] = matrix3[:,100:-1:73]
+dvec_classichaar3 = dwt(matrix3, wavelet(WT.haar))
+
+# Figure 15 (Approximation error plot up to the top 5000 coefficients)
+DVEC = [ dvec_classichaar[:], dvec_classichaar2[:], dvec_classichaar3[:], 
+    dvec_haar[:], dvec_eghwt[:] ]
+T = [ "Classic Haar (direct)", "Classic Haar (zero pad)", 
+    "Classic Haar (even refl)", "Graph Haar (PTV cost)", "eGHWT (PTV cost)" ]
+L = [ (:dashdot, :orange), (:dashdot, :red), (:dashdot, :green),
+      (:solid, :blue), (:solid, :black) ]
+approx_error3(DVEC, T, L, 5000)
+savefig("barbara_face100_nondyadic_haar_error.pdf")
+savefig("barbara_face100_nondyadic_haar_error.png")