stick with previous way and comment the alternative for future use

Author: haotian127

src/LP-NGWP.jl

@@ -230,17 +230,19 @@ end
 
 function const_meyer_wavelets(𝚽, Uf; idx = 1:size(Uf, 1))
     N = size(𝚽, 1)
-    # # assemble smooth orthogonal projector custom to nodes `idx`
-    # P = Uf[idx, :]' * Uf[idx, idx]
-    # if diag(P) == χ(idx, N)
-    #     B = 𝚽[:, idx]
-    # else
-    #     # folding the eigenspace, i.e., 𝚽's column space
-    #     Y = 𝚽 * P
-    #     # find its column space's orthogonal basis
-    #     B = svd(Y).U
-    # end
-    B = 𝚽 * Uf[:, idx]
+    # 1. nice interpretable way by smooth orthogonal projector
+    # assemble smooth orthogonal projector custom to nodes `idx`
+    P = Uf[idx, :]' * Uf[idx, idx]
+    if diag(P) == χ(idx, N)
+        B = 𝚽[:, idx]
+    else
+        # folding the eigenspace, i.e., 𝚽's column space
+        Y = 𝚽 * P
+        # find its column space's orthogonal basis
+        B = svd(Y).U
+    end
+    # 2. alternative faster way by orthogonal unfolding operator
+    # B = 𝚽[:, idx] * Uf[idx, idx]
     # perform varimax rotation to get the meyer_wavelets
     Wavelets = varimax(B)
     return Wavelets