refactor: improve error handling and logging in nullspace functions

Author: vitktor

src/decompositions/irred-decomp.jl

@@ -15,13 +15,13 @@ function ws_nullspace(
     if all(f -> isapprox(f, zero(f), atol=tol), Bₐ)
         return ws
     end
-    dim(ws) == 1 && return nothing
+    dim(ws) == 1 && return nothing # implied by lines 15-17
     Cs = zero_combinations(Bₐ; tol=tol, logging=logging)
     isempty(Cs) && return nothing
     return WeightSpace(
         weight(ws),
         PolySpace{field_type(ws)}([div_by_smallest_coeff(sum(c .* B)) for c in Cs])
-        )
+    )
 end
 
 function common_nullspace_as_weight_vectors(
@@ -38,16 +38,16 @@ function common_nullspace_as_weight_vectors(
     for X in Xs
         wsₙ = ws_nullspace(X, a, wsₙ)
         if isnothing(wsₙ)
-            logging && printstyled("Found no common nullspace, smth is wrong!\n", color=:red)
+            error("Found no common nullspace, while looking for $(nhwv) highest weight vectors; dim(ws) = $(dim(ws))")
             return WeightVector[] # error("Didn't find any highest weight vectors, had to find $(nhwv)")
         end
     end
 
     if dim(wsₙ) != nhwv
         error("Found $(dim(wsₙ)) highest weight vectors, expected $(nhwv)")
     end
-    logging && printstyled("Found: ", color=:green)
-    logging && printstyled("$(join(map(repr, basis(space(wsₙ))), ", "))\n", color=:green)
+    # logging && printstyled("Found: ", color=:green)
+    # logging && printstyled("$(join(map(repr, basis(space(wsₙ))), ", "))\n", color=:green)
 
     return [WeightVector(weight(wsₙ), div_by_smallest_coeff(v)) for v in basis(space(wsₙ))] # TODO: remove div_by_smallest_coeff
 end
@@ -109,11 +109,19 @@ function irreducibles(
 
     logging && printstyled("Applying Clebsch-Gordan to decompose $(to_string(V)) into irreducibles...\n", color=:blue)
     sym_ws, cg_decomp = sym_weight_struct(highest_weight(Vb), algebra(action(ρ)), power(V), ws)
+    @assert Set(weights(sym_ws)) == keys(cg_decomp)
+
     logging && printstyled(to_string(V), " = ", str_irr_decomp("V", cg_decomp), "\n", color=:blue)
 
     Xs = positive_root_elements(algebra(action(ρ)))
     hw_vectors = common_nullspace_as_weight_vectors(Xs, action(ρ), sym_ws, cg_decomp; logging=logging)
-    return [IrreducibleRepresentation(action(ρ), hwv) for hwv in hw_vectors]
+
+    irrs = [IrreducibleRepresentation(action(ρ), hwv) for hwv in hw_vectors]
+    logging && println("sum = ", sum([dim(space(irr)) for irr in irrs]))
+    logging && println("dim(V): ", dim(V), ", dim(Vb): ", dim(Vb), ", power(V): ", power(V))
+    logging && println("dim(cg_decomp): ", weyl_dim(cg_decomp, algebra(action(ρ))))
+    @assert sum([dim(space(irr)) for irr in irrs]) == dim(V)
+    return irrs
 end
 
 function irreducibles(
@@ -148,6 +156,8 @@ function irreducibles(
         ρᵢ = GroupRepresentation(action(ρ), Vᵢ)
         append!(irreds, irreducibles(ρᵢ; logging=logging, decomp_count=decomp_count))
     end
+
+    @assert sum([dim(space(irr)) for irr in irreds]) == dim(V)
     return irreds
 end
 
@@ -173,6 +183,7 @@ function irreducibles(
 
     logging && printstyled("Applying Clebsch-Gordan to decompose $(to_string(V)) into irreducibles...\n", color=:blue)
     cg_decomp = tensor(hw₁, hw₂, algebra(action(ρ)))
+
     if length(cg_decomp) == 1
         logging && printstyled("The tensor product is irreducible, returning the tensor product of the highest weight vectors\n", color=:red)
         return [IrreducibleRepresentation(action(ρ), hw_vector(V₁)*hw_vector(V₂))]
@@ -184,7 +195,10 @@ function irreducibles(
 
     Xs = positive_root_elements(algebra(action(ρ)))
     hw_vectors = common_nullspace_as_weight_vectors(Xs, action(ρ), tensor_ws, cg_decomp; logging=logging)
-    return [IrreducibleRepresentation(action(ρ), hwv) for hwv in hw_vectors]
+
+    irrs = [IrreducibleRepresentation(action(ρ), hwv) for hwv in hw_vectors]
+    @assert sum([dim(space(irr)) for irr in irrs]) == dim(V)
+    return irrs
 end
 
 function irreducibles(
@@ -216,12 +230,13 @@ function irreducibles(
 
     irreds = IrreducibleRepresentation{A}[]
     for irr₁ in irrs₁, irr₂ in irrs₂
-        V = TensorProduct([space(irr₁), space(irr₂)])
-        ρV = GroupRepresentation(action(ρ), V)
+        Vt = TensorProduct([space(irr₁), space(irr₂)])
+        ρV = GroupRepresentation(action(ρ), Vt)
         irrs = irreducibles(ρV; logging=logging, decomp_count=decomp_count)
-        @assert sum([dim(space(irr)) for irr in irrs]) == dim(V)
+        @assert sum([dim(space(irr)) for irr in irrs]) == dim(Vt)
         append!(irreds, irrs)
     end
+    @assert sum([dim(space(irr)) for irr in irreds]) == dim(V)
     return irreds
 end
 

src/utils/basic.jl

@@ -3,7 +3,7 @@ export rref, monomials, coeffs_matrix, polynomials
 export num_mons
 export superscript, subscript
 export multiexponents, div_by_smallest_coeff
-export eye, a2p, p2a, rand_unit
+export eye, a2p, p2a, rand_unit, xx
 
 a2p(M::AbstractMatrix{<:Number}) = [M; ones(eltype(M), 1, size(M, 2))]
 p2a(M::AbstractMatrix{<:Number}) = (M./M[end:end,:])[1:end-1,:]
@@ -17,6 +17,7 @@ Base.rand(::Type{Rational{T}}) where T <: Integer = rand(T(-10):T(10)) // rand(T
 Base.rand(::Type{Complex{Rational{T}}}) where T <: Integer = rand(Rational{T}) + rand(Rational{T})*im
 Base.rand(T::DataType, n::Int) = [rand(T) for _ in 1:n]
 
+xx(v) = [0 -v[3] v[2]; v[3] 0 -v[1]; -v[2] v[1] 0]
 eye(n::Integer) = eye(Float64, n)
 eye(T::Type, n::Integer) = Matrix{T}(I(n))
 
@@ -161,11 +162,16 @@ function rref(
     F::Vector{<:AbstractPolynomial{T}};
     tol::Real=1e-5
 ) where T
+    printstyled("Computing rref, merging monomials...\n", color=:blue)
     mons = monomials(F)
+    printstyled("Monomials merged\n", color=:blue)
+    printstyled("Constructing coefficients matrix...\n", color=:blue)
     M = Matrix(transpose(coeffs_matrix(F, mons)))
+    printstyled("Coefficients matrix constructed\n", color=:blue)
     rref!(M, tol)
     sparsify!(M, tol)
     N = filter(row -> !iszero(row), eachrow(M))
+    printstyled("RREF computed\n", color=:blue)
     return polynomials(N, mons)
 end
 
@@ -183,17 +189,21 @@ function zero_combinations(
         M = rand(T, size(M, 2), size(M, 1)) * M
     end
     logging && print(crayon"#f4d03f", "Computing nullspace of $(size(M)) matrix...\n")
-    N = Matrix(transpose(nullspace(M; atol=tol)))
-    logging && print(crayon"#f4d03f", "Nullspace is $(size(N, 1))-dimensional\n")
+    Nc = Matrix(transpose(nullspace(M; atol=tol)))
+    logging && print(crayon"#f4d03f", "Nullspace is $(size(Nc, 1))-dimensional\n")
+    N = copy(Nc)
     sparsify!(N, tol)
     rref!(N, tol)
     for n in N
         if norm(n) > 1e3
+            display(Nc)
             display(N)
             error("Nullspace contains large coefficients")
         end
     end
     sparsify!(N, tol)
+    # println("Nc: ", "max = $(maximum(abs, Nc; init=-Inf))", ", min = $(minimum(abs, Nc; init=Inf))")
+    # println("N: ", "max = $(maximum(abs, N; init=-Inf))", ", min = $(minimum(abs, N; init=Inf))")
     return eachrow(N)
 end