Function for calculating rank

helper.jl

@@ -71,3 +71,11 @@ function gram_schmidt(vecs, ws, precision=1e-10, verbose=false)
     end
     return vecs[mask]
 end
+
+"Rank of a basis set (pre-normalized) under c-product"
+function c_rank(vecs, ws, threshold=1e-8)
+    basis = hcat(vecs...)
+    N = transpose(basis) * spdiagm(ws) * basis
+    abs_singular_values = eigvals(N) .|> abs .|> sqrt
+    return count(abs_singular_values .> threshold)
+end