Merge pull request #7 from haotian127/master

update HL's dissertation figures and a couple of small fixes

Author: BoundaryValueProblems

Project.toml

@@ -1,7 +1,7 @@
 name = "MultiscaleGraphSignalTransforms"
 uuid = "cc44fc10-d0fb-5545-af72-ba1cb661a341"
 authors = ["Naoki Saito <saito@math.ucdavis.edu>", "Yiqun Shao <shaoyiqun1992@gmail.com>", "Haotian Li <htsli@ucdavis.edu>"]
-version = "1.5.1"
+version = "1.5.2"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"

src/LP-NGWP.jl

@@ -247,3 +247,49 @@ function const_meyer_wavelets(𝚽, Uf; idx = 1:size(Uf, 1))
     Wavelets = varimax(B)
     return Wavelets
 end
+
+"""
+    lp_ngwp_analysis(G::GraphSig, 𝚽::Matrix{Float64}, W_dual::SparseMatrixCSC{Float64, Int64},
+            GP_dual::GraphPart; ϵ::Float64 = 0.3)
+
+perform the LP-NGWP transform of the graph signal(s) `G.f`.
+
+# Input Arguments
+- `G::GraphSig`: a GraphSig object
+- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors
+- `W_dual::SparseMatrixCSC{Float64, Int64}`: weight matrix of the dual graph
+- `GP_dual::GraphPart`: GraphPart object of the dual graph
+
+# Output Argument
+- `wavelet_packet::Array{Float64,3}`: the lapped NGWP dictionary. The first
+    index is for selecting wavelets at a fixed level; the second index is for
+    selecting the level `j`; the third index is for selecting elements in the
+    wavelet vector.
+
+"""
+function lp_ngwp_analysis(G::GraphSig, 𝚽::Matrix{Float64}, W_dual::SparseMatrixCSC{Float64, Int64},
+                 GP_dual::GraphPart; ϵ::Float64 = 0.3)
+    rs = GP_dual.rs
+    inds = GP_dual.inds
+    (N, jmax) = Base.size(inds)
+    fcols = size(G.f, 2)
+
+    GP_dual.tag = zeros(Int, N, jmax)
+    GP_dual.tag[:, 1] = Vector{Int}(0:(N - 1))
+
+    Uf = Matrix{Float64}(I, N, N)
+    used_node = Set()
+
+    dmatrix = zeros(N, jmax, fcols)
+    dmatrix[:, 1, :] = G.f
+    for j = 2:jmax
+        regioncount = count(!iszero, rs[:, j]) - 1
+        keep_folding!(Uf, used_node, W_dual, GP_dual; ϵ = ϵ, j = j - 1)
+        for r = 1:regioncount
+            indr = rs[r, j]:(rs[r + 1, j] - 1)
+            GP_dual.tag[indr, j] = Vector{Int}(0:(length(indr) - 1))
+            dmatrix[indr, j, :] = const_meyer_wavelets(𝚽, Uf; idx = inds[indr, j])' * G.f
+        end
+    end
+    return dmatrix
+end

src/MultiscaleGraphSignalTransforms.jl

@@ -52,7 +52,7 @@ export natural_eigdist
 export SunFlowerGraph, dualgraph
 export pc_ngwp, pairclustering
 export vm_ngwp, varimax
-export lp_ngwp, rising_cutoff, find_pairinds, pair_inds_shadding
+export lp_ngwp, rising_cutoff, find_pairinds, pair_inds_shadding, lp_ngwp_analysis
 export LPHGLET_Synthesis, LPHGLET_Analysis_All, HGLET_dictionary, LPHGLET_dictionary
 export unitary_folding_operator, keep_folding!
 export ngwp_analysis, ngwp_bestbasis, NGWP_jkl

src/helpers.jl

@@ -121,52 +121,74 @@ SCATTER\\_GPLOT!(X; ...) adds a plot to `current` one.
 
 """
 function scatter_gplot(X; marker = nothing, ms = 4, plotOrder = :normal, c = :viridis)
-    dim = size(X,2)
-    if marker != nothing && plotOrder != :normal
-        if plotOrder == :s2l
+    N, dim = size(X)
+    if marker != nothing
+        if size(marker) == (N,) || size(marker) == (N, 1)
+            marker = marker[:]  # reshape N x 1 matrix to a vector of length N
+        else
+            error("marker only accepts a vector of length $(N) or a matrix of size $(N) x 1.")
+        end
+        if plotOrder == :normal
+            idx = 1:N
+        elseif plotOrder == :s2l
             idx = sortperm(marker)
         elseif plotOrder == :l2s
-            idx = sortperm(marker, rev=true)
+            idx = sortperm(marker, rev = true)
         else
-            print("Error: plotOrder only supports for :normal, :s2l, or :l2s.")
+            error("plotOrder only supports for :normal, :s2l, or :l2s.")
         end
-        X = X[idx,:]
+        X = X[idx, :]
         marker = marker[idx]
         if length(ms) > 1
             ms = ms[idx]
         end
     end
     if dim == 2
-        scatter(X[:,1],X[:,2], marker_z = marker, ms = ms, c = c, legend = false, mswidth = 0, cbar = true, aspect_ratio = 1, grid = false)
+        scatter(X[:, 1], X[:, 2], marker_z = marker, ms = ms, c = c,
+                legend = false, mswidth = 0, cbar = true, aspect_ratio = 1,
+                grid = false)
     elseif dim == 3
-        scatter(X[:,1],X[:,2],X[:,3], marker_z = marker, ms = ms, c = c, legend = false, mswidth = 0, cbar = true, aspect_ratio = 1, grid = false)
+        scatter(X[:, 1], X[:, 2], X[:, 3], marker_z = marker, ms = ms, c = c,
+                legend = false, mswidth = 0, cbar = true, aspect_ratio = 1,
+                grid = false)
     else
-        print("Dimension Error: scatter_gplot only supports for 2-dim or 3-dim scatter plots.")
+        error("Dimension Error: scatter_gplot only supports for 2-dim or 3-dim scatter plots.")
     end
 end
 
 function scatter_gplot!(X; marker = nothing, ms = 4, plotOrder = :normal, c = :viridis)
-    dim = size(X,2)
-    if marker != nothing && plotOrder != :normal
-        if plotOrder == :s2l
+    N, dim = size(X)
+    if marker != nothing
+        if size(marker) == (N,) || size(marker) == (N, 1)
+            marker = marker[:]  # reshape N x 1 matrix to a vector of length N
+        else
+            error("marker only accepts a vector of length $(N) or a matrix of size $(N) x 1.")
+        end
+        if plotOrder == :normal
+            idx = 1:N
+        elseif plotOrder == :s2l
             idx = sortperm(marker)
         elseif plotOrder == :l2s
-            idx = sortperm(marker, rev=true)
+            idx = sortperm(marker, rev = true)
         else
-            print("Error: plotOrder only supports for :normal, :s2l, or :l2s.")
+            error("plotOrder only supports for :normal, :s2l, or :l2s.")
         end
-        X = X[idx,:]
+        X = X[idx, :]
         marker = marker[idx]
         if length(ms) > 1
             ms = ms[idx]
         end
     end
     if dim == 2
-        scatter!(X[:,1],X[:,2], marker_z = marker, ms = ms, c = c, legend = false, mswidth = 0, cbar = true, aspect_ratio = 1, grid = false)
+        scatter!(X[:, 1], X[:, 2], marker_z = marker, ms = ms, c = c,
+                legend = false, mswidth = 0, cbar = true, aspect_ratio = 1,
+                grid = false)
     elseif dim == 3
-        scatter!(X[:,1],X[:,2],X[:,3], marker_z = marker, ms = ms, c = c, legend = false, mswidth = 0, cbar = true, aspect_ratio = 1, grid = false)
+        scatter!(X[:, 1], X[:, 2], X[:, 3], marker_z = marker, ms = ms, c = c,
+                legend = false, mswidth = 0, cbar = true, aspect_ratio = 1,
+                grid = false)
     else
-        print("Dimension Error: scatter_gplot! only supports for 2-dim or 3-dim scatter plots.")
+        error("Dimension Error: scatter_gplot! only supports for 2-dim or 3-dim scatter plots.")
     end
 end
 

test/dissertations/htli/figs/Grid23x22_fdigit3_rngwf_top100.png

@@ -8,7 +8,7 @@ SGWT_dual = (SGWT * SGWT') \ SGWT
 f = digit_img[:]
 
 important_idx = sortperm((SGWT' * f).^2; rev = true)
-plot(layout = Plots.grid(10,10), size = (2300, 2200))
+plot(layout = Plots.grid(10,10), size = (2300, 2800))
 for i = 1:100
     grid_vector_plot!(important_idx[i], i, SGWT)
     if i == 90

test/dissertations/htli/figs/Grid23x22_fdigit3_sgwt_top100.png

@@ -9,12 +9,15 @@ f = digit_img[:]
 
 important_idx = sortperm((rNGWF' * f).^2; rev = true)
 red_box_inds = [41, 43, 66, 69, 71, 73, 75, 84, 85, 86, 89, 99]
-plot(layout = Plots.grid(10,10), size = (2300, 2200))
+orange_box_inds = [40, 53, 54, 55, 72]
+plot(layout = Plots.grid(10,10), size = (2300, 2800))
 for i = 1:100
     plot!(clims = (-0.002, 0.02))
     grid_vector_plot!(important_idx[i], i, rNGWF)
     if i in red_box_inds
         plot_square!(Nx, Ny; subplot = i)
+    elseif i in orange_box_inds
+        plot_square!(Nx, Ny; subplot = i, c = :orange)
     end
 end
 plt = current()

test/dissertations/htli/figs/Path512_HGLET_j2k1l5.png

@@ -4,26 +4,33 @@ gr(dpi = 200)
 ## (a) P64
 include("setups/path64.jl")
 ρ1 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
-plt = ROT_TSD_ratio_histogram(ρ1)
-savefig(plt, "../figs/Path64_ROT_TSD.png")
+# plt = ROT_TSD_ratio_histogram(ρ1)
+# savefig(plt, "../figs/Path64_ROT_TSD.png")
 
 ## (b) P₇ x P₃
 include("setups/grid7x3.jl")
 ρ2 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
-plt = ROT_TSD_ratio_histogram(ρ2)
-savefig(plt, "../figs/Grid7x3_ROT_TSD.png")
+# plt = ROT_TSD_ratio_histogram(ρ2)
+# savefig(plt, "../figs/Grid7x3_ROT_TSD.png")
 
 ## (c) Erdos Rényi
 include("setups/er.jl")
 ρ3 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
-plt = ROT_TSD_ratio_histogram(ρ3)
-savefig(plt, "../figs/Erdos_Renyi_ROT_TSD.png")
+# plt = ROT_TSD_ratio_histogram(ρ3)
+# savefig(plt, "../figs/Erdos_Renyi_ROT_TSD.png")
 
 ## (d) weighted RGC100
 include("setups/rgc100.jl")
 ρ4 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q; edge_length = edge_length)
-plt = ROT_TSD_ratio_histogram(ρ4)
-savefig(plt, "../figs/wRGC100_ROT_TSD.png")
+# plt = ROT_TSD_ratio_histogram(ρ4)
+# savefig(plt, "../figs/wRGC100_ROT_TSD.png")
+
+## boxplot
+plt = boxplot(["(a) Path" "(b) Grid" "(c) ER" "(d) RGC100"], [ρ1, ρ2, ρ3, ρ4];
+    legend = false, frame = :box, ylim = [0.9, 1.9], tickfontsize = 11,
+    outliers = true, grid = false, range = 3, lw = 1, size = (800, 600),
+    ylab = "ρ", yguidefontsize = 14)
+savefig(plt, "../figs/ROT_TSD_boxplots.png")
 
 ## Table 4.1
 display_basic_stats([ρ1, ρ2, ρ3, ρ4])

test/dissertations/htli/figs/Path512_LPHGLET_j2k1l5.png

@@ -1,10 +1,11 @@
 cd(@__DIR__); include("setups/simpletree.jl");
-gr(dpi=200)
+pyplot(dpi=200)
 
 for l = 97:100
+    local plt
     plot(size = (200, 500), framestyle = :none)
     gplot!(W, X, width = 1)
-    scatter_gplot!(X; marker = 𝚽[:, l])
+    scatter_gplot!(X; marker = 𝚽[:, l], ms = 5)
     plt = plot!(cbar = false, clims = (-0.3, 0.3))
     savefig(plt, "../figs/SimpleTree_eigenvector$(l-1)_red.png")
 end

test/dissertations/htli/figs/RGC100_HGLET_LPHGLET_j1k0l0.png

@@ -4,12 +4,12 @@ pyplot(dpi = 200)
 ## (a)
 plot(size = (200,500))
 gplot!(W, X, width = 1)
-scatter_gplot!(X[ib1, :]; c = :pink)
-scatter_gplot!(X[ib2, :]; c = :orange)
-scatter_gplot!(X[ib3, :]; c = :green)
-scatter_gplot!(X[ib4, :]; c = :yellow)
-scatter_gplot!(X[ijc, :]; c = :red)
-scatter_gplot!(X[ir, :]; c = :grey)
+scatter_gplot!(X[ib1, :]; c = :pink, ms = 5)
+scatter_gplot!(X[ib2, :]; c = :orange, ms = 5)
+scatter_gplot!(X[ib3, :]; c = :green, ms = 5)
+scatter_gplot!(X[ib4, :]; c = :yellow, ms = 5)
+scatter_gplot!(X[ijc, :]; c = :red, ms = 5)
+scatter_gplot!(X[ir, :]; c = :grey, ms = 5)
 plot!(framestyle = :none)
 plot!([0, 0], [0.3, 2], line = :arrow, c = :black, lw = 1)
 plt = annotate!(0, 3, text("root", 9))
@@ -19,12 +19,12 @@ savefig(plt, "../figs/SimpleTree.png")
 ## (b)
 plot(size = (200, 500))
 gplot!(Ws, Xs, width = 1)
-scatter_gplot!(Xs[11:11, :]; c = :pink, ms = 4)
-scatter_gplot!(Xs[10:10, :]; c = :orange, ms = 4)
-scatter_gplot!(Xs[13:13, :]; c = :green, ms = 4)
-scatter_gplot!(Xs[12:12, :]; c = :yellow, ms = 4)
-scatter_gplot!(Xs[[2,4,6,8], :]; c = :red, ms = 4)
-scatter_gplot!(Xs[[1,3,5,7,9], :]; c = :grey, ms = 4)
+scatter_gplot!(Xs[11:11, :]; c = :pink, ms = 5)
+scatter_gplot!(Xs[10:10, :]; c = :orange, ms = 5)
+scatter_gplot!(Xs[13:13, :]; c = :green, ms = 5)
+scatter_gplot!(Xs[12:12, :]; c = :yellow, ms = 5)
+scatter_gplot!(Xs[[2,4,6,8], :]; c = :red, ms = 5)
+scatter_gplot!(Xs[[1,3,5,7,9], :]; c = :grey, ms = 5)
 plt = plot!(framestyle = :none, xlims = [minimum(X[:, 1]), maximum(X[:, 1])],
             ylims = [minimum(X[:, 2]), maximum(X[:, 2])])
 savefig(plt, "../figs/SimpleTree_simplified.png")

test/dissertations/htli/figs/RGC100_HGLET_LPHGLET_j1k0l11.png

@@ -5,30 +5,30 @@ pyplot(dpi = 200)
 println("l = $(ieb1[end] - 1) and λₗ = $(𝛌[ieb1[end] - 1])")
 plot(size = (200, 500), framestyle = :none)
 gplot!(W, X, width = 1)
-scatter_gplot!(X; marker = 𝚽[:, ieb1[end]])
+scatter_gplot!(X; marker = 𝚽[:, ieb1[end]], ms = 5)
 plt = plot!(cbar = false, clims = (-0.3, 0.3))
 savefig(plt, "../figs/SimpleTree_eigenvector$(ieb1[end]-1)_pink.png")
 
 ## (b)
 println("l = $(ieb2[end] - 1) and λₗ = $(𝛌[ieb2[end] - 1])")
 plot(size = (200, 500), framestyle = :none)
 gplot!(W, X, width = 1)
-scatter_gplot!(X; marker = 𝚽[:, ieb2[end]])
+scatter_gplot!(X; marker = 𝚽[:, ieb2[end]], ms = 5)
 plt = plot!(cbar = false, clims = (-0.3, 0.3))
 savefig(plt, "../figs/SimpleTree_eigenvector$(ieb2[end]-1)_orange.png")
 
 ## (c)
 println("l = $(ieb3[end] - 1) and λₗ = $(𝛌[ieb3[end] - 1])")
 plot(size = (200, 500), framestyle = :none)
 gplot!(W, X, width = 1)
-scatter_gplot!(X; marker = 𝚽[:, ieb3[end]])
+scatter_gplot!(X; marker = 𝚽[:, ieb3[end]], ms = 5)
 plt = plot!(cbar = false, clims = (-0.3, 0.3))
 savefig(plt, "../figs/SimpleTree_eigenvector$(ieb3[end]-1)_green.png")
 
 ## (d)
 println("l = $(ieb4[end] - 1) and λₗ = $(𝛌[ieb4[end] - 1])")
 plot(size = (200, 500), framestyle = :none)
 gplot!(W, X, width = 1)
-scatter_gplot!(X; marker = 𝚽[:, ieb4[end]])
+scatter_gplot!(X; marker = 𝚽[:, ieb4[end]], ms = 5)
 plt = plot!(cbar = false, clims = (-0.3, 0.3))
 savefig(plt, "../figs/SimpleTree_eigenvector$(ieb4[end]-1)_yellow.png")

test/dissertations/htli/figs/RGC100_HGLET_LPHGLET_j2k0l2.png

@@ -0,0 +1,20 @@
+cd(@__DIR__); include("setups/path512.jl")
+pyplot(dpi = 200)
+
+## construct full HGLET and Lapped-HGLET dictionaries
+@time HGLET_dic = HGLET_dictionary(GP, G_Sig; method = :L)
+@time LPHGLET_dic = LPHGLET_dictionary(GP, G_Sig; method = :L, ϵ = 0.3)
+
+j = 3; k = 2; l = 6;
+WH = HGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
+WlH = LPHGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
+
+plt = plot(WH[:, l], c = :black, grid = false, frame = :box, lw = 0.8,
+     legendfontsize = 11, legend = false, size = (500, 400), ylim = [-0.15, 0.15])
+     xticks!([1; 64:64:N], vcat(string(1), [string(k) for k in 64:64:N]))
+savefig(plt, "../figs/Path512_HGLET_j$(j-1)k$(k-1)l$(l-1).png")
+
+plt = plot(WlH[:, l], c = :black, grid = false, frame = :box, lw = 0.8,
+     legendfontsize = 11, legend = false, size = (500, 400), ylim = [-0.15, 0.15])
+     xticks!([1; 64:64:N], vcat(string(1), [string(k) for k in 64:64:N]))
+savefig(plt, "../figs/Path512_LPHGLET_j$(j-1)k$(k-1)l$(l-1).png")

test/dissertations/htli/figs/RGC100_HGLET_LPHGLET_j2k0l4.png

@@ -10,17 +10,38 @@ gr(dpi = 200)
 T = [(2, 1, 1), (2, 1, 12), (3, 1, 3), (3, 1, 5),
     (4, 5, 6), (6, 23, 3), (7, 5, 2), (8, 15, 1)]
 
-for (j, k, l) in T
+XLs = [
+    [-50, 50], [-50, 50], [-110, 40], [-110, 40],
+    [-20, 50], [50, 120], [-80, -20], [-40, -10]
+]
+YLs = [
+    [-50, 50], [-50, 50], [-50, 100], [-50, 100],
+    [-50, 20], [-70, 0], [-90, -30], [-160, -130]
+]
+
+for i in 1:length(T)
+    (j, k, l) = T[i]
     𝚽H = HGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
     𝚽lH = LPHGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
     p1 = gplot(W, X; width = 1)
-        scatter_gplot!(X; marker = 𝚽H[:, l], plotOrder = :s2l)
+        scatter_gplot!(X; marker = 𝚽H[:, l], plotOrder = :s2l, ms = 5)
         plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = false)
         plot!(left_margin = 3mm, bottom_margin = 3mm, clims = (-0.05, 0.05))
     p2 = gplot(W, X; width = 1)
-        scatter_gplot!(X; marker = 𝚽lH[:, l], plotOrder = :s2l)
+        scatter_gplot!(X; marker = 𝚽lH[:, l], plotOrder = :s2l, ms = 5)
         plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = false)
         plot!(left_margin = 3mm, bottom_margin = 3mm, clims = (-0.05, 0.05))
-    plt = plot(p1, p2, size = (800, 400), layout = Plots.grid(1, 2))
+    plt = plot(p1, p2, size = (640, 300), layout = Plots.grid(1, 2),
+        xlims = XLs[i], ylims = YLs[i])
+        # xlims = [-120, 120], ylims = [-150, 150])
+    # gplot(W, X; width = 1)
+    # ord = :s2l 
+    # if i == 8
+    #     ord = :l2s
+    # end
+    # scatter_gplot!(X; marker = 𝚽lH[:, l] - 𝚽H[:, l], plotOrder = ord, ms = 4)
+    # plot!(xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = true)
+    # plot!(left_margin = 3mm, bottom_margin = 3mm, right_margin = 3mm)
+    # plt = plot!(size = (500, 500))
     savefig(plt, "../figs/RGC100_HGLET_LPHGLET_j$(j-1)k$(k-1)l$(l-1).png")
 end

test/dissertations/htli/figs/RGC100_HGLET_LPHGLET_j3k4l5.png

@@ -18,5 +18,5 @@ savefig(plt, "../figs/SunFlower_barbara_feye.png")
 DVEC = getall_expansioncoeffs2(G_Sig, GP_dual, GP_dual_Lsym, VM_NGWP, PC_NGWP, LP_NGWP,
                                VM_NGWP_Lsym, PC_NGWP_Lsym, LP_NGWP_Lsym, 𝚽, 𝚽sym)
 approx_error_plot2(DVEC);
-plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 11)
+plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 11, size = (600, 600))
 savefig(plt, "../figs/SunFlower_barbara_feye_DAG_approx.png")

test/dissertations/htli/figs/RGC100_HGLET_LPHGLET_j5k22l2.png

@@ -18,5 +18,5 @@ savefig(plt, "../figs/SunFlower_barbara_ftrouser.png")
 DVEC = getall_expansioncoeffs2(G_Sig, GP_dual, GP_dual_Lsym, VM_NGWP, PC_NGWP, LP_NGWP,
                                VM_NGWP_Lsym, PC_NGWP_Lsym, LP_NGWP_Lsym, 𝚽, 𝚽sym)
 approx_error_plot2(DVEC);
-plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 11)
+plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 11, size = (600, 600))
 savefig(plt, "../figs/SunFlower_barbara_ftrouser_DAG_approx.png")