add figure scripts of NGWP JFAA paper and rearrange paperscripts folder

Author: haotian127

test/paperscripts/NGWP_JFAA2021/datasets/barbara_gray_matrix.jld

@@ -0,0 +1,110 @@
+# script for Fig.1, Fig.2, Fig.4, Fig.6
+
+using MultiscaleGraphSignalTransforms, LightGraphs, Plots, LaTeXStrings, MultivariateStats
+pyplot(dpi = 200)
+
+Nx, Ny = 7, 3
+G = LightGraphs.grid([Nx, Ny]); N = nv(G);
+L = Matrix(laplacian_matrix(G))
+Q = incidence_matrix(G; oriented = true)
+𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1, :])';
+
+#################### Fig. 1(a) eigenvectors by nondecreasing eigenvalue ordering
+plot(layout = Plots.grid(3, 7))
+for i in 1:N
+    heatmap!(reshape(𝚽[:, i],(Nx, Ny))', c = :viridis, cbar = false,
+                clims = (-0.4,0.4), frame = :none, ratio = 1,
+                title = latexstring("\\phi_{", i-1, "}"), titlefont = 12,
+                subplot = i)
+end
+plt = current()
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/grid7x3_evsp_title.png"))
+
+#################### Fig. 1(b) eigenvectors by natural frequency ordering
+# find correct 2D index
+grid2eig_ind = [1,2,3,6,8,12,15,4,5,7,9,13,16,18,10,11,14,17,19,20,21];
+eig2grid_ind = sortperm(grid2eig_ind);
+eig2dct = Array{Int64,3}(undef, Nx, Ny, 2);
+for i = 1:Nx; for j = 1:Ny; eig2dct[i,j,1] = i-1; eig2dct[i,j,2] = j-1; end; end
+eig2dct = reshape(eig2dct, (N, 2)); eig2dct = eig2dct[eig2grid_ind, :];
+
+plot(layout = Plots.grid(3, 7))
+for i in 1:N
+    k = grid2eig_ind[i]
+    heatmap!(reshape(𝚽[:,k],(Nx,Ny))', c = :viridis, cbar = false,
+                clims = (-0.4,0.4), frame = :none, ratio = 1,
+                title = latexstring("\\varphi_{", string(eig2dct[k,1]),
+                ",", string(eig2dct[k,2]), "}"), titlefont = 12, subplot = i)
+end
+plt = current()
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/grid7x3_dct_title2.png"))
+
+# DAG pseudo-metric
+distDAG = eigDAG_Distance(𝚽, Q, N)
+
+# MDS embedding into R²
+D = distDAG
+E = transform(fit(MDS, D, maxoutdim=2, distances=true))
+
+# set up all heatmap plots' positions
+dx = 0.01; dy = dx;
+xej = zeros(Nx, N); yej=zeros(Ny, N);
+a = 5.0; b = 9.0;
+for k = 1:N
+    xej[:,k] = LinRange(E[1,k] - Ny * a * dx, E[1, k] + Ny * a * dx, Nx)
+    yej[:,k] = LinRange(E[2,k] - a * dy, E[2, k] + a * dy, Ny)
+end
+
+#################### Fig. 2
+plot()
+for k=1:N
+    heatmap!(xej[:, k], yej[:, k], reshape(𝚽[:, k], (Nx, Ny))', c = :viridis,
+                colorbar = false, ratio = 1, annotations = (xej[4, k],
+                yej[3, k] + b*dy, text(latexstring("\\varphi_{",
+                string(eig2dct[k, 1]), ",", string(eig2dct[k, 2]), "}"), 10)))
+end
+plt = plot!(xlim = [-1.4, 1.3], ylim = [-1.4, 1.3], grid = false, clims = (-0.4, 0.4))
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/Grid7x3_DAG_MDS.png"))
+
+#################### Fig. 4
+# first level partition
+p1x = [-0.2, 1.0, NaN]; p1y = [1.3, -1.0, NaN];
+plot!(p1x, p1y, c = :red, legend = false, width = 3)
+# second level partition
+p2x = [-1.0, 0.2, NaN, 0.4, 1.2, NaN]; p2y = [-0.8, 0.45, NaN, 0.25, 0.2, NaN];
+plot!(p2x, p2y, c=:orange, legend = false, width = 2)
+plt = current()
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/Grid7x3_DAG_2levels_partition.png"))
+
+
+## Build Dual Graph
+Gstar_Sig = dualgraph(distDAG)
+GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
+GP_primal = pairclustering(𝚽, GP_dual)
+
+@time VM_NGWP = vm_ngwp(𝚽, GP_dual)
+
+## level 2 VM-NGWP vectors
+j = 3; W_VM = VM_NGWP[:, j, :]'
+
+wav_kl = [[0 0];[0 1];[0 2];[1 0];[1 1];[1 2];[1 3];[2 0];[2 1];[2 2];[3 0];
+            [3 1];[3 2];[3 3];[2 3];[2 4];[2 5];[2 6];[3 4];[3 5];[3 6]];
+wav_kl = wav_kl[eig2grid_ind,:];
+
+# reorder_ind = [2,3,1,5,7,4,6, 16,17,15,  9,11,8,10, 18,20,21,19, 13,14,12]
+reorder_ind = [1,3,2,5,7,4,6, 9,10,8,16,18,15,17, 11,13,14,12,20,21,19]
+W_VM = W_VM[:,reorder_ind[eig2grid_ind]];
+sgn = ones(N); sgn[grid2eig_ind[[4,6,8,10,14]]] .= -1; W_VM = W_VM * Diagonal(sgn);
+
+#################### Fig. 6
+plot()
+for k=1:N
+    heatmap!(xej[:,k],yej[:,k],reshape(W_VM[:,k],(Nx,Ny))',c=:viridis,colorbar=false,ratio=1,annotations=(xej[4,k], yej[3,k]+b*dy, text(latexstring("\\psi_{", string(wav_kl[k,1]), ",", string(wav_kl[k,2]), "}"),10)))
+end
+plot!(aspect_ratio = 1, xlim = [-1.4, 1.3], ylim = [-1.4, 1.3], grid = false, clims=(-0.34,0.34))
+# first level partition
+p1x = [-0.2, 1.0, NaN]; p1y = [1.3, -1.0, NaN]; plot!(p1x, p1y, c = :red, legend = false, width = 3)
+# second level partition
+p2x = [-1.0, 0.2, NaN, 0.4, 1.2, NaN]; p2y = [-0.8, 0.45, NaN, 0.25, 0.2, NaN]; plot!(p2x, p2y, c=:orange, legend = false, width = 2)
+plt = current()
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/Grid7x3_DAG_VM_NGWP_lvl2_wavelets.png"))

test/paperscripts/NGWP_JFAA2021/datasets/new_toronto.jld

@@ -0,0 +1,31 @@
+# script for Fig.5
+
+using MultiscaleGraphSignalTransforms, LightGraphs, Plots; gr(dpi = 200)
+import WaveletsExt: wiggle
+
+## Build Graph
+N = 512; G = path_graph(N)
+X = zeros(N,2); X[:, 1] = 1:N
+L = Matrix(laplacian_matrix(G))
+𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1,:])'
+W = 1.0 * adjacency_matrix(G)
+
+## Build NGWPs
+Gstar_Sig = GraphSig(W)
+G_Sig = GraphSig(W, xy = X)
+GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
+GP_primal = pairclustering(𝚽, GP_dual)
+
+@time VM_NGWP = vm_ngwp(𝚽, GP_dual)
+
+#################### Fig.5
+j = 5
+for k in [1, 2, 5]
+   WW = sort_wavelets(VM_NGWP[GP_dual.rs[k, j]:(GP_dual.rs[k + 1, j] - 1), j, :]')
+   if k == 2
+      WW[:, end] *= -1
+   end
+   plt = wiggle(WW; sc = 0.75)
+   # savefig(plt, joinpath(@__DIR__, "../paperfigs/Path512_VM_NGWP_j$(j-1)k$(k-1).png"))
+end
+current()

test/paperscripts/NGWP_JFAA2021/datasets/new_toronto_graph.lgz

@@ -0,0 +1,82 @@
+# script for Fig.8(b)(c), Fig.9, Fig.10(b)(c), Fig.11
+
+using MultiscaleGraphSignalTransforms, LightGraphs, Plots; gr(dpi = 200)
+
+## Build weighted graph
+G, L, X = SunFlowerGraph(N = 400); N = nv(G)
+𝛌, 𝚽 = eigen(Matrix(L))
+sgn = (maximum(𝚽, dims = 1)[:] .> -minimum(𝚽, dims = 1)[:]) .* 2 .- 1
+𝚽 = 𝚽 * Diagonal(sgn)
+Q = incidence_matrix(G; oriented = true)
+W = 1.0 * adjacency_matrix(G)
+edge_weight = [e.weight for e in edges(G)]
+
+## Build Dual Graph by DAG metric
+distDAG = eigDAG_Distance(𝚽, Q, N; edge_weight = edge_weight)
+Gstar_Sig = dualgraph(distDAG)
+G_Sig = GraphSig(W, xy = X)
+GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
+GP_primal = pairclustering(𝚽, GP_dual)
+
+@time VM_NGWP = vm_ngwp(𝚽, GP_dual) #54.986524 seconds (1.19 M allocations: 25.127 GiB, 2.36% gc time)
+@time PC_NGWP = pc_ngwp(𝚽, GP_dual, GP_primal) #0.611176 seconds (225.41 k allocations: 844.488 MiB, 11.71% gc time)
+
+
+#################### Fig. 8(b) barbara eye graph signal
+using MAT
+f = matread(joinpath(@__DIR__, "../datasets",
+            "sunflower_barbara_voronoi.mat"))["f_eye_voronoi"]
+G_Sig.f = reshape(f, (N, 1))
+scatter_gplot(X; marker = f, ms = LinRange(4.0, 14.0, N), c = :greys);
+plt = plot!(xlim = [-1.2,1.2], ylim = [-1.2,1.2], frame = :none)
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_feye.png"))
+
+#################### Fig. 8(c) barbara eye relative l2 approximation error by various methods
+DVEC = getall_expansioncoeffs(G_Sig, GP_dual, VM_NGWP, PC_NGWP, 𝚽)
+approx_error_plot(DVEC);
+plt = plot!(xguidefontsize = 16, yguidefontsize = 16, legendfontsize = 12)
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_feye_DAG_approx.png"))
+
+#################### Fig. 9 barbara eye 16 most important VM-NGWP vectors (ignore the DC vector)
+dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
+dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
+important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
+for i in 2:17
+    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
+    w = VM_NGWP[dr, dc, :]
+    println("(j, k, l) = ", NGWP_jkl(GP_dual, dr, dc))
+    scatter_gplot(X; marker = w, ms = LinRange(4.0, 14.0, N), c = :greys)
+    plt = plot!(xlim = [-1.2, 1.2], ylim = [-1.2, 1.2], frame = :none,
+                cbar = false, clims = (-0.15, 0.15))
+    # savefig(plt, joinpath(@__DIR__,
+    #     "../paperfigs/SunFlower_barbara_feye_DAG_VM_NGW_important_basis_vector$(lpad(i,2,"0")).png"))
+end
+
+#################### Fig. 10(b) barbara pants graph signal
+f = matread(joinpath(@__DIR__, "../datasets",
+            "sunflower_barbara_voronoi.mat"))["f_trouser_voronoi"]
+scatter_gplot(X; marker = f, ms = LinRange(4.0, 14.0, N), c = :greys);
+plt = plot!(xlim = [-1.2,1.2], ylim = [-1.2,1.2], frame = :none)
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_ftrouser.png"))
+
+#################### Fig. 10(c) barbara eye relative l2 approximation error by various methods
+G_Sig.f = reshape(f, (N, 1))
+DVEC = getall_expansioncoeffs(G_Sig, GP_dual, VM_NGWP, PC_NGWP, 𝚽)
+approx_error_plot(DVEC);
+plt = plot!(xguidefontsize = 16, yguidefontsize = 16, legendfontsize = 12)
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower_barbara_ftrouser_DAG_approx.png"))
+
+#################### Fig. 11 barbara pants 16 most important VM-NGWP vectors (ignore the DC vector)
+dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
+dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
+important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
+for i in 2:17
+    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
+    w = VM_NGWP[dr, dc, :]
+    println("(j, k, l) = ", NGWP_jkl(GP_dual, dr, dc))
+    scatter_gplot(X; marker = w, ms = LinRange(4.0, 14.0, N), c = :greys)
+    plt = plot!(xlim = [-1.2, 1.2], ylim = [-1.2, 1.2], frame = :none,
+                cbar = false, clims = (-0.15, 0.15))
+    # savefig(plt, joinpath(@__DIR__,
+    #     "../paperfigs/SunFlower_barbara_ftrouser_DAG_VM_NGW_important_basis_vector$(lpad(i,2,"0")).png"))
+end

test/paperscripts/NGWP_JFAA2021/datasets/sunflower_barbara.mat

@@ -0,0 +1,34 @@
+## script for Fig.7, Fig.8(a), Fig.10(a)
+using VoronoiDelaunay, VoronoiCells, GeometricalPredicates
+using MultiscaleGraphSignalTransforms, Plots, LightGraphs, JLD; gr(dpi = 200)
+
+barbara = JLD.load(joinpath(@__DIR__, "..", "datasets", "barbara_gray_matrix.jld"), "barbara")
+G, L, X = SunFlowerGraph(); N = nv(G)
+#################### Fig. 7(a) sunflower graph
+gplot(1.0*adjacency_matrix(G),X; width=1); scatter_gplot!(X; c = :red, ms = LinRange(1,9,N)); plt = plot!(frame = :none)
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/SunFlower.png"))
+
+## Voronoi tessellation
+width_x = maximum(abs.(X[:,1])) * 2; width_y = maximum(abs.(X[:,2])) * 2;
+width = VoronoiDelaunay.max_coord - VoronoiDelaunay.min_coord
+center_coord = (VoronoiDelaunay.min_coord + VoronoiDelaunay.max_coord)/2
+X_transform = zeros(N,2)
+for i in 1:N
+    X_transform[i,:] = X[i,:] ./ [width_x/width, width_y/width] + [center_coord, center_coord]
+end
+
+pts = [Point2D(X_transform[i,1], X_transform[i,2]) for i in 1:N]
+tess = DelaunayTessellation(N)
+push!(tess, pts)
+#################### Fig. 7(b) voronoi tessellation
+xx, yy = getplotxy(voronoiedges(tess))
+plt = plot(xx, yy, xlim=[1,2], ylim=[1,2], linestyle=:auto, linewidth=1, linecolor=:blue, grid=false, label="", aspect_ratio=1, frame=:box)
+# savefig(plt, joinpath(@__DIR__, "../paperfigs/Sunflower_Barbara_Voronoi_cells.png"))
+
+#################### Fig. 8(a) barbara sunflower eye
+heatmap(barbara, yflip=true, ratio=1, c=:greys); scatter_gplot!(transform2D(X; s = 20, t = [395, 100]); ms = 2, c = :red); sample_location_plt = plot!(cbar = false, frame = :none)
+# savefig(sample_location_plt, joinpath(@__DIR__, "../paperfigs/barb_sunflower_eye.png"))
+
+#################### Fig. 10(a) barbara sunflower pants
+heatmap(barbara, yflip=true, ratio=1, c=:greys); scatter_gplot!(transform2D(X; s = 20, t = [280, 320]); ms = 2, c = :red); sample_location_plt = plot!(cbar = false, frame = :none)
+# savefig(sample_location_plt, joinpath(@__DIR__, "../paperfigs/barb_sunflower_pants.png"))