add weight structure computation on variables

Author: vitktor

src/on_varieties/lie-symmetries/lie-algebras/actions.jl

@@ -62,40 +62,34 @@ function act(elem::LieAlgebraElem{ScalingLieAlgebra}, f::Union{Expression, Monom
     return expand(dot(differentiate(f, variables(action)), -X*variables(action)))
 end
 
-function weight_structure(a::ScalingLieAction, vars::Vector{Variable})
-    
+function inv_weight_space(a::AbstractLieAlgebraAction, vars::Vector{Variable})
+    inv_vars = setdiff(vars, variables(a))
+    isempty(inv_vars) && return nothing
+    ws = WeightSpace(zeros(Int, rank(algebra(a))), zeros(ComplexF64, length(vars), length(inv_vars)))
+    for (i, var) in enumerate(inv_vars)
+        ws[vars_dict[var], i] = 1
+    end
+    return ws
 end
 
-function hw_spaces(a::ScalingLieAction, vars::Vector{Variable})
+function weight_structure(a::ScalingLieAction, vars::Vector{Variable})
     @assert variables(a) ⊆ vars
     vars_dict = Dict(zip(vars, 1:length(vars)))
-    hw_spaces_dict = Dict{Vector{Int}, WeightSpace}()
+    hw_struct = WeightStructure()
     for (i, var) in enumerate(variables(a))
         weight = exponents(a)[:, i]
-        if haskey(hw_spaces_dict, weight)
-            sp = zeros(ComplexF64, length(vars))
-            sp[vars_dict[var]] = 1
-            hw_spaces_dict[weight].space = hcat(hw_spaces_dict[weight].space, sp)
-        else
-            hw_space = WeightSpace(weight, zeros(ComplexF64, length(vars)))
-            hw_space.space[vars_dict[var]] = 1
-            hw_spaces_dict[weight] = hw_space
-        end
-    end
-    inv_vars = setdiff(vars, variables(a))
-    isempty(inv_vars) && return collect(values(hw_spaces_dict))
-    inv_hw_space = WeightSpace(zeros(Int, rank(algebra(a))), zeros(ComplexF64, length(vars), length(inv_vars)))
-    for (i, var) in enumerate(inv_vars)
-        inv_hw_space.space[vars_dict[var], i] = 1
-    end
-    if haskey(hw_spaces_dict, weight(inv_hw_space))
-        hw_spaces_dict[weight(inv_hw_space)].space = hcat(hw_spaces_dict[weight(inv_hw_space)].space, inv_hw_space.space)
-    else
-        hw_spaces_dict[weight(inv_hw_space)] = inv_hw_space
+        hwv = WeightVector(weight, zeros(ComplexF64, length(vars)))
+        hwv[vars_dict[var]] = 1
+        push!(hw_struct, hwv)
     end
-    return collect(values(hw_spaces_dict))
+    inv_ws = inv_weight_space(a, vars)
+    !isnothing(inv_ws) && push!(hw_struct, inv_ws)
+    return hw_struct
 end
 
+hw_spaces(a::ScalingLieAction, vars::Vector{Variable}) = weight_structure(a, vars)
+
+
 struct LieAlgebraAction <: AbstractLieAlgebraAction
     alg::LieAlgebra
     var_groups::Vector{Vector{Variable}}
@@ -129,38 +123,27 @@ function act(elem::LieAlgebraElem{LieAlgebra}, f::Union{Expression, Monomial}, a
     return expand(sum([dot(differentiate(f, vars), -X*vars) for vars in var_groups(action)]))
 end
 
-function weight_structure(a::LieAlgebraAction, vars::Vector{Variable})
-    
-end
-
-function hw_spaces(a::LieAlgebraAction, vars::Vector{Variable})
+function weight_structure(a::LieAlgebraAction, vars::Vector{Variable}; as_hw_spaces::Bool=false)
     @assert variables(a) ⊆ vars
     vars_dict = Dict(zip(vars, 1:length(vars)))
-    hw_spaces_vec = WeightSpace[]
-    for hw_space in hw_spaces(algebra(a))
-        hws = WeightSpace(hw_space.weight, zeros(ComplexF64, length(vars), dim(hw_space)*nvar_groups(a)))
+    ws = WeightStructure()
+    alg_ws = as_hw_spaces ? hw_spaces(algebra(a)) : weight_structure(algebra(a))
+    for w_space in alg_ws
+        new_w_space = WeightSpace(weight(w_space), zeros(ComplexF64, length(vars), dim(w_space)*nvar_groups(a)))
         for (i, vars_group) in enumerate(var_groups(a))
             for (j, var) in enumerate(vars_group)
-                hws.space[vars_dict[var], (i-1)*dim(hw_space)+1:i*dim(hw_space)] = hw_space.space[j, :]
+                new_w_space[vars_dict[var], (i-1)*dim(w_space)+1:i*dim(w_space)] = w_space[j, :]
             end
         end
-        push!(hw_spaces_vec, hws)
-    end
-    inv_vars = setdiff(vars, variables(a))
-    isempty(inv_vars) && return hw_spaces_vec
-    inv_hw_space = WeightSpace(zeros(Int, rank(algebra(a))), zeros(ComplexF64, length(vars), length(inv_vars)))
-    for (i, var) in enumerate(inv_vars)
-        inv_hw_space.space[vars_dict[var], i] = 1
-    end
-    hws_dict = Dict(zip([weight(hws) for hws in hw_spaces_vec], hw_spaces_vec))
-    if haskey(hws_dict, weight(inv_hw_space))
-        hws_dict[weight(inv_hw_space)].space = hcat(hws_dict[weight(inv_hw_space)].space, inv_hw_space.space)
-    else
-        push!(hw_spaces_vec, inv_hw_space)
+        push!(ws, new_w_space)
     end
-    return hw_spaces_vec
+    inv_ws = inv_weight_space(a, vars)
+    !isnothing(inv_ws) && push!(ws, inv_ws)
+    return ws
 end
 
+hw_spaces(a::LieAlgebraAction, vars::Vector{Variable}) = weight_structure(a, vars, as_hw_spaces=true)
+
 
 struct SumLieAlgebraAction <: AbstractLieAlgebraAction
     alg::SumLieAlgebra

src/on_varieties/lie-symmetries/lie-reprs/representation.jl

@@ -1,5 +1,5 @@
 export LieAlgebraRepresentation,
-    space,
+    w_space,
     isotypic_components,
     highest_weights
 
@@ -30,11 +30,11 @@ end
 
 action(ρ::LieAlgebraRepresentation) = ρ.action
 algebra(ρ::LieAlgebraRepresentation) = algebra(ρ.action)
-space(ρ::LieAlgebraRepresentation) = ρ.V
+w_space(ρ::LieAlgebraRepresentation) = ρ.V
 nisotypic(ρ::LieAlgebraRepresentation) = length(ρ.isotypic)
 isotypic_components(ρ::LieAlgebraRepresentation) = ρ.isotypic
 irreducibles(ρ::LieAlgebraRepresentation) = vcat([irreducibles(iso) for iso in ρ.isotypic]...)
-dim(ρ::LieAlgebraRepresentation) = dim(space(ρ))
+dim(ρ::LieAlgebraRepresentation) = dim(w_space(ρ))
 highest_weights(ρ::LieAlgebraRepresentation) = [highest_weight(ic) for ic in isotypic_components(ρ)]
 
 # called by Shift+Enter
@@ -64,8 +64,8 @@ function nullspace_as_weight_vectors(
     tol::Real=1e-5
 ) where {T <: Number}
     wvs = WeightVector[]
-    for weight_space in weight_spaces(ws)
-        N = space(weight_space)*nullspace(M*space(weight_space))
+    for weight_space in ws
+        N = matrix(weight_space)*nullspace(M*matrix(weight_space))
         sparsify!(N, tol)
         vs = M2VV(N; copy=false)
         for v in vs

src/on_varieties/lie-symmetries/scalings.jl

@@ -3,7 +3,7 @@ export Grading,
     grading,
     scaling_symmetries
 
-using AbstractAlgebra: ZZ, matrix, GF, lift, hnf, snf_with_transform 
+using AbstractAlgebra: ZZ, matrix as aa_matrix, GF, lift, hnf, snf_with_transform 
 using LinearAlgebra: diag
 
 mutable struct Grading{Tv<:Integer,Ti<:Integer}
@@ -162,7 +162,7 @@ function _snf_scaling_symmetries(F::System)::Tuple{Vector{Int}, Vector{Matrix{In
     if size(K, 1) > size(K, 2)
         K = [K zeros(eltype(K), size(K, 1), size(K,1)-size(K,2))]
     end
-    K = matrix(ZZ, K)
+    K = aa_matrix(ZZ, K)
 
     S, U, _ = snf_with_transform(K)
     U, S = Matrix(U), Int.(diag(Matrix(S)))
@@ -182,10 +182,10 @@ function _hnf_reduce!(
 )
     for (i, (sᵢ, Uᵢ)) in enumerate(zip(s, U))
         if sᵢ == 0
-            U[i] = Matrix(hnf(matrix(ZZ, Uᵢ)))
+            U[i] = Matrix(hnf(aa_matrix(ZZ, Uᵢ)))
         else
             try
-                U[i] = Int.(lift.(Matrix(hnf(matrix(GF(sᵢ), Uᵢ)))))
+                U[i] = Int.(lift.(Matrix(hnf(aa_matrix(GF(sᵢ), Uᵢ)))))
             catch
                 U[i] = mod.(Uᵢ, sᵢ)
             end
@@ -198,12 +198,12 @@ function _hnf_reduce(grading::Grading{Tv,Ti}) where {Tv<:Integer,Ti<:Integer}
     U₀ = grading.free_part
     red_grading = Grading{Tv,Ti}()
     if !isnothing(U₀)
-        red_grading.free_part = Matrix(hnf(matrix(ZZ, U₀)))
+        red_grading.free_part = Matrix(hnf(aa_matrix(ZZ, U₀)))
         red_grading.nscalings = size(U₀, 1)
     end
     for (sᵢ, Uᵢ) in grading.mod_part
         try
-            Uᵢ = lift.(Matrix(hnf(matrix(GF(sᵢ), Uᵢ))))
+            Uᵢ = lift.(Matrix(hnf(aa_matrix(GF(sᵢ), Uᵢ))))
         catch
             Uᵢ = mod.(Uᵢ, sᵢ)
         end

src/on_varieties/lie-symmetries/weights.jl

@@ -22,20 +22,33 @@ end
 
 weight(wv::WeightVector) = wv.weight
 vector(wv::WeightVector) = wv.vector
+Base.setindex!(wv::WeightVector, n::Number, i::Integer) = wv.vector[i] = n
 
 mutable struct WeightSpace
     weight::Vector{Int}
-    space::Matrix{ComplexF64}
+    matrix::Matrix{ComplexF64}
 end
 
 WeightSpace(wv::WeightVector) = WeightSpace(weight(wv), V2M(vector(wv)))
-WeightSpace(weight::AbstractVector, space::AbstractVector) = WeightSpace(weight, V2M(space))
+WeightSpace(weight::AbstractVector, v::AbstractVector) = WeightSpace(weight, V2M(v))
 
 weight(ws::WeightSpace) = ws.weight
 weight_size(ws::WeightSpace) = length(weight(ws))
-space(ws::WeightSpace) = ws.space
-Base.size(ws::WeightSpace) = size(ws.space, 1)
-dim(ws::WeightSpace) = size(ws.space, 2)
+matrix(ws::WeightSpace) = ws.matrix
+Base.size(ws::WeightSpace) = size(ws.matrix, 1)
+dim(ws::WeightSpace) = size(ws.matrix, 2)
+Base.getindex(ws::WeightSpace, inds...) = ws.matrix[inds...]
+Base.setindex!(ws::WeightSpace, M::AbstractVecOrMat{T}, inds...) where {T <: Number} = ws.matrix[inds...] = M
+
+function Base.hcat(ws1::WeightSpace, ws2::WeightSpace)
+    if weight(ws1) != weight(ws2)
+        throw(ArgumentError("Weights of the two spaces must be equal"))
+    end
+    if size(ws1) != size(ws2)
+        throw(ArgumentError("Sizes of the two spaces must be equal"))
+    end
+    return WeightSpace(weight(ws1), hcat(matrix(ws1), matrix(ws2)))
+end
 
 function Base.show(io::IO, ::MIME"text/plain", ws::WeightSpace)
     println(io, "WeightSpace of dimension $(dim(ws))")
@@ -48,10 +61,10 @@ end
 
 function Base.iterate(ws::WeightSpace, state=1)
     state > dim(ws) && return nothing
-    return (WeightVector(weight(ws), space(ws)[:, state]), state+1)
+    return (WeightVector(weight(ws), matrix(ws)[:, state]), state+1)
 end
 
-function intersect_spaces(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where {T <: Number}
+function Base.:∩(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where {T <: Number}
     N = nullspace(hcat(A, B))
     N = hcat([div_by_lowest_magnitude(N[:,i], 1e-8) for i in 1:size(N, 2)]...)
     sparsify!(N, 1e-8)
@@ -60,30 +73,23 @@ end
 
 function Base.:∩(ws1::WeightSpace, ws2::WeightSpace)
     new_weight = vcat(weight(ws1), weight(ws2))
-    new_space = intersect_spaces(space(ws1), space(ws2))
-    iszero(size(new_space, 2)) && return nothing
-    return WeightSpace(new_weight, new_space)
-end
-
-function Base.:∩(ws1::Vector{WeightSpace}, ws2::Vector{WeightSpace})
-    new_spaces = [ws1[i] ∩ ws2[j] for i in 1:length(ws1), j in 1:length(ws2)]
-    return [ws for ws in new_spaces if !isnothing(ws)]
+    new_matrix = ∩(matrix(ws1), matrix(ws2))
+    iszero(size(new_matrix, 2)) && return nothing
+    return WeightSpace(new_weight, new_matrix)
 end
 
-Base.:∩(ws::Vector{Vector{WeightSpace}}) = reduce(∩, ws)
-
 
 struct WeightStructure
-    weight_spaces::Vector{WeightSpace}
-    dict::Dict{Vector{Int}, Matrix{ComplexF64}}
+    weights::Vector{Vector{Int}}
+    dict::Dict{Vector{Int}, WeightSpace} # weight => space
 end
 
 WeightStructure(
     weights::Vector{Vector{Int}},
     weight_spaces::Vector{Matrix{T}} where T <: Number
 ) = WeightStructure(
-    [WeightSpace(w, M) for (w, M) in zip(weights, weight_spaces)],
-    Dict(zip(weights, weight_spaces))
+    weights,
+    Dict(zip(weights, [WeightSpace(w, M) for (w, M) in zip(weights, weight_spaces)]))
 )
 
 WeightStructure(
@@ -103,20 +109,30 @@ WeightStructure(
     weight_vectors::Vector{Vector{T}} where T <: Number
 ) = WeightStructure(weights, [V2M(v) for v in weight_vectors])
 
+function WeightStructure(w_spaces::Vector{WeightSpace})
+    ws = WeightStructure()
+    for w_space in w_spaces
+        push!(ws, w_space)
+    end
+    return ws
+end
+
 weight_spaces(
     ws::WeightStructure;
     as_matrices::Bool=false
-) = as_matrices ? [space(wsp) for wsp in ws.weight_spaces] : ws.weight_spaces
+) = as_matrices ? [matrix(ws[w]) for w in weights(ws)] : [ws[w] for w in weights(ws)]
+
 weight_spaces(
     ws::WeightStructure,
     inds...;
     as_matrices::Bool=false
 ) = weight_spaces(ws; as_matrices=as_matrices)[inds...]
-weight_size(ws::WeightStructure) = weight_size(ws[1])
-space_dim(ws::WeightStructure) = size(ws[1])
-weights(ws::WeightStructure) = [weight(ws) for ws in weight_spaces(ws)]
-weight(ws::WeightStructure, i::Integer) = weight(weight_spaces(ws)[i])
-nweights(ws::WeightStructure) = length(weight_spaces(ws))
+
+weights(ws::WeightStructure) = ws.weights
+weight_size(ws::WeightStructure) = length(first(weights(ws)))
+space_dim(ws::WeightStructure) = size(ws[first(weights(ws))])
+weight(ws::WeightStructure, i::Integer) = weights(ws)[i]
+nweights(ws::WeightStructure) = length(weights(ws))
 dims(ws::WeightStructure) = [dim(M) for M in weight_spaces(ws)]
 
 function Base.show(io::IO, ws::WeightStructure)
@@ -126,30 +142,46 @@ function Base.show(io::IO, ws::WeightStructure)
     print(io, " max weight space dimension: ", maximum(dims(ws)))
 end
 
-Base.getindex(ws::WeightStructure, i::Integer) = ws.weight_spaces[i]
+Base.length(ws::WeightStructure) = nweights(ws)
 Base.getindex(ws::WeightStructure, weight::Vector{Int}) = ws.dict[weight]
+Base.getindex(ws::WeightStructure, i::Integer) = ws[weight(ws, i)]
+Base.setindex!(ws::WeightStructure, ws_new::WeightSpace, weight::Vector{Int}) = ws.dict[weight] = ws_new
 Base.haskey(ws::WeightStructure, weight::Vector{Int}) = haskey(ws.dict, weight)
 
-function Base.push!(ws::WeightStructure, wv::WeightVector)
-    if haskey(ws, weight(wv))
-        ws.dict[weight(wv)] = hcat(ws[weight(wv)], space(wv))
+function Base.iterate(ws::WeightStructure, state=1)
+    state > nweights(ws) && return nothing
+    return (ws[state], state+1)
+end
+
+function Base.push!(ws::WeightStructure, w_space::WeightSpace)
+    if haskey(ws, weight(w_space))
+        ws[weight(w_space)] = hcat(ws[weight(w_space)], w_space)
     else
-        wsp = WeightSpace(wv)
-        push!(ws.weight_spaces, wsp)
-        ws.dict[weight(wv)] = space(wsp)
+        push!(weights(ws), weight(w_space))
+        ws[weight(w_space)] = w_space
     end
+    return ws
 end
 
-function sum_weight_structure(ws::WeightStructure, d::Integer)
-    new_weight_spaces = [zeros(ComplexF64, space_dim(ws)*d, size(M, 2)*d) for M in weight_spaces(ws)]
-    for (M, new_M) in zip(weight_spaces(ws), new_weight_spaces)
-        for j in 1:d
-            new_M[(1:space_dim(ws)) .+ (j-1)*space_dim(ws), (1:size(M,2)) .+ (j-1)*size(M,2)] = M
-        end
-    end
-    return WeightStructure(weights(ws), new_weight_spaces)
+Base.push!(ws::WeightStructure, wv::WeightVector) = push!(ws, WeightSpace(wv))
+
+function Base.:∩(ws1::WeightStructure, ws2::WeightStructure)
+    new_spaces = [w_space₁ ∩ w_space₂ for w_space₁ in ws1, w_space₂ in ws2]
+    return WeightStructure([ws for ws in new_spaces if !isnothing(ws)])
 end
 
+Base.:∩(ws::Vector{WeightStructure}) = reduce(∩, ws)
+
+# function sum_weight_structure(ws::WeightStructure, d::Integer)
+#     new_weight_spaces = [zeros(ComplexF64, space_dim(ws)*d, size(M, 2)*d) for M in weight_spaces(ws)]
+#     for (M, new_M) in zip(weight_spaces(ws), new_weight_spaces)
+#         for j in 1:d
+#             new_M[(1:space_dim(ws)) .+ (j-1)*space_dim(ws), (1:size(M,2)) .+ (j-1)*size(M,2)] = M
+#         end
+#     end
+#     return WeightStructure(weights(ws), new_weight_spaces)
+# end
+
 function sym_weight_structure(ws::WeightStructure, d::Integer, mexps::Vector{<:SparseVector})
     d = Int(d)
     d == 0 && return WeightStructure([zeros(Int, weight_size(ws))], [[1;;]])
@@ -241,5 +273,4 @@ coefficients(we::WeightExpression) = we.coeffs
 basis(we::WeightExpression) = we.basis
 weight(we::WeightExpression) = we.weight
 expression(we::WeightExpression) = dot(coefficients(we), to_expressions(basis(we)))
-
 WeightVector(we::WeightExpression) = WeightVector(weight(we), coefficients(we))