From 0a474ce80a54c6a9feb02df88ab288a8e0b73170 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Mon, 22 May 2023 17:58:17 +0200 Subject: [PATCH 01/15] Code from BasisLieHighestWeight, Pull request oscar-system/Oscar.jl#2115 --- .../BasisLieHighestWeight/docs/doc.main | 3 + .../docs/src/basisLieHighestWeight.md | 16 + .../src/BasisLieHighestWeight.jl | 507 ++++++++++++++++++ .../BasisLieHighestWeight/src/LieAlgebras.jl | 54 ++ .../src/MonomialOrder.jl | 18 + .../BasisLieHighestWeight/src/NewMonomial.jl | 101 ++++ .../src/RootConversion.jl | 366 +++++++++++++ .../BasisLieHighestWeight/src/TensorModels.jl | 58 ++ .../src/VectorSpaceBases.jl | 75 +++ .../BasisLieHighestWeight/src/WeylPolytope.jl | 158 ++++++ .../BasisLieHighestWeight/test/MBOld.jl | 301 +++++++++++ .../BasisLieHighestWeight/test/runtests.jl | 86 +++ 12 files changed, 1743 insertions(+) create mode 100644 experimental/BasisLieHighestWeight/docs/doc.main create mode 100644 experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md create mode 100644 experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl create mode 100644 experimental/BasisLieHighestWeight/src/LieAlgebras.jl create mode 100644 experimental/BasisLieHighestWeight/src/MonomialOrder.jl create mode 100644 experimental/BasisLieHighestWeight/src/NewMonomial.jl create mode 100644 experimental/BasisLieHighestWeight/src/RootConversion.jl create mode 100644 experimental/BasisLieHighestWeight/src/TensorModels.jl create mode 100644 experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl create mode 100644 experimental/BasisLieHighestWeight/src/WeylPolytope.jl create mode 100644 experimental/BasisLieHighestWeight/test/MBOld.jl create mode 100644 experimental/BasisLieHighestWeight/test/runtests.jl diff --git a/experimental/BasisLieHighestWeight/docs/doc.main b/experimental/BasisLieHighestWeight/docs/doc.main new file mode 100644 index 000000000000..d4c956dfab8b --- /dev/null +++ b/experimental/BasisLieHighestWeight/docs/doc.main @@ -0,0 +1,3 @@ +[ + "basisLieHighestWeight.md" +] diff --git a/experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md b/experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md new file mode 100644 index 000000000000..7fc64f8ce43c --- /dev/null +++ b/experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md @@ -0,0 +1,16 @@ +```@meta +CurrentModule = Oscar.BasisLieHighestWeight +``` + +```@setup oscar +using Oscar.BasisLieHighestWeight +``` + +```@contents +Pages = ["basisLieHighestWeight.md"] +``` + +# Monomial bases for Lie algebras +```@docs +basisLieHighestWeight2 +``` diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl new file mode 100644 index 000000000000..40669c734a1f --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -0,0 +1,507 @@ +module BasisLieHighestWeight +export basis_lie_highest_weight +export is_fundamental + +using Polymake + +include("./NewMonomial.jl") + +fromGap = Oscar.GAP.gap_to_julia + +@doc """ + basisLieHighestWeight(type::String, rank::Int, highest_weight::Vector{Int}; + operators::Union{String, Vector{Int}} = "regular", + monomial_order::Union{String, Function} = "GRevLex", cache_size::Int = 0, + parallel::Bool = false, return_no_minkowski::Bool = false, + return_operators::Bool = false) + +Computes a monomial basis for the highest weight module with highest weight +``highest_weight`` (in terms of the fundamental weights), for a simple Lie algebra of type +``type`` and rank ``rank``. + +# Parameters +- `type`: type of liealgebra we want to investigate, one of "A", "B", "C", "D", "E", "F", "G" +- `rank`: rank of liealgebra +- `highest_weight`: highest-weight +- `operators`: list of operators, either "regular" or integer array. The functionality of choosing a random longest word + is currently not implemented, because we used https://github.com/jmichel7/Gapjm.jl to work with coxeter + groups need a method to obtain all non left descending elements to extend a word +- `monomial_order`: monomial order in which our basis gets defined with regards to our operators +- `cache_size`: number of computed monomials we want to cache, default is 0 +- `parallel`: currently not implemented, because we used Distributed.jl, but if true parts of the algorithms can be + parallelized +- `return_no_minkowski`: if true return monomials for which Monkowski-property did not suffice to find all monomials +- `return_operators`: if true return the GAP objects operators + +# Examples +```jldoctest +julia> BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1, 1], return_no_minkowski = true, return_operators = true) +(Set(ZZMPolyRingElem[x1*x2, x2, 1, x1*x3, x3^2, x1, x3, x2*x3]), Set([[1, 0], [0, 1]]), GAP: [ v.1, v.2, v.3 ]) + + +julia> BasisLieHighestWeight.basis_lie_highest_weight("A", 3, [2, 2, 3], monomial_order = "Lex") +Set{ZZMPolyRingElem} with 1260 elements: + x3*x5^2*x6^2 + x2*x3*x5^2*x6 + x4^2*x5^2*x6 + x1^2*x3^3*x5 + x2*x3*x4*x5^3*x6^2 + x1*x3*x4*x5^3*x6^2 + x1^2*x3*x4*x6 + x1*x3*x4^3 + x4^2*x5^3*x6^4 + x1*x2*x3*x5^2 + x3^2*x4^4*x5^2*x6 + x2^2*x3*x6^2 + x1*x2^2*x3^2*x5 + x1*x3*x4*x5^2 + x1^2*x2*x6 + x1*x3^2*x4*x5*x6^3 + x1^2*x2*x4*x5^2*x6^2 + x4^4*x5 + x1^2*x2*x3^2*x6 + x1*x3^2*x5^2 + x2*x3*x4*x5^3 + ⋮ + +julia> BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1, 0], operators = [1,2,1]) +Set{ZZMPolyRingElem} with 3 elements: + 1 + x3 + x2*x3 + +julia> BasisLieHighestWeight.basis_lie_highest_weight("C", 3, [1, 1, 1], monomial_order = "Lex") +Set{ZZMPolyRingElem} with 512 elements: + x1*x5*x6*x8 + x6^4 + x3*x4^2*x8 + x3*x4*x6*x7 + x8^3 + x3*x6^2 + x2*x3 + x5*x6^2*x9 + x6*x8^2*x9 + x1*x6*x7 + x5*x6*x9^2 + x6^2*x7^2*x8 + x5*x7*x8 + x4*x6^2*x7*x8^2 + x4^2*x5*x7 + x1*x5^2*x6 + x1*x6*x8 + x3*x4*x5 + x2*x4*x6^2*x7 + x4*x6*x7 + x1*x4*x7*x8^2 + ⋮ +``` +""" +function basis_lie_highest_weight(type::String, rank::Int, highest_weight::Vector{Int}; + operators::Union{String, Vector{Int}} = "regular", + monomial_order::Union{String, Function} = "GRevLex", cache_size::Int = 0, + parallel::Bool = false, return_no_minkowski::Bool = false, + return_operators::Bool = false) + """ + Pseudocode: + + basis_lie_highest_weight(highest_weight) + return compute_monomials(highest_weight) + + compute_monomials(highest_weight) + if highest_weight was already computed + return old results + if highest_weight = [0, ..., 0] or [0, ..., 1, ..., 0] + return add_by_hand(highest_weight, {}) + else + set_mon = {} + go through all partitions lambda_1 + lambda_2 = highest_weight + add compute_monomials(lambda_1) (+) compute_monomials(lambda_1) to set_mon + if set_mon too small + add_by_hand(highest_weight, set_mon) + return set_mon + + add_by_hand(highest_weight, set_mon) + add_known_monomials(set_mon) + go through all weightspaces that are not full + add_new_monomials(weightspace, set_mon) + return set_mon + + add_known_monomials(set_mon) + add all monomials from set_mon to basis + + add_new_monomials(weightspace, set_mon) + calculate monomials with weight in weightspace + go through them one by one in monomial_order until basis is full + return set_mon + """ + # The function precomputes objects that are independent of the highest weight and that can be used in all recursion + # steps. Then it starts the recursion and returns the result. + + # initialization of objects that can be precomputed + lie_algebra, chevalley_basis = create_lie_lgebra(type, rank) # lie_algebra of type, rank and its chevalley_basis + # operators that are represented by our monomials. x_i is connected to operators[i] + operators = get_operators(type, rank, operators, lie_algebra, chevalley_basis) + weights = weights_for_operators(lie_algebra, chevalley_basis[3], operators) # weights of the operators + weights = (weight->Int.(weight)).(weights) + weights_eps = [w_to_eps(type, rank, w) for w in weights] # other root system + ZZx, x = PolynomialRing(ZZ, length(operators)) # for our monomials + monomial_order_lt = get_monomial_order_lt(monomial_order, ZZx, x) # less than function to sort monomials by order + + # save computations from recursions + calc_highest_weight = Dict{Vector{Int}, Set{ZZMPolyRingElem}}([0 for i in 1:rank] => Set([ZZx(1)])) + # we save all highest weights, for which the Minkowski-sum did not suffice to gain all monomials + no_minkowski = Set{Vector{Int}}() + + # start recursion over highest_weight + monomial_basis = compute_monomials(type, rank, lie_algebra, ZZx, x, highest_weight, operators, weights, + weights_eps, monomial_order_lt, calc_highest_weight, cache_size, parallel, no_minkowski) + + # output + if return_no_minkowski && return_operators + return monomial_basis, no_minkowski, operators + elseif return_no_minkowski + return monomial_basis, no_minkowski + elseif return_operators + return monomial_basis, operators + else + return monomial_basis + end +end + +function sub_simple_refl(word::Vector{Int}, lie_algebra::GAP.Obj)::GAP.Obj + """ + substitute simple reflections (i,i+1), saved in dec by i, with E_{i,i+1} + """ + root_system = GAP.Globals.RootSystem(lie_algebra) + canonical_generators = fromGap(GAP.Globals.CanonicalGenerators(root_system)[1], recursive = false) + operators = GAP.Obj([canonical_generators[i] for i in word], recursive = false) + return operators +end + +function get_operators(type::String, rank::Int, operators::Union{String, Vector{Int}}, lie_algebra::GAP.Obj, + chevalley_basis::GAP.Obj)::GAP.Obj + """ + handles user input for operators + "regular" for all operators + "longest-word" for random longest-word in Weyl-group (currently not implemented) + operators::Vector{Int} for explicit longest-word + """ + # create standard operators + if operators == "regular" # use operators as specified by GAP + operators = chevalley_basis[1] + return operators + # The functionality longest-word required Coxetergroups from Gapjm.jl (https://github.com/jmichel7/Gapjm.jl and was + # temporarily deleted + # choose a random longest word. Created by extending by random not leftdescending reflections until total length is + # reached + #elseif operators == "longest-word" + # operators = longest_weyl_word(t,n) + # operators = sub_simple_refl(operators, lie_algebra, n) + # return operators + end + + # use user defined operators + # wrong input + if !(typeof(operators) == Vector{Int}) + println("operators needs to be of type Vector{Int}") + return -1 + end + if !(all([(1 <= i <= rank) for i in operators])) + println("all values of operators need to between 1 and the rank of the lie algebra.") + end + # If one of the conditions is met, the algorithms works. Otherwise a warning is printed (and can be ignored). + #if !(is_longest_weyl_word(type, rank, operators)) && !(Set(operators) == [i for i=1:n]) + # println("WARNING: operators may be incorrect input.") + #end + operators = sub_simple_refl(operators, lie_algebra) + return operators +end + +function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, x::Vector{ZZMPolyRingElem}, + highest_weight::Vector{Int}, operators::GAP.Obj, weights::Vector{Vector{Int64}}, + weights_eps::Vector{Vector{Int64}}, monomial_order_lt::Function, + calc_highest_weight::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, cache_size::Int, + parallel::Bool, no_minkowski::Set{Vector{Int}})::Set{ZZMPolyRingElem} + """ + This function calculates the monomial basis M_{highest_weight} recursively. The recursion saves all computed + results in calc_highest_weight and we first check, if we already encountered this highest weight in a prior step. + If this is not the case, we need to perform computations. The recursion works by using the Minkowski-sum. + If M_{highest_weight} is the desired set of monomials (identified by the exponents as lattice points), it is know + that for lambda_1 + lambda_2 = highest_weight we have M_{lambda_1} + M_{lambda_2} subseteq M_{highest_weight}. + The complexity grows exponentially in the size of highest_weight. Therefore, it is very helpful to obtain a part of + M_{highest_weight} by going through all partitions of highest_weight and using the Minkowski-property. The base + cases of the recursion are the fundamental weights highest_weight = [0, ..., 1, ..., 0]. In this case, or if the + Minkowski-property did not find enough monomials, we need to perform the computations "by hand". + """ + # simple cases + # we already computed the highest_weight result in a prior recursion step + if haskey(calc_highest_weight, highest_weight) + return calc_highest_weight[highest_weight] + elseif highest_weight == [0 for i in 1:rank] # we mathematically know the solution + return Set(ZZx(1)) + end + + # calculation required + # gap_dim is number of monomials that we need to find, i.e. |M_{highest_weight}|. + # if highest_weight is a fundamental weight, partition into smaller summands is possible. This is the basecase of + # the recursion. + gap_dim = GAP.Globals.DimensionOfHighestWeightModule(lie_algebra, GAP.Obj(highest_weight)) # fundamental weights + if is_fundamental(highest_weight) || sum(abs.(highest_weight)) == 0 + push!(no_minkowski, highest_weight) + set_mon = add_by_hand(type, rank, lie_algebra, ZZx, x, highest_weight, operators, weights, weights_eps, + monomial_order_lt, gap_dim, Set{ZZMPolyRingElem}(), cache_size, parallel) + push!(calc_highest_weight, highest_weight => set_mon) + return set_mon + else + # use Minkowski-Sum for recursion + set_mon = Set{ZZMPolyRingElem}() + i = 0 + sub_weights = compute_sub_weights(highest_weight) + l = length(sub_weights) + # go through all partitions lambda_1 + lambda_2 = highest_weight until we have enough monomials or used all + # partitions + while length(set_mon) < gap_dim && i < l + i += 1 + lambda_1 = sub_weights[i] + lambda_2 = highest_weight .- lambda_1 + mon_lambda_1 = compute_monomials(type, rank, lie_algebra, ZZx, x, lambda_1, operators, weights, weights_eps, + monomial_order_lt, calc_highest_weight, cache_size, parallel, + no_minkowski) + mon_lambda_2 = compute_monomials(type, rank, lie_algebra, ZZx, x, lambda_2, operators, weights, weights_eps, + monomial_order_lt, calc_highest_weight, cache_size, parallel, + no_minkowski) + # Minkowski-sum: M_{lambda_1} + M_{lambda_2} \subseteq M_{highest_weight}, if monomials get identified with + # points in ZZ^n + mon_sum = Set([p*q for p in mon_lambda_1 for q in mon_lambda_2]) + union!(set_mon, mon_sum) + end + # check if we found enough monomials + if length(set_mon) < gap_dim + push!(no_minkowski, highest_weight) + set_mon = add_by_hand(type, rank, lie_algebra, ZZx, x, highest_weight, operators, weights, weights_eps, + monomial_order_lt, gap_dim, set_mon, cache_size, parallel) + end + push!(calc_highest_weight, highest_weight => set_mon) + return set_mon + end +end + +@doc """ + is_fundamental(highest_weight::Vector{Int})::Bool + + returns true if ``highest_weight`` is fundamental, i.e. [0, ..., 1, ..., 0] + +# Examples +```jldoctest +julia> BasisLieHighestWeight.is_fundamental([0, 1, 0]) +true + +julia> BasisLieHighestWeight.is_fundamental([0, 1, 1]) +false +``` +""" +function is_fundamental(highest_weight::Vector{Int})::Bool + one = false + for i in highest_weight + if i > 0 + if one || i > 1 + return false + else + one = true + end + end + end + return false +end + +function compute_sub_weights(highest_weight::Vector{Int})::Vector{Vector{Int}} + """ + returns list of weights w != 0 with 0 <= w <= highest_weight elementwise + """ + sub_weights = [] + foreach(Iterators.product((0:x for x in highest_weight)...)) do i + push!(sub_weights, [i...]) + end + popfirst!(sub_weights) # [0, ..., 0] + pop!(sub_weights) # highest_weight + sort!(sub_weights, by=x->sum((x).^2)) + return sub_weights +end + +function add_known_monomials!(weight::Vector{Int}, set_mon_in_weightspace::Dict{Vector{Int64}, + Set{ZZMPolyRingElem}}, number_of_operators::Int, weights::Vector{Vector{Int64}}, + matrices_of_operators::Vector{SMat{ZZRingElem}}, calc_monomials::Dict{ZZMPolyRingElem, + Tuple{TVec, Vector{Int}}}, x::Vector{ZZMPolyRingElem}, + space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}}, + v0::SRow{ZZRingElem}, cache_size::Int) + """ + By using the Minkowski-sum, we know that all monomials in set_mon_in_weightspace are in our basis. Since we want to + extend the weightspace with missing monomials, we need to calculate and add the vector of each monomial to our + basis. + """ + for mon in set_mon_in_weightspace[weight] + # calculate the vector vec associated with mon + if cache_size == 0 + d = sz(matrices_of_operators[1]) + vec = calc_vec(v0, mon, matrices_of_operators) + else + vec = calc_new_mon!(x , mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, + cache_size) + end + + # check if vec extends the basis + if !haskey(space, weight) + space[weight] = nullSpace() + end + add_and_reduce!(space[weight], vec) + end +end + +function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, x::Vector{ZZMPolyRingElem}, + matrices_of_operators::Vector{SMat{ZZRingElem}}, number_of_operators::Int, + weights::Vector{Vector{Int}}, monomial_order_lt::Function, weight::Vector{Int}, + dim_weightspace::Int, weights_eps::Vector{Vector{Int}}, + set_mon_in_weightspace::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, space::Dict{Vector{Int64}, + Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}}, v0::SRow{ZZRingElem}, + cache_size::Int, set_mon::Set{ZZMPolyRingElem}) + """ + If a weightspace is missing monomials, we need to calculate them by trial and error. We would like to go through all + monomials in the order monomial_order_lt and calculate the corresponding vector. If it extends the basis, we add it + to the result and else we try the next one. We know, that all monomials that work lay in the weyl-polytope. + Therefore, we only inspect the monomials that lie both in the weyl-polytope and the weightspace. Since the weyl- + polytope is bounded these are finitely many and we can sort them and then go trough them, until we found enough. + """ + + # get monomials from weyl-polytope that are in the weightspace, sorted by monomial_order_lt + poss_mon_in_weightspace = convert_lattice_points_to_monomials(ZZx, + get_lattice_points_of_weightspace(weights_eps, w_to_eps(type, rank, weight), type)) + poss_mon_in_weightspace = sort(poss_mon_in_weightspace, lt=monomial_order_lt) + + # check which monomials should get added to the basis + i=0 + if weight == 0 # check if [0 0 ... 0] already in basis + i += 1 + end + number_mon_in_weightspace = length(set_mon_in_weightspace[weight]) + # go through possible monomials one by one and check if it extends the basis + while number_mon_in_weightspace < dim_weightspace + i += 1 + + mon = poss_mon_in_weightspace[i] + if mon in set_mon + continue + end + + # calculate the vector vec associated with mon + if cache_size == 0 + d = sz(matrices_of_operators[1]) + vec = calc_vec(v0, mon, matrices_of_operators) + else + vec = calc_new_mon!(x , mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, + cache_size) + end + #println("vec:" , vec) + + # check if vec extends the basis + if !haskey(space, weight) + space[weight] = nullSpace() + end + vec_red = add_and_reduce!(space[weight], vec) + if isempty(vec_red) # v0 == 0 + continue + end + + # save monom + number_mon_in_weightspace += 1 + push!(set_mon, mon) + end +end + + +function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, x::Vector{ZZMPolyRingElem}, + highest_weight::Vector{Int}, operators::GAP.Obj, weights::Vector{Vector{Int64}}, + weights_eps::Vector{Vector{Int64}}, monomial_order_lt::Function, gap_dim::Int, + set_mon::Set{ZZMPolyRingElem}, cache_size::Int, parallel::Bool)::Set{ZZMPolyRingElem} + """ + This function calculates the missing monomials by going through each non full weightspace and adding possible + monomials manually by computing their corresponding vectors and checking if they enlargen the basis. + """ + #println("add_by_hand: ", highest_weight) + # initialization + # matrices g_i for (g_1^a_1 * ... * g_k^a_k)*v + matrices_of_operators = tensorMatricesForOperators(lie_algebra, highest_weight, operators) + number_of_operators = length(matrices_of_operators) + e = [1*(1:number_of_operators .== i) for i in 1:number_of_operators] # e_i + space = Dict(0*weights[1] => nullSpace()) # span of basis vectors to keep track of the basis + v0 = sparse_row(ZZ, [(1,1)]) # starting vector v + # saves the calculated vectors to decrease necessary matrix multiplicatons + calc_monomials = Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}(ZZx(1) => (v0, 0 * weights[1])) + push!(set_mon, ZZx(1)) + # required monomials of each weightspace + weightspaces = get_dim_weightspace(type, rank, lie_algebra, highest_weight) + + # sort the monomials from the minkowski-sum by their weightspaces + set_mon_in_weightspace = Dict{Vector{Int}, Set{ZZMPolyRingElem}}() + for (weight, _) in weightspaces + set_mon_in_weightspace[weight] = Set{ZZMPolyRingElem}() + end + for mon in set_mon + weight = calc_weight(mon, weights) + push!(set_mon_in_weightspace[weight], mon) + end + + # only inspect weightspaces with missing monomials + weights_with_full_weightspace = Set{Vector{Int}}() + for (weight, dim_weightspace) in weightspaces + if (length(set_mon_in_weightspace[weight]) == dim_weightspace) + push!(weights_with_full_weightspace, weight) + end + end + delete!(weightspaces, weights_with_full_weightspace) + + # use parallel computations if parallel = true. The weightspaces could be calculated completely indepent (except for + # the caching). This is not implemented, since I used the package Distributed.jl for this, which is not in the + # Oscar dependencies. But I plan to reimplement this. + # insert known monomials into basis + for (weight, _) in weightspaces + add_known_monomials!(weight, set_mon_in_weightspace, number_of_operators, weights, matrices_of_operators, + calc_monomials, x, space, e, v0, cache_size) + end + + # calculate new monomials + for (weight, dim_weightspace) in weightspaces + add_new_monomials!(type, rank, ZZx, x, matrices_of_operators, number_of_operators, weights, monomial_order_lt, + weight, dim_weightspace, weights_eps, set_mon_in_weightspace, calc_monomials, space, e, v0, + cache_size, set_mon) + end + return set_mon +end + +function get_dim_weightspace(type::String, rank::Int, lie_algebra::GAP.Obj, + highest_weight::Vector{Int})::Dict{Vector{Int}, Int} + """ + Calculates dictionary with weights as keys and dimension of corresponding weightspace as value. GAP computes the + dimension for all positive weights. The dimension is constant on orbits of the weylgroup, and we can therefore + calculate the dimension of each weightspace. + """ + # calculate dimension for dominant weights with GAP + root_system = GAP.Globals.RootSystem(lie_algebra) + dominant_weights, dominant_weights_dim = fromGap(GAP.Globals.DominantCharacter(root_system, + GAP.Obj(highest_weight))) + dominant_weights = convert(Vector{Vector{Int}}, dominant_weights) + weightspaces = Dict{Vector{Int}, Int}() + + # calculate dimension for the rest by checking which positive weights lies in the orbit. + for i in 1:length(dominant_weights) + orbit_weights = orbit_weylgroup(type, rank, lie_algebra, dominant_weights[i]) + dim_weightspace = dominant_weights_dim[i] + for weight in orbit_weights + weightspaces[highest_weight - weight] = dim_weightspace + end + end + return weightspaces +end + +end +export BasisLieHighestWeight diff --git a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl new file mode 100644 index 000000000000..d5255631a525 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl @@ -0,0 +1,54 @@ +using Oscar + +fromGap = Oscar.GAP.gap_to_julia + + +function create_lie_lgebra(type::String, rank::Int)::Tuple{GAP.Obj, GAP.Obj} + """ + Creates the Lie-algebra as a GAP object that gets used for a lot of other computations with GAP + """ + lie_algebra = GAP.Globals.SimpleLieAlgebra(GAP.Obj(type), rank, GAP.Globals.Rationals) + return lie_algebra, GAP.Globals.ChevalleyBasis(lie_algebra) +end + + +gapReshape(A) = sparse_matrix(QQ, hcat(A...)) + +# temporary workaround for issue 2128 +function multiply_scalar(A::SMat{T}, d) where T + for i in 1:nrows(A) + scale_row!(A, i, T(d)) + end + return A +end + +function matricesForOperators(lie_algebra::GAP.Obj, highest_weight::Vector{Int}, + ops::GAP.Obj)::Vector{SMat{ZZRingElem}} + """ + used to create tensorMatricesForOperators + """ + M = Oscar.GAP.Globals.HighestWeightModule(lie_algebra, Oscar.GAP.julia_to_gap(highest_weight)) + matrices_of_operators = Oscar.GAP.Globals.List(ops, o -> Oscar.GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o)) + matrices_of_operators = gapReshape.( Oscar.GAP.gap_to_julia(matrices_of_operators)) + denominators = map(y->denominator(y[2]), union(union(matrices_of_operators...)...)) + common_denominator = lcm(denominators)# // 1 + matrices_of_operators = (A->change_base_ring(ZZ, multiply_scalar(A, common_denominator))).(matrices_of_operators) + return matrices_of_operators +end + + +function weights_for_operators(lie_algebra::GAP.Obj, cartan::GAP.Obj, operators::GAP.Obj)::Vector{Vector{Int}} + """ + Calculates the weight wts[i] for each operator ops[i] + """ + cartan = fromGap(cartan, recursive=false) + operators = fromGap(operators, recursive=false) + asVec(v) = fromGap(GAP.Globals.ExtRepOfObj(v)) + if any(iszero.(asVec.(operators))) + error("ops should be non-zero") + end + nzi(v) = findfirst(asVec(v) .!= 0) + return [ + [asVec(h*v)[nzi(v)] / asVec(v)[nzi(v)] for h in cartan] for v in operators + ] +end diff --git a/experimental/BasisLieHighestWeight/src/MonomialOrder.jl b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl new file mode 100644 index 000000000000..5185fd8f3ebf --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl @@ -0,0 +1,18 @@ +function get_monomial_order_lt(monomial_order::Union{String, Function}, ZZx::ZZMPolyRing, + x::Vector{ZZMPolyRingElem})::Function + """ + Returns the desired monomial_order function less than + """ + #if isa(monomial_order, Function) + # choosen_monomial_order = monomial_order + if monomial_order == "GRevLex" + choosen_monomial_order = degrevlex(x) + elseif monomial_order == "RevLex" + choosen_monomial_order = revlex(x) + elseif monomial_order == "Lex" + choosen_monomial_order = lex(x) + else + println("no suitable order picked") + end + return (mon1, mon2) -> (cmp(choosen_monomial_order, mon1, mon2) == -1) +end diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl new file mode 100644 index 000000000000..cba90354ce9f --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -0,0 +1,101 @@ +include("./VectorSpaceBases.jl") +include("./TensorModels.jl") +include("./LieAlgebras.jl") +include("./MonomialOrder.jl") +include("./WeylPolytope.jl") + +fromGap = Oscar.GAP.gap_to_julia + + +function calc_weight(mon::ZZMPolyRingElem, weights::Vector{Vector{Int}})::Vector{Int} + """ + calculates weight associated with monomial mon + """ + degree_mon = degrees(mon) + weight = [0 for i in 1:length(weights[1])] + for i in 1:length(degree_mon) + weight .+= degree_mon[i] * weights[i] + end + return weight +end + +function calc_vec(v0::SRow{ZZRingElem}, mon::ZZMPolyRingElem, + matrices_of_operators::Vector{SMat{ZZRingElem}})::SRow{ZZRingElem} + """ + calculates vector associated with monomial mon + """ + vec = v0 + degree_mon = degrees(mon) + for i in length(degree_mon):-1:1 + for j in 1:degree_mon[i] + vec = mul(vec, transpose(matrices_of_operators[i])) # currently there is no sparse matrix * vector mult + end + end + return vec +end + +function highest_calc_sub_monomial(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + number_of_operators::Int)::ZZMPolyRingElem + """ + returns the key in calc_monomials that can be extended by the least amount of left-operations to mon + """ + sub_mon = copy(mon) + + for i in 1:number_of_operators + while is_divisible_by(sub_mon, x[i]) + if haskey(calc_monomials, sub_mon) + return sub_mon + else + sub_mon /= x[i] + end + end + end + return sub_mon +end + +function calc_new_mon!(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, weights::Vector{Vector{Int}}, + matrices_of_operators::Vector{SMat{ZZRingElem}}, number_of_operators::Int, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}}, + cache_size::Int)::SRow{ZZRingElem} + # calculate vector of mon by extending a previous calculated vector to a + # monom that differs only by left-multiplication, save results in calc_monomials + sub_mon = highest_calc_sub_monomial(x, mon, calc_monomials, number_of_operators) + #println("sub_mon: ", sub_mon) + sub_mon_cur = copy(sub_mon) + (vec, weight) = calc_monomials[sub_mon] + for i in number_of_operators:-1:1 + for k in degrees(sub_mon)[i]:(degrees(mon)[i]-1) + sub_mon_cur *= x[i] + weight += weights[i] + if !haskey(space, weight) + space[weight] = nullSpace() + end + + vec = mul(vec, transpose(matrices_of_operators[i])) # currently there is no sparse matrix * vector mult + if length(calc_monomials) < cache_size + calc_monomials[sub_mon_cur] = (vec, weight) + end + + # check if the extended monomial can be deleted from calculated_monomials, i.e. the other possible + # extensions by left multiplication with some x[i] are already contained + can_be_deleted = true + k = number_of_operators + for l in 1:number_of_operators + if degrees(sub_mon_cur - x[i])[l] != 0 + k = l + end + end + for l in 1:k + can_be_deleted = can_be_deleted && haskey(calc_monomials, sub_mon_cur - x[i] + x[l]) + end + if can_be_deleted && sub_mon_cur != x[i] + delete!(calc_monomials, sub_mon_cur - x[i]) + end + end + end + #println(length(calc_monomials)) + return vec +end + diff --git a/experimental/BasisLieHighestWeight/src/RootConversion.jl b/experimental/BasisLieHighestWeight/src/RootConversion.jl new file mode 100644 index 000000000000..6e1191b2d89f --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/RootConversion.jl @@ -0,0 +1,366 @@ +#using Gapjm +using Oscar + +fromGap = Oscar.GAP.gap_to_julia + +function w_to_eps(type::String, rank::Int, weight::Vector{Int})::Vector{Int} + """ + converts weight in rootsystem w_i to eps_i + """ + if type in ["A", "B", "C", "D", "E", "F", "G"] + return alpha_to_eps(type, rank, w_to_alpha(type, rank, weight)) + else + println("Type needs to be one of A-D") + end +end + +function eps_to_w(type::String, rank::Int, weight::Vector{Int})::Vector{Int} + """ + converts weight in rootsystem eps_i to w_i + """ + if type in ["A", "B", "C", "D", "E", "F", "G"] + return round.(alpha_to_w(type, rank, eps_to_alpha(type, rank, weight))) + else + println("Type needs to be one of A-D") + end +end + +function alpha_to_eps(type::String, rank::Int, weight::Vector{Int})::Vector{Int} + """ + converts weight in rootsystem alpha_i to eps_i + """ + if type == "A" + return alpha_to_eps_A(rank, weight) + elseif type in ["B", "C", "D"] + return alpha_to_eps_BCD(type, rank, weight) + elseif type == "E" && rank in [6, 7, 8] + return alpha_to_eps_E(rank, weight) + elseif type == "F" && rank == 4 + return alpha_to_eps_F(weight) + elseif type == "G" && rank == 2 + return alpha_to_eps_G(weight) + else + println("This rank of lie algebra is not supported.") + end +end + +function eps_to_alpha(type::String, rank::Int, weight::Vector{Int})::Vector{Int} + """ + converts weight in rootsystem eps_i to alpha_i + """ + if type == "A" + return eps_to_alpha_A(rank, weight) + elseif type in ["B", "C", "D"] + return eps_to_alpha_BCD(type, rank, weight) + elseif type == "E" && rank in [6, 7, 8] + return eps_to_alpha_E(rank, weight) + elseif type == "F" && rank == 4 + return eps_to_alpha_F(weight) + elseif type == "G" && rank == 2 + return eps_to_alpha_G(weight) + else + println("This rank of lie algebra is not supported.") + end +end + +function w_to_alpha(type, rank, weight::Vector{Int})::Vector{Int} + C = get_CartanMatrix(type, rank) + return [i for i in C*weight] +end + +function alpha_to_w(type::String, rank::Int, weight::Vector{Int})::Vector{Int} + C_inv = get_inverse_CartanMatrix(type, rank) + return [i for i in C_inv*weight] +end + +function get_CartanMatrix(type::String, rank::Int) + L = GAP.Globals.SimpleLieAlgebra(GAP.Obj(type), rank, GAP.Globals.Rationals) + R = GAP.Globals.RootSystem(L) + C_list = fromGap(GAP.Globals.CartanMatrix(R)) + C = zeros(rank, rank) + for i in 1:rank + for j in 1:rank + C[i,j] = C_list[i][j] + end + end + return C +end + +function get_inverse_CartanMatrix(type::String, rank::Int) + return inv(get_CartanMatrix(type, rank)) +end + +function alpha_to_eps_BCD(type::String, rank::Int, weight::Vector{Int})::Vector{Int} + """ + for B-D + """ + eps = [0.0 for i in 1:rank] + for i in 1:(rank-1) + eps[i] += weight[i] + eps[i+1] -= weight[i] + end + if type == "B" + eps[rank] += weight[rank] + elseif type == "C" + eps[rank] += 2*weight[rank] + elseif type == "D" + eps[rank - 1] += weight[rank] + eps[rank] += weight[rank] + end + return eps +end + +function eps_to_alpha_BCD(type::String, rank::Int, weight::Vector{Int})::Vector{Int} + """ + for B-D + """ + alpha = [0.0 for i in 1:rank] + for i in 1:(rank-2) + alpha[i] = weight[i] + weight[i+1] += weight[i] + end + if type == "B" + alpha[rank - 1] = weight[rank - 1] + alpha[rank] += weight[rank-1] + weight[rank] + elseif type == "C" + alpha[rank - 1] = weight[rank - 1] + alpha[rank] += 0.5*weight[rank - 1] + 0.5*weight[rank] # requires eps to be always even + elseif type == "D" + alpha[rank - 1] += (weight[rank - 1] - weight[rank])/2 + alpha[rank] += (weight[rank - 1] + weight[rank])/2 + end + return alpha +end + +function alpha_to_eps_E(rank::Int, weight::Vector{Int})::Vector{Int} + """ + for E + """ + if rank == 6 + return alpha_to_eps_E6(weight) + elseif rank == 7 + return alpha_to_eps_E7(weight) + elseif rank == 8 + return alpha_to_eps_E8(weight) + end +end + +function eps_to_alpha_E(rank::Int, weight) + """ + for E + """ + if rank == 6 + return eps_to_alpha_E6(weight) + elseif rank == 7 + return eps_to_alpha_E7(weight) + elseif rank == 8 + return eps_to_alpha_E8(weight) + end +end + +function alpha_to_eps_E6(weight::Vector{Int})::Vector{Int} + """ + for E6, potentially wrong order or roots (1-2-3-5-6, 3-4) + """ + eps = [0.0 for i in 1:6] + for i in 1:4 + eps[i] += weight[i] + eps[i + 1] += - weight[i] + end + eps[4] += weight[5] + eps[5] += weight[5] + for i in 1:5 + eps[i] += -0.5*weight[6] + end + eps[6] += 0.5*sqrt(3)*weight[6] + return eps +end + +function eps_to_alpha_E6(weight) + """ + for E6 + """ + alpha = [0.0 for i in 1:6] + for j in 1:3 + for i in 1:j + alpha[j] += weight[i] + end + alpha[j] += j*(sqrt(3) / 3) *weight[6] + end + for i in 1:4 + alpha[4] += 0.5*weight[i] + alpha[5] += 0.5*weight[i] + end + alpha[4] += -0.5*weight[5] + (sqrt(3) / 2)*weight[6] + alpha[5] += 0.5*weight[5] + 5*(sqrt(3) / 6)*weight[6] + alpha[6] = +2*(sqrt(3) / 3)*weight[6] + #println("eps_to_alpha_E6: ", alpha) + return alpha +end + +function alpha_to_eps_E7(weight::Vector{Int})::Vector{Int} + """ + for E7, potentially wrong order of roots (1-2-3-4-6-7, 4-5) + """ + eps = [0.0 for i in 1:7] + for i in 1:5 + eps[i] += weight[i] + eps[i + 1] += - weight[i] + end + eps[5] += weight[6] + eps[6] += weight[6] + for i in 1:6 + eps[i] += -0.5*weight[7] + end + eps[7] += 0.5*sqrt(2)*weight[7] + return eps +end + +function eps_to_alpha_E7(weight::Vector{Int})::Vector{Int} + """ + for E7 + """ + alpha = [0.0 for i in 1:7] + for j in 1:4 + for i in 1:j + alpha[j] += weight[i] + end + alpha[j] += j*(sqrt(2) / 2) *weight[7] + end + for i in 1:5 + alpha[5] += 0.5*weight[i] + alpha[6] += 0.5*weight[i] + end + alpha[5] += -0.5*weight[6] + sqrt(2)*weight[7] + alpha[6] += 0.5*weight[6] + 3*(sqrt(2) / 2)*weight[7] + alpha[7] = sqrt(2)*weight[7] + return alpha +end + +function alpha_to_eps_E8(weight::Vector{Int})::Vector{Int} + """ + for E8 + """ + eps = [0.0 for i in 1:8] + for i in 1:6 + eps[i] += weight[i] + eps[i+1] += - weight[i] + end + eps[6] += weight[7] + eps[7] += weight[7] + for i in 1:8 + eps[i] += -0.5*weight[8] + end + return eps +end + +function eps_to_alpha_E8(weight::Vector{Int})::Vector{Int} + """ + for E8 + """ + alpha = [0.0 for i in 1:8] + for j in 1:5 + for i in 1:j + alpha[j] += weight[i] + end + alpha[j] += -j*weight[8] + end + for i in 1:6 + alpha[6] += 0.5*weight[i] + alpha[7] += 0.5*weight[i] + end + alpha[6] += -0.5*weight[7] - 2.5*weight[8] + alpha[7] += 0.5*weight[7] - 3.5*weight[8] + alpha[8] = -2*weight[8] + return alpha +end + +function alpha_to_eps_F(weight::Vector{Int})::Vector{Int} + """ + for F + """ + eps = [0.0 for i in 1:4] + eps[1] = weight[1] - 0.5*weight[4] + eps[2] = - weight[1] + weight[2] - 0.5*weight[4] + eps[3] = - weight[2] + weight[3] - 0.5*weight[4] + eps[4] = - 0.5*weight[4] + return eps +end + +function eps_to_alpha_F(weight::Vector{Int})::Vector{Int} + """ + for F + """ + alpha = [0 for i in 1:4] + alpha[1] = weight[1] - weight[4] + alpha[2] = weight[1] + weight[2] - 2*weight[4] + alpha[3] = weight[1] + weight[2] + weight[3] - 3*weight[4] + alpha[4] = -2*weight[4] + return alpha +end + +function alpha_to_eps_G(weight::Vector{Int})::Vector{Int} + """ + for G_2 + """ + eps = [0.0 for i in 1:3] + eps[1] = weight[1] - weight[2] + eps[2] = - weight[1] + 2*weight[2] + eps[3] = - weight[2] + choose_representant_eps(eps) + return eps +end + +function eps_to_alpha_G(weight::Vector{Int})::Vector{Int} + """ + for G_2 + """ + alpha = [0.0 for i in 1:2] + if length(weight) >= 3 + weight .-= weight[3] + end + alpha[1] = weight[1] + alpha[2] = (weight[1] + weight[2]) / 3 + return alpha +end + +function choose_representant_eps(weight::Vector{Int}) + # choose representant eps_1 + ... + eps_m = 0 + if any(<(0), weight) # non negative + weight .-= min(weight ...) + end +end + +function alpha_to_eps_A(rank::Int, weight::Vector{Int})::Vector{Int} + """ + for A + """ + eps = [0 for i in 1:(rank + 1)] + for i in 1:rank + eps[i] += weight[i] + eps[i + 1] -= weight[i] + end + choose_representant_eps(eps) + return eps +end + +function eps_to_alpha_A(rank::Int, weight::Vector{Int})::Vector{Int} + """ + for A + """ + if length(weight) == rank + append!(weight, 0) + end + alpha = [0.0 for i in 1:(rank + 1)] + for i in 1:(rank + 1) + for j in 1:i + alpha[i] += weight[j] + end + end + m = alpha[rank + 1] / (rank + 1) + for i in 1:rank + alpha[i] -= i*m + end + pop!(alpha) + return alpha +end diff --git a/experimental/BasisLieHighestWeight/src/TensorModels.jl b/experimental/BasisLieHighestWeight/src/TensorModels.jl new file mode 100644 index 000000000000..0b56cc2635f9 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/TensorModels.jl @@ -0,0 +1,58 @@ +using Oscar + +function kron(A::SMat{ZZRingElem}, B::SMat{ZZRingElem})::SMat{ZZRingElem} + """ + Computes the Kronecker-product of A and B + """ + res = sparse_matrix(ZZ, nrows(A)*nrows(B), ncols(A)*ncols(B)) + for i in 1:nrows(B) + for j in 1:nrows(A) + new_row_tuples = Vector{Tuple{Int, ZZRingElem}}([(1,ZZ(0))]) + for (index_A, element_A) in union(getindex(A, j)) + for (index_B, element_B) in union(getindex(B, i)) + push!(new_row_tuples, ((index_A-1)*ncols(B)+index_B, element_A*element_B)) + end + end + new_row = sparse_row(ZZ, new_row_tuples) + setindex!(res, new_row, (j-1)*nrows(B)+i) + end + end + return res +end + +# temprary fix sparse in Oscar does not work +function tensorProduct(A::SMat{ZZRingElem}, B::SMat{ZZRingElem})::SMat{ZZRingElem} + temp_mat = kron(A, spid(sz(B))) + kron(spid(sz(A)), B) + res = sparse_matrix(ZZ, nrows(A)*nrows(B), ncols(A)*ncols(B)) + for i in 1:nrows(temp_mat) + setindex!(res, getindex(temp_mat, i), i) + end + return res +end + + +spid(n::Int) = identity_matrix(SMat, ZZ, n)::SMat{ZZRingElem} +sz(A::SMat{ZZRingElem}) = nrows(A)::Int #size(A)[1] +#tensorProduct(A, B) = kron(A, spid(sz(B))) + kron(spid(sz(A)), B) +tensorProducts(As, Bs) = (AB->tensorProduct(AB[1], AB[2])).(zip(As, Bs)) +tensorPower(A, n) = (n == 1) ? A : tensorProduct(tensorPower(A, n-1), A) +tensorPowers(As, n) = (A->tensorPower(A, n)).(As) + +function tensorMatricesForOperators(lie_algebra::GAP.Obj, highest_weight::Vector{Int}, + operators::GAP.Obj)::Vector{SMat{ZZRingElem}} + """ + Calculates the matrices g_i corresponding to the operator ops[i]. + """ + matrices_of_operators = [] + for i in 1:length(highest_weight) + if highest_weight[i] <= 0 + continue + end + wi = Int.(1:length(highest_weight) .== i) # i-th fundamental weight + _matrices_of_operators = matricesForOperators(lie_algebra, wi, operators) + _matrices_of_operators = tensorPowers(_matrices_of_operators, highest_weight[i]) + matrices_of_operators = matrices_of_operators == [] ? _matrices_of_operators : + tensorProducts(matrices_of_operators, _matrices_of_operators) + end + return matrices_of_operators +end diff --git a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl new file mode 100644 index 000000000000..f4cf0694a3a5 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl @@ -0,0 +1,75 @@ +# manages the linear (in-)dependence of integer vectors +# this file is only of use to basis_lie_highest_weight for the basis vectors + +using Oscar + +TVec = SRow{ZZRingElem} # TVec is datatype of basisvectors in basisLieHighestWeight +Short = UInt8 # for exponents of monomials; max. 255 + +struct VSBasis + basis_vectors::Vector{TVec} # vector of basisvectors + pivot::Vector{Int} # vector of pivotelements, i.e. pivot[i] is first nonzero element of basis_vectors[i] +end + +nullSpace() = VSBasis([], []) # empty Vektorraum + +reduce_col(a::TVec, b::TVec, i::Int) = (b[i]*a - a[i]*b)::TVec # create zero entry in a + +function normalize(v::TVec)::Tuple{TVec, Int64} + """ + divides vector by gcd of nonzero entries, returns vector and first nonzero index + used: add_and_reduce! + """ + if is_empty(v) + return (v, 0) + end + pivot = first(v)[1] # first nonzero element of vector + return divexact(v, gcd(map(y->y[2], union(v)))), pivot +end + +function add_and_reduce!(sp::VSBasis, v::TVec)::TVec + """ + for each pivot of sp.basis_vectors we make entry of v zero and return the result and insert it into sp + 0 => linear dependent + * => linear independent, new column element of sp.basis_vectors since it increases basis + invariants: the row of a pivotelement in any column in basis_vectors is 0 (except the pivotelement) + elements of basis_vectors are integers, gcd of each column is 1 + """ + # initialize objects + basis_vectors = sp.basis_vectors + pivot = sp.pivot + v, newPivot = normalize(v) + + # case v zero vector + if newPivot == 0 + return v + end + + # use pivots of basis basis_vectors to create zeros in v + for j in 1:length(basis_vectors) + i = pivot[j] + if i != newPivot + continue + end + v = reduce_col(v, basis_vectors[j], i) + v, newPivot = normalize(v) + if newPivot == 0 + #return 0 + return v + end + end + + # new pivot element of v + pos = findfirst(pivot .> newPivot) + if (pos === nothing) + pos = length(pivot) + 1 + end + + # save result + insert!(basis_vectors, pos, v) + insert!(pivot, pos, newPivot) + return v +end + + + diff --git a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl new file mode 100644 index 000000000000..d830702b0228 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl @@ -0,0 +1,158 @@ +#using Gapjm +using Oscar + +include("./RootConversion.jl") + +fromGap = Oscar.GAP.gap_to_julia + +function weylpolytope(type::String, rank::Int, lie_algebra::GAP.Obj, + highest_weight::Vector{Int})::Polymake.BigObjectAllocated + """ + returns weyl-polytope in homogeneous coordinates, i.e. convex hull of orbit of weyl-group of + type type and rank rank on highest weight vector highest_weight + """ + vertices = orbit_weylgroup(type, rank, lie_algebra, highest_weight) + vertices = transpose(hcat(vertices ...)) + vertices = [w for w in eachrow(vertices)] + vertices = transpose(hcat(vertices ...)) + vertices_hom = [ones(Int64, size(vertices)[1]) vertices] # in homogeneous coordinates + weylpoly = polytope.Polytope(POINTS=vertices_hom) + return weylpoly +end + +function orbit_weylgroup(type::String, rank::Int, lie_algebra::GAP.Obj, weight_vector::Vector{Int}) + """ + operates weyl-group of type type and rank rank on vector weight_vector and returns list of vectors in orbit + input and output weights in terms of w_i + """ + # initialization + weyl_group = GAP.Globals.WeylGroup(GAP.Globals.RootSystem(lie_algebra)) + orbit_iterator = GAP.Globals.WeylOrbitIterator(weyl_group, GAP.Obj(weight_vector)) + vertices = [] + + # operate with the weylgroup on weight_vector + GAP.Globals.IsDoneIterator(orbit_iterator) + while !(GAP.Globals.IsDoneIterator(orbit_iterator)) + w = GAP.Globals.NextIterator(orbit_iterator) + push!(vertices, fromGap(w)) + end + + # return + vertices = convert(Vector{Vector{Int}}, vertices) + return vertices +end + +function get_points_polytope(polytope) + """ + returns all points (interior and vertices) of a polytope in regular (i.e. not homogenoues coordinates). + """ + interior_points = convert(Matrix{Int64}, polytope.INTERIOR_LATTICE_POINTS) + vertices_points = convert(Matrix{Int64}, polytope.VERTICES) + points = [interior_points; vertices_points][:, 2:end] + return points +end + + +function convert_lattice_points_to_monomials(ZZx, lattice_points_weightspace) + return [finish(push_term!(MPolyBuildCtx(ZZx), ZZ(1), convert(Vector{Int}, convert(Vector{Int64}, lattice_point)))) + for lattice_point in lattice_points_weightspace] +end + +function get_lattice_points_of_weightspace(weights, weight, type) + """ + calculates all lattice points in a given weightspace for a lie algebra of type type + input: + weights: the operator weights in eps_i + weight: lambda - mu + + output: all lattice points with weight weight + """ + if type in ["A", "G"] + return get_lattice_points_of_weightspace_A_G_n(weights, weight) + else + return get_lattice_points_of_weightspace_Xn(weights, weight) + end +end + +function get_lattice_points_of_weightspace_A_G_n(weights, weight) + """ + calculates all monomials in a given weightspace for lie algebras that have type A or G + input: + weights: the operator weights in eps_i + weight: lambda - mu + + output: all monomials with weight weight + + works by calculating all integer solutions to the following linear program: + [ 1 | | ] [ x ] + [ 1 weights[1]... weights[k]] * [ | ] = weight + [... | | ] [ res ] + [ 1 | | ] [ | ] + where res[i] >= 0 for all i + + example: + weights = [[1, 0, 2], [-1, 1, 1], [0, -1, 0]] (i.e. a_1 = eps_1 - eps_2, a_2 = eps_2 - eps_3, a_12 = eps_1 - eps_3) + weight = [2, 1, 0] + -> poly = polytope.Polytope(INEQUALITIES=[0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1], + EQUATIONS=[-2 1 1 0 2; -1 1 -1 1 1; 0 1 0 -1 0]) + => returns [[1 0 0], [1 1 0]] + """ + # build linear (in-)equalities + weights = [reshape(w, 1, :) for w in weights] + n = length(weights) + ineq = zeros(Int64, n, n+2) + for i in 1:n + ineq[i, 2+i] = 1 + end + equ = cat([-i for i in vec(weight)], [1 for i=1:length(weight)], dims = (2,2)) + equ = cat(equ, [transpose(w) for w in weights] ..., dims = (2,2)) + + # find integer solutions of linear (in-)equation as lattice points of polytope + poly = polytope.Polytope(INEQUALITIES=ineq, EQUATIONS=equ) + + # convert lattice-points to Oscar monomials + lattice_points_weightspace = lattice_points(Polyhedron(poly)) + lattice_points_weightspace = [lattice_point[2:end] for lattice_point in lattice_points_weightspace] + return lattice_points_weightspace + +end + +function get_lattice_points_of_weightspace_Xn(weights, weight) + """ + calculates all lattice points in a given weightspace for lie algebras that don't have type A + input: + weights: the operator weights in eps_i + weight: lambda - mu + + output: all lattice points with weight weight + + works by calculating all integer solutions to the following linear program: + [ | | ] [ x ] + [weights[1]...weights[k]] * [ | ] = weight + [ | | ] [ res ] + [ | | ] [ | ] + where res[i] >= 0 for all i + + example: + weights = [[1, 0, 2], [-1, 1, 1], [0, -1, 0]] (i.e. a_1 = eps_1 - eps_2, a_2 = eps_2 - eps_3, a_12 = eps_1 - eps_3) + weight = [2, 1, 0] + -> poly = polytope.Polytope(INEQUALITIES=[0 1 0 0; 0 0 1 0; 0 0 0 1], EQUATIONS=[-2 1 0 2; -1 -1 1 1; 0 0 -1 0]) + => returns + """ + # build linear (in-)equalities + weights = [reshape(w, 1, :) for w in weights] + n = length(weights) + ineq = zeros(Int64, n, n+1) + for i in 1:n + ineq[i, 1+i] = 1 + end + equ = [-i for i in vec(weight)] + equ = cat(equ, [transpose(w) for w in weights] ..., dims = (2,2)) + + # find integer solutions of linear (in-)equation as lattice points of polytope + poly = polytope.Polytope(INEQUALITIES=ineq, EQUATIONS=equ) + + # convert lattice-points to Oscar monomials + lattice_points_weightspace = lattice_points(Polyhedron(poly)) + return lattice_points_weightspace +end diff --git a/experimental/BasisLieHighestWeight/test/MBOld.jl b/experimental/BasisLieHighestWeight/test/MBOld.jl new file mode 100644 index 000000000000..b990e3d39127 --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/MBOld.jl @@ -0,0 +1,301 @@ + # This file is an old version of the algorithm that can compute (not all cases) of + # BasisLieHighestWeight.basis_lie_highest_weight and is used only in runtests.jl to check that the newer algorithm matches + # There is code doubling, but I am not sure how the src part is going to change when its integrated with the other + # lie algebra work and would prefer to change this file after doing the integration to make sure that everything stays + # correct. + + +module MBOld + +export basisLieHighestWeight + +using Oscar + + +G = Oscar.GAP.Globals +forGap = Oscar.GAP.julia_to_gap +fromGap = Oscar.GAP.gap_to_julia + +TVec = SRow{ZZRingElem} +Short = UInt8 + +struct VSBasis + A::Vector{TVec} + pivot::Vector{Int} +end + + +nullSpace() = VSBasis([], []) + + +function normalize(v::TVec) + """ + divides vector by gcd of nonzero entries, returns vector and first nonzero index + used: addAndReduce! + """ + if isempty(v) + return (0, 0) + end + + pivot = first(v)[1] + + return divexact(v, gcd(map(y->y[2], union(v)))), pivot +end + + +reduceCol(a, b, i::Int) = b[i]*a - a[i]*b + + +function addAndReduce!(sp::VSBasis, v::TVec) + """ + for each pivot of sp.A we make entry of v zero and return the result + 0 => linear dependent + * => linear independent, new column element of sp.A since it increases basis + invariants: the row of a pivotelement in any column in A is 0 (except the pivotelement) + elements of A are integers, gcd of each column is 1 + """ + A = sp.A + pivot = sp.pivot + v, newPivot = normalize(v) + if newPivot == 0 + return 0 + end + + for j = 1:length(A) + i = pivot[j] + if i != newPivot + continue + end + v = reduceCol(v, A[j], i) + v, newPivot = normalize(v) + if newPivot == 0 + return 0 + end + end + + pos = findfirst(pivot .> newPivot) + if (pos === nothing) + pos = length(pivot) + 1 + end + + insert!(A, pos, v) + insert!(pivot, pos, newPivot) + + return v +end + + +#### Lie algebras + +G = Oscar.GAP.Globals +forGap = Oscar.GAP.julia_to_gap +fromGap = Oscar.GAP.gap_to_julia + + +function lieAlgebra(t::String, n::Int) + L = G.SimpleLieAlgebra(forGap(t), n, G.Rationals) + return L, G.ChevalleyBasis(L) +end + +# temporary workaround for issue 2128 +function multiply_scalar(A::SMat{T}, d) where T + for i in 1:nrows(A) + scale_row!(A, i, T(d)) + end + return A + #return identity_matrix(SMat, QQ, size(A)[1])*A +end + +gapReshape(A) = sparse_matrix(QQ, hcat(A...)) + +function matricesForOperators(L, hw, ops) + """ + used to create tensorMatricesForOperators + """ + M = G.HighestWeightModule(L, forGap(hw)) + mats = G.List(ops, o -> G.MatrixOfAction(G.Basis(M), o)) + mats = gapReshape.(fromGap(mats)) + denominators = map(y->denominator(y[2]), union(union(mats...)...)) + #d = convert(QQ, lcm(denominators)) + d = lcm(denominators)# // 1 + mats = (A->change_base_ring(ZZ, multiply_scalar(A, d))).(mats) + return mats +end + +function weightsForOperators(L, cartan, ops) + cartan = fromGap(cartan, recursive=false) + ops = fromGap(ops, recursive=false) + asVec(v) = fromGap(G.ExtRepOfObj(v)) + if any(iszero.(asVec.(ops))) + error("ops should be non-zero") + end + nzi(v) = findfirst(asVec(v) .!= 0) + return [ + [asVec(h*v)[nzi(v)] / asVec(v)[nzi(v)] for h in cartan] for v in ops + ] +end + +#### tensor model + +function kron(A, B) + res = sparse_matrix(ZZ, nrows(A)*nrows(B), ncols(A)*ncols(B)) + for i in 1:nrows(B) + for j in 1:nrows(A) + new_row_tuples = Vector{Tuple{Int, ZZRingElem}}([(1,ZZ(0))]) + for (index_A, element_A) in union(getindex(A, j)) + for (index_B, element_B) in union(getindex(B, i)) + push!(new_row_tuples, ((index_A-1)*ncols(B)+index_B, element_A*element_B)) + end + end + new_row = sparse_row(ZZ, new_row_tuples) + setindex!(res, new_row, (j-1)*nrows(B)+i) + end + end + return res +end + +# temprary fix sparse in Oscar does not work +function tensorProduct(A, B) + temp_mat = kron(A, spid(sz(B))) + kron(spid(sz(A)), B) + res = sparse_matrix(ZZ, nrows(A)*nrows(B), ncols(A)*ncols(B)) + for i in 1:nrows(temp_mat) + setindex!(res, getindex(temp_mat, i), i) + end + return res +end + +spid(n) = identity_matrix(SMat, ZZ, n) +sz(A) = nrows(A) +tensorProducts(As, Bs) = (AB->tensorProduct(AB[1], AB[2])).(zip(As, Bs)) +tensorPower(A, n) = (n == 1) ? A : tensorProduct(tensorPower(A, n-1), A) +tensorPowers(As, n) = (A->tensorPower(A, n)).(As) + + +function tensorMatricesForOperators(L, hw, ops) + """ + Calculates the matrices g_i corresponding to the operator ops[i]. + """ + #println("hw: ", hw) + mats = [] + + for i in 1:length(hw) + if hw[i] <= 0 + continue + end + wi = Int.(1:length(hw) .== i) # i-th fundamental weight + _mats = matricesForOperators(L, wi, ops) + _mats = tensorPowers(_mats, hw[i]) + if size(mats)[1] > 0 + A = _mats[1] + B = mats[1] + end + mats = mats == [] ? _mats : tensorProducts(mats, _mats) + end + return mats +end + + +""" + basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; parallel::Bool = true) :: Tuple{Vector{Vector{Short}},Vector{TVec}} + +Compute a monomial basis for the highest weight module with highest weight ``hw`` (in terms of the fundamental weights), for a simple Lie algebra of type ``t`` and rank ``n``. + +Example +====== +```jldoctest +julia> dim, monomials, vectors = PolyBases.MB.basisLieHighestWeight("A", 2, [1,0]) +(3, Vector{UInt8}[[0x00, 0x00, 0x00], [0x01, 0x00, 0x00], [0x00, 0x00, 0x01]], SparseArrays.SparseVector{Int64, Int64}[ [1] = 1, [2] = 1, [3] = 1]) +``` +""" + +function basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; roots = []) #--- :: Tuple{Int64,Vector{Vector{Short}},Vector{TVec}} + L, CH = lieAlgebra(t, n) + ops = CH[1] # positive root vectors + # .. reorder.. + wts = weightsForOperators(L, CH[3], ops) + wts = (v->Int.(v)).(wts) + + mats = tensorMatricesForOperators(L, hw, ops) + hwv = sparse_row(ZZ, [(1,1)]) + + monomials = compute(hwv, mats, wts) + ZZx, x = PolynomialRing(ZZ, length(monomials[1])) + monomials = [finish(push_term!(MPolyBuildCtx(ZZx), ZZ(1), convert(Vector{Int}, mon))) for mon in monomials] + monomials = Set(monomials) + return monomials +end + +nullMon(m) = zeros(Short, m) + + +function compute(v0, mats, wts::Vector{Vector{Int}}) + m = length(mats) + monomials = [nullMon(m)] + lastPos = 0 + id(mon) = sum((1 << (sum(mon[1:i])+i-1) for i in 1:m-1) , init = 1) + e = [Short.(1:m .== i) for i in 1:m] + maxid(deg) = id(deg.*e[1]) + + blacklists = [falses(maxid(0))] + blacklists[end][id(monomials[1])] = 1 + newMons(deg) = (push!(blacklists, falses(maxid(deg)))) + checkMon(mon) = blacklists[end-1][id(mon)] + setMon(mon) = (blacklists[end][id(mon)] = true) + + vectors = [v0] + weights = [0 * wts[1]] + space = Dict(weights[1] => nullSpace()) + addAndReduce!(space[weights[1]], v0) + + deg = 0 + while length(monomials) > lastPos + + startPos = lastPos + 1 + newPos = length(monomials) + deg = deg + 1 + newMons(deg) + for i in 1:m, di in deg:-1:1 + for p in startPos:newPos + + if !all(monomials[p][1:i-1] .== 0) + continue + end + + if monomials[p][i]+1 > di + startPos = p+1 + continue + end + if monomials[p][i]+1 < di + break + end + + mon = monomials[p] + e[i] + + if any(i != j && mon[j] > 0 && !checkMon(mon-e[j]) for j in 1:m) + continue + end + + wt = weights[p] + wts[i] + if !haskey(space, wt) + space[wt] = nullSpace() + end + + vec = mul(vectors[p], transpose(mats[i])) + vec = addAndReduce!(space[wt], vec) + if vec == 0 + continue + end + + setMon(mon) + push!(monomials, mon) + push!(weights, wt) + push!(vectors, vec) + end + end + lastPos = newPos + end + + return monomials +end + +end diff --git a/experimental/BasisLieHighestWeight/test/runtests.jl b/experimental/BasisLieHighestWeight/test/runtests.jl new file mode 100644 index 000000000000..2c1b167a535e --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/runtests.jl @@ -0,0 +1,86 @@ +using Oscar +using Test +# using TestSetExtensions + +include("MBOld.jl") + +G = Oscar.GAP.Globals +forGap = Oscar.GAP.julia_to_gap +fromGap = Oscar.GAP.gap_to_julia + +""" +We are testing our code in multiple ways. First, we calculated two small examples per hand and compare those. Then we +check basic properties of the result. For example we know the size of our monomial basis. These properties get partially +used in the algorithm and could therefore be true for false results. We have another basic algorithm that solves the +problem without the recursion, weightspaces and saving of computations. The third test compares the results we can +compute with the weaker version. +""" + +function compare_algorithms(dynkin::Char, n::Int64, lambda::Vector{Int64}) + # old algorithm + mons_old = MBOld.basisLieHighestWeight(string(dynkin), n, lambda) # basic algorithm + + # new algorithm + mons_new = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda) + L = G.SimpleLieAlgebra(forGap(string(dynkin)), n, G.Rationals) + gap_dim = G.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension + + # comparison + # convert set of monomials over different ring objects to string representation to compare for equality + @test issetequal(string.(mons_old), string.(mons_new)) # compare if result of old and new algorithm match + @test gap_dim == length(mons_new) # check if dimension is correct +end + +function check_dimension(dynkin::Char, n::Int64, lambda::Vector{Int64}, monomial_order::String) + w = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda, monomial_order=monomial_order) + L = G.SimpleLieAlgebra(forGap(string(dynkin)), n, G.Rationals) + gap_dim = G.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension + @test gap_dim == length(w) # check if dimension is correct +end + +@testset "Test basisLieHighestWeight" begin + @testset "Known examples" begin + mons = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0]) + @test issetequal(string.(mons), Set(["1", "x3", "x1"])) + mons = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0], operators=[1,2,1]) + @test issetequal(string.(mons), Set(["1", "x2*x3", "x3"])) + end + @testset "Compare with simple algorithm and check dimension" begin + @testset "Dynkin type $dynkin" for dynkin in ('A', 'B', 'C', 'D') + @testset "n = $n" for n in 1:4 + if (!(dynkin == 'B' && n < 2) && !(dynkin == 'C' && n < 2) && !(dynkin == 'D' && n < 4)) + for i in 1:n # w_i + lambda = zeros(Int64,n) + lambda[i] = 1 + compare_algorithms(dynkin, n, lambda) + end + + if (n > 1) + lambda = [1, (0 for i in 1:n-2)..., 1] # w_1 + w_n + compare_algorithms(dynkin, n, lambda) + end + + if (n < 4) + lambda = ones(Int64,n) # w_1 + ... + w_n + compare_algorithms(dynkin, n, lambda) + end + end + end + end + end + @testset "Check dimension" begin + @testset "Monomial order $monomial_order" for monomial_order in ("Lex", "RevLex", "GRevLex") + # the functionality longest-word was temporarily removed because it required coxeter groups from + # https://github.com/jmichel7/Gapjm.jl + #@testset "Operators $ops" for ops in ("regular", "longest-word") + check_dimension('A', 3, [1,1,1], monomial_order) + #check_dimension('B', 3, [2,1,0], monomial_order, ops) + #check_dimension('C', 3, [1,1,1], monomial_order, ops) + #check_dimension('D', 4, [3,0,1,1], monomial_order, ops) + #check_dimension('F', 4, [2,0,1,0], monomial_order, ops) + #check_dimension('G', 2, [1,0], monomial_order, ops) + #check_dimension('G', 2, [2,2], monomial_order, ops) + #end + end + end +end From f362ac99a8462ad7f8d9d0d02cf88ed26e924bfa Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Fri, 14 Jul 2023 23:18:33 +0200 Subject: [PATCH 02/15] Removed part of fromGap, gens(ZZx) instaed of x --- .../src/BasisLieHighestWeight.jl | 57 +++++++++++-------- .../BasisLieHighestWeight/src/LieAlgebras.jl | 36 ++++++++++-- .../BasisLieHighestWeight/src/NewMonomial.jl | 9 --- .../src/RootConversion.jl | 13 +---- .../src/VectorSpaceBases.jl | 3 - .../BasisLieHighestWeight/src/WeylPolytope.jl | 24 +------- .../BasisLieHighestWeight/test/runtests.jl | 12 ++-- 7 files changed, 71 insertions(+), 83 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 40669c734a1f..06cbfef62003 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -2,9 +2,21 @@ module BasisLieHighestWeight export basis_lie_highest_weight export is_fundamental +using ..Oscar +using ..Oscar: GAPWrap using Polymake +# TODO Change operators from GAP.Obj to Oscar objects +# TODO Doctests +# TODO + +include("./VectorSpaceBases.jl") include("./NewMonomial.jl") +include("./TensorModels.jl") +include("./LieAlgebras.jl") +include("./MonomialOrder.jl") +include("./RootConversion.jl") +include("./WeylPolytope.jl") fromGap = Oscar.GAP.gap_to_julia @@ -101,7 +113,7 @@ function basis_lie_highest_weight(type::String, rank::Int, highest_weight::Vecto monomial_order::Union{String, Function} = "GRevLex", cache_size::Int = 0, parallel::Bool = false, return_no_minkowski::Bool = false, return_operators::Bool = false) - """ + """ Pseudocode: basis_lie_highest_weight(highest_weight) @@ -138,14 +150,14 @@ function basis_lie_highest_weight(type::String, rank::Int, highest_weight::Vecto # steps. Then it starts the recursion and returns the result. # initialization of objects that can be precomputed - lie_algebra, chevalley_basis = create_lie_lgebra(type, rank) # lie_algebra of type, rank and its chevalley_basis + lie_algebra, chevalley_basis = create_lie_algebra(type, rank) # lie_algebra of type, rank and its chevalley_basis # operators that are represented by our monomials. x_i is connected to operators[i] operators = get_operators(type, rank, operators, lie_algebra, chevalley_basis) weights = weights_for_operators(lie_algebra, chevalley_basis[3], operators) # weights of the operators weights = (weight->Int.(weight)).(weights) weights_eps = [w_to_eps(type, rank, w) for w in weights] # other root system - ZZx, x = PolynomialRing(ZZ, length(operators)) # for our monomials - monomial_order_lt = get_monomial_order_lt(monomial_order, ZZx, x) # less than function to sort monomials by order + ZZx, _ = PolynomialRing(ZZ, length(operators)) # for our monomials + monomial_order_lt = get_monomial_order_lt(monomial_order, ZZx, gens(ZZx)) # less than function to sort monomials by order # save computations from recursions calc_highest_weight = Dict{Vector{Int}, Set{ZZMPolyRingElem}}([0 for i in 1:rank] => Set([ZZx(1)])) @@ -153,7 +165,7 @@ function basis_lie_highest_weight(type::String, rank::Int, highest_weight::Vecto no_minkowski = Set{Vector{Int}}() # start recursion over highest_weight - monomial_basis = compute_monomials(type, rank, lie_algebra, ZZx, x, highest_weight, operators, weights, + monomial_basis = compute_monomials(type, rank, lie_algebra, ZZx, highest_weight, operators, weights, weights_eps, monomial_order_lt, calc_highest_weight, cache_size, parallel, no_minkowski) # output @@ -217,7 +229,7 @@ function get_operators(type::String, rank::Int, operators::Union{String, Vector{ return operators end -function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, x::Vector{ZZMPolyRingElem}, +function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, highest_weight::Vector{Int}, operators::GAP.Obj, weights::Vector{Vector{Int64}}, weights_eps::Vector{Vector{Int64}}, monomial_order_lt::Function, calc_highest_weight::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, cache_size::Int, @@ -248,7 +260,7 @@ function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::Z gap_dim = GAP.Globals.DimensionOfHighestWeightModule(lie_algebra, GAP.Obj(highest_weight)) # fundamental weights if is_fundamental(highest_weight) || sum(abs.(highest_weight)) == 0 push!(no_minkowski, highest_weight) - set_mon = add_by_hand(type, rank, lie_algebra, ZZx, x, highest_weight, operators, weights, weights_eps, + set_mon = add_by_hand(type, rank, lie_algebra, ZZx, highest_weight, operators, weights, weights_eps, monomial_order_lt, gap_dim, Set{ZZMPolyRingElem}(), cache_size, parallel) push!(calc_highest_weight, highest_weight => set_mon) return set_mon @@ -264,10 +276,10 @@ function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::Z i += 1 lambda_1 = sub_weights[i] lambda_2 = highest_weight .- lambda_1 - mon_lambda_1 = compute_monomials(type, rank, lie_algebra, ZZx, x, lambda_1, operators, weights, weights_eps, + mon_lambda_1 = compute_monomials(type, rank, lie_algebra, ZZx, lambda_1, operators, weights, weights_eps, monomial_order_lt, calc_highest_weight, cache_size, parallel, no_minkowski) - mon_lambda_2 = compute_monomials(type, rank, lie_algebra, ZZx, x, lambda_2, operators, weights, weights_eps, + mon_lambda_2 = compute_monomials(type, rank, lie_algebra, ZZx, lambda_2, operators, weights, weights_eps, monomial_order_lt, calc_highest_weight, cache_size, parallel, no_minkowski) # Minkowski-sum: M_{lambda_1} + M_{lambda_2} \subseteq M_{highest_weight}, if monomials get identified with @@ -278,7 +290,7 @@ function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::Z # check if we found enough monomials if length(set_mon) < gap_dim push!(no_minkowski, highest_weight) - set_mon = add_by_hand(type, rank, lie_algebra, ZZx, x, highest_weight, operators, weights, weights_eps, + set_mon = add_by_hand(type, rank, lie_algebra, ZZx, highest_weight, operators, weights, weights_eps, monomial_order_lt, gap_dim, set_mon, cache_size, parallel) end push!(calc_highest_weight, highest_weight => set_mon) @@ -311,7 +323,7 @@ function is_fundamental(highest_weight::Vector{Int})::Bool end end end - return false + return one end function compute_sub_weights(highest_weight::Vector{Int})::Vector{Vector{Int}} @@ -331,7 +343,7 @@ end function add_known_monomials!(weight::Vector{Int}, set_mon_in_weightspace::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, number_of_operators::Int, weights::Vector{Vector{Int64}}, matrices_of_operators::Vector{SMat{ZZRingElem}}, calc_monomials::Dict{ZZMPolyRingElem, - Tuple{TVec, Vector{Int}}}, x::Vector{ZZMPolyRingElem}, + Tuple{TVec, Vector{Int}}}, space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}}, v0::SRow{ZZRingElem}, cache_size::Int) """ @@ -345,7 +357,7 @@ function add_known_monomials!(weight::Vector{Int}, set_mon_in_weightspace::Dict{ d = sz(matrices_of_operators[1]) vec = calc_vec(v0, mon, matrices_of_operators) else - vec = calc_new_mon!(x , mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, + vec = calc_new_mon!(gens(ZZx) , mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, cache_size) end @@ -357,7 +369,7 @@ function add_known_monomials!(weight::Vector{Int}, set_mon_in_weightspace::Dict{ end end -function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, x::Vector{ZZMPolyRingElem}, +function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, matrices_of_operators::Vector{SMat{ZZRingElem}}, number_of_operators::Int, weights::Vector{Vector{Int}}, monomial_order_lt::Function, weight::Vector{Int}, dim_weightspace::Int, weights_eps::Vector{Vector{Int}}, @@ -373,7 +385,7 @@ function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, x::Vector polytope is bounded these are finitely many and we can sort them and then go trough them, until we found enough. """ - # get monomials from weyl-polytope that are in the weightspace, sorted by monomial_order_lt + # get monomials that are in the weightspace, sorted by monomial_order_lt poss_mon_in_weightspace = convert_lattice_points_to_monomials(ZZx, get_lattice_points_of_weightspace(weights_eps, w_to_eps(type, rank, weight), type)) poss_mon_in_weightspace = sort(poss_mon_in_weightspace, lt=monomial_order_lt) @@ -398,10 +410,9 @@ function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, x::Vector d = sz(matrices_of_operators[1]) vec = calc_vec(v0, mon, matrices_of_operators) else - vec = calc_new_mon!(x , mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, + vec = calc_new_mon!(gens(ZZx), mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, cache_size) end - #println("vec:" , vec) # check if vec extends the basis if !haskey(space, weight) @@ -419,7 +430,7 @@ function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, x::Vector end -function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, x::Vector{ZZMPolyRingElem}, +function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, highest_weight::Vector{Int}, operators::GAP.Obj, weights::Vector{Vector{Int64}}, weights_eps::Vector{Vector{Int64}}, monomial_order_lt::Function, gap_dim::Int, set_mon::Set{ZZMPolyRingElem}, cache_size::Int, parallel::Bool)::Set{ZZMPolyRingElem} @@ -427,7 +438,6 @@ function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPoly This function calculates the missing monomials by going through each non full weightspace and adding possible monomials manually by computing their corresponding vectors and checking if they enlargen the basis. """ - #println("add_by_hand: ", highest_weight) # initialization # matrices g_i for (g_1^a_1 * ... * g_k^a_k)*v matrices_of_operators = tensorMatricesForOperators(lie_algebra, highest_weight, operators) @@ -466,12 +476,12 @@ function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPoly # insert known monomials into basis for (weight, _) in weightspaces add_known_monomials!(weight, set_mon_in_weightspace, number_of_operators, weights, matrices_of_operators, - calc_monomials, x, space, e, v0, cache_size) + calc_monomials, space, e, v0, cache_size) end # calculate new monomials for (weight, dim_weightspace) in weightspaces - add_new_monomials!(type, rank, ZZx, x, matrices_of_operators, number_of_operators, weights, monomial_order_lt, + add_new_monomials!(type, rank, ZZx, matrices_of_operators, number_of_operators, weights, monomial_order_lt, weight, dim_weightspace, weights_eps, set_mon_in_weightspace, calc_monomials, space, e, v0, cache_size, set_mon) end @@ -487,8 +497,9 @@ function get_dim_weightspace(type::String, rank::Int, lie_algebra::GAP.Obj, """ # calculate dimension for dominant weights with GAP root_system = GAP.Globals.RootSystem(lie_algebra) - dominant_weights, dominant_weights_dim = fromGap(GAP.Globals.DominantCharacter(root_system, - GAP.Obj(highest_weight))) + result = GAP.Globals.DominantCharacter(root_system, GAP.Obj(highest_weight)) + dominant_weights = [map(Int, item) for item in result[1]] + dominant_weights_dim = map(Int, result[2]) dominant_weights = convert(Vector{Vector{Int}}, dominant_weights) weightspaces = Dict{Vector{Int}, Int}() diff --git a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl index d5255631a525..a25386409c0c 100644 --- a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl +++ b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl @@ -1,9 +1,7 @@ -using Oscar - fromGap = Oscar.GAP.gap_to_julia -function create_lie_lgebra(type::String, rank::Int)::Tuple{GAP.Obj, GAP.Obj} +function create_lie_algebra(type::String, rank::Int)::Tuple{GAP.Obj, GAP.Obj} """ Creates the Lie-algebra as a GAP object that gets used for a lot of other computations with GAP """ @@ -23,12 +21,12 @@ function multiply_scalar(A::SMat{T}, d) where T end function matricesForOperators(lie_algebra::GAP.Obj, highest_weight::Vector{Int}, - ops::GAP.Obj)::Vector{SMat{ZZRingElem}} + operators::GAP.Obj)::Vector{SMat{ZZRingElem}} """ used to create tensorMatricesForOperators """ M = Oscar.GAP.Globals.HighestWeightModule(lie_algebra, Oscar.GAP.julia_to_gap(highest_weight)) - matrices_of_operators = Oscar.GAP.Globals.List(ops, o -> Oscar.GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o)) + matrices_of_operators = Oscar.GAP.Globals.List(operators, o -> Oscar.GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o)) matrices_of_operators = gapReshape.( Oscar.GAP.gap_to_julia(matrices_of_operators)) denominators = map(y->denominator(y[2]), union(union(matrices_of_operators...)...)) common_denominator = lcm(denominators)# // 1 @@ -39,15 +37,41 @@ end function weights_for_operators(lie_algebra::GAP.Obj, cartan::GAP.Obj, operators::GAP.Obj)::Vector{Vector{Int}} """ - Calculates the weight wts[i] for each operator ops[i] + Calculates the weight weights[i] for each operator operators[i] """ + """cartan = [Vector{Int}(x) for x in GAP.Globals.ExtRepOfObj.(cartan)] + operators = [Vector{Int}(x) for x in GAP.Globals.ExtRepOfObj.(operators)]# + if any(iszero.(operators)) + error("ops should be non-zero") + end + println([findfirst(v .!= 0) for v in operators]) + + return [ + [(dot(h, v))[findfirst(v .!= 0)] / (v)[findfirst(v .!= 0)] for h in cartan] for v in operators + ] + + + """ + # TODO delete fromGap. Multiplication of cartan and operators is not regular matrix multiplication cartan = fromGap(cartan, recursive=false) operators = fromGap(operators, recursive=false) asVec(v) = fromGap(GAP.Globals.ExtRepOfObj(v)) + #println(cartan) + #println(operators) if any(iszero.(asVec.(operators))) error("ops should be non-zero") end nzi(v) = findfirst(asVec(v) .!= 0) + #println([nzi(v) for v in operators]) + for h in cartan + for v in operators + #println("-") + #println(asVec(v)) + #println(asVec(h)) + #println(asVec(h*v)) + end + end + return [ [asVec(h*v)[nzi(v)] / asVec(v)[nzi(v)] for h in cartan] for v in operators ] diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl index cba90354ce9f..6d0e9ce21885 100644 --- a/experimental/BasisLieHighestWeight/src/NewMonomial.jl +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -1,12 +1,3 @@ -include("./VectorSpaceBases.jl") -include("./TensorModels.jl") -include("./LieAlgebras.jl") -include("./MonomialOrder.jl") -include("./WeylPolytope.jl") - -fromGap = Oscar.GAP.gap_to_julia - - function calc_weight(mon::ZZMPolyRingElem, weights::Vector{Vector{Int}})::Vector{Int} """ calculates weight associated with monomial mon diff --git a/experimental/BasisLieHighestWeight/src/RootConversion.jl b/experimental/BasisLieHighestWeight/src/RootConversion.jl index 6e1191b2d89f..34d2be40a1f8 100644 --- a/experimental/BasisLieHighestWeight/src/RootConversion.jl +++ b/experimental/BasisLieHighestWeight/src/RootConversion.jl @@ -1,8 +1,3 @@ -#using Gapjm -using Oscar - -fromGap = Oscar.GAP.gap_to_julia - function w_to_eps(type::String, rank::Int, weight::Vector{Int})::Vector{Int} """ converts weight in rootsystem w_i to eps_i @@ -76,13 +71,7 @@ end function get_CartanMatrix(type::String, rank::Int) L = GAP.Globals.SimpleLieAlgebra(GAP.Obj(type), rank, GAP.Globals.Rationals) R = GAP.Globals.RootSystem(L) - C_list = fromGap(GAP.Globals.CartanMatrix(R)) - C = zeros(rank, rank) - for i in 1:rank - for j in 1:rank - C[i,j] = C_list[i][j] - end - end + C = Matrix{Int}(GAP.Globals.CartanMatrix(R)) return C end diff --git a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl index f4cf0694a3a5..a5f3bdb417a8 100644 --- a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl +++ b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl @@ -1,8 +1,5 @@ # manages the linear (in-)dependence of integer vectors # this file is only of use to basis_lie_highest_weight for the basis vectors - -using Oscar - TVec = SRow{ZZRingElem} # TVec is datatype of basisvectors in basisLieHighestWeight Short = UInt8 # for exponents of monomials; max. 255 diff --git a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl index d830702b0228..49f7c485b5fd 100644 --- a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl +++ b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl @@ -1,25 +1,3 @@ -#using Gapjm -using Oscar - -include("./RootConversion.jl") - -fromGap = Oscar.GAP.gap_to_julia - -function weylpolytope(type::String, rank::Int, lie_algebra::GAP.Obj, - highest_weight::Vector{Int})::Polymake.BigObjectAllocated - """ - returns weyl-polytope in homogeneous coordinates, i.e. convex hull of orbit of weyl-group of - type type and rank rank on highest weight vector highest_weight - """ - vertices = orbit_weylgroup(type, rank, lie_algebra, highest_weight) - vertices = transpose(hcat(vertices ...)) - vertices = [w for w in eachrow(vertices)] - vertices = transpose(hcat(vertices ...)) - vertices_hom = [ones(Int64, size(vertices)[1]) vertices] # in homogeneous coordinates - weylpoly = polytope.Polytope(POINTS=vertices_hom) - return weylpoly -end - function orbit_weylgroup(type::String, rank::Int, lie_algebra::GAP.Obj, weight_vector::Vector{Int}) """ operates weyl-group of type type and rank rank on vector weight_vector and returns list of vectors in orbit @@ -34,7 +12,7 @@ function orbit_weylgroup(type::String, rank::Int, lie_algebra::GAP.Obj, weight_v GAP.Globals.IsDoneIterator(orbit_iterator) while !(GAP.Globals.IsDoneIterator(orbit_iterator)) w = GAP.Globals.NextIterator(orbit_iterator) - push!(vertices, fromGap(w)) + push!(vertices, Vector{Int}(w)) end # return diff --git a/experimental/BasisLieHighestWeight/test/runtests.jl b/experimental/BasisLieHighestWeight/test/runtests.jl index 2c1b167a535e..2e1a353586d2 100644 --- a/experimental/BasisLieHighestWeight/test/runtests.jl +++ b/experimental/BasisLieHighestWeight/test/runtests.jl @@ -3,10 +3,7 @@ using Test # using TestSetExtensions include("MBOld.jl") - -G = Oscar.GAP.Globals forGap = Oscar.GAP.julia_to_gap -fromGap = Oscar.GAP.gap_to_julia """ We are testing our code in multiple ways. First, we calculated two small examples per hand and compare those. Then we @@ -22,19 +19,20 @@ function compare_algorithms(dynkin::Char, n::Int64, lambda::Vector{Int64}) # new algorithm mons_new = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda) - L = G.SimpleLieAlgebra(forGap(string(dynkin)), n, G.Rationals) - gap_dim = G.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension + L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) + gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension # comparison # convert set of monomials over different ring objects to string representation to compare for equality @test issetequal(string.(mons_old), string.(mons_new)) # compare if result of old and new algorithm match @test gap_dim == length(mons_new) # check if dimension is correct + print(".") end function check_dimension(dynkin::Char, n::Int64, lambda::Vector{Int64}, monomial_order::String) w = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda, monomial_order=monomial_order) - L = G.SimpleLieAlgebra(forGap(string(dynkin)), n, G.Rationals) - gap_dim = G.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension + L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) + gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension @test gap_dim == length(w) # check if dimension is correct end From 073b3ec5a107b5649727f20455519f1da08061c5 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 12:45:00 +0200 Subject: [PATCH 03/15] Simplified paramters with structures and recomputation --- .../src/BasisLieHighestWeight.jl | 239 +++++++++++------- .../src/MonomialOrder.jl | 4 +- .../BasisLieHighestWeight/src/NewMonomial.jl | 33 ++- .../BasisLieHighestWeight/src/WeylPolytope.jl | 12 +- .../BasisLieHighestWeight/test/runtests.jl | 14 +- 5 files changed, 187 insertions(+), 115 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 06cbfef62003..2952bd316561 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -8,7 +8,36 @@ using Polymake # TODO Change operators from GAP.Obj to Oscar objects # TODO Doctests -# TODO +# TODO BirationalSequence needs better names for weights and operators +# TODO Adapt arguments in function outside of BasisLieHighestWeight.jl + +struct LieAlgebra + lie_type::String + rank::Int + lie_algebra_gap::GAP.Obj +end + +struct BirationalSequence + operators::GAP.Obj # TODO Integer + operators_vectors::Vector{Vector{Any}} + weights::Vector{Vector{Int}} + weights_eps::Vector{Vector{Int}} +end + +struct MonomialBasis + set_mon::Set{ZZMPolyRingElem} + dimension::Int + no_minkowski::Set{Vector{Int}} + # polytope::polytope +end + +struct BasisLieHighestWeightStructure + lie_algebra::LieAlgebra + birational_sequence::BirationalSequence + highest_weight::Vector{Int} + monomial_order::Union{String, Function} + monomial_basis::MonomialBasis +end include("./VectorSpaceBases.jl") include("./NewMonomial.jl") @@ -23,8 +52,7 @@ fromGap = Oscar.GAP.gap_to_julia @doc """ basisLieHighestWeight(type::String, rank::Int, highest_weight::Vector{Int}; operators::Union{String, Vector{Int}} = "regular", - monomial_order::Union{String, Function} = "GRevLex", cache_size::Int = 0, - parallel::Bool = false, return_no_minkowski::Bool = false, + monomial_order::Union{String, Function} = "GRevLex", cache_size::Int = 0, return_no_minkowski::Bool = false, return_operators::Bool = false) Computes a monomial basis for the highest weight module with highest weight @@ -40,10 +68,6 @@ Computes a monomial basis for the highest weight module with highest weight groups need a method to obtain all non left descending elements to extend a word - `monomial_order`: monomial order in which our basis gets defined with regards to our operators - `cache_size`: number of computed monomials we want to cache, default is 0 -- `parallel`: currently not implemented, because we used Distributed.jl, but if true parts of the algorithms can be - parallelized -- `return_no_minkowski`: if true return monomials for which Monkowski-property did not suffice to find all monomials -- `return_operators`: if true return the GAP objects operators # Examples ```jldoctest @@ -108,12 +132,16 @@ Set{ZZMPolyRingElem} with 512 elements: ⋮ ``` """ -function basis_lie_highest_weight(type::String, rank::Int, highest_weight::Vector{Int}; - operators::Union{String, Vector{Int}} = "regular", - monomial_order::Union{String, Function} = "GRevLex", cache_size::Int = 0, - parallel::Bool = false, return_no_minkowski::Bool = false, - return_operators::Bool = false) - """ +function basis_lie_highest_weight( + type::String, + rank::Int, + highest_weight::Vector{Int}; + operators::Union{String, Vector{Int}} = "regular", + monomial_order::Union{String, Function} = "GRevLex", + cache_size::Int = 0, + )::BasisLieHighestWeightStructure + + """ Pseudocode: basis_lie_highest_weight(highest_weight) @@ -150,14 +178,19 @@ function basis_lie_highest_weight(type::String, rank::Int, highest_weight::Vecto # steps. Then it starts the recursion and returns the result. # initialization of objects that can be precomputed - lie_algebra, chevalley_basis = create_lie_algebra(type, rank) # lie_algebra of type, rank and its chevalley_basis + lie_algebra_gap, chevalley_basis = create_lie_algebra(type, rank) # lie_algebra of type, rank and its chevalley_basis + lie_algebra = LieAlgebra(type, rank, lie_algebra_gap) # operators that are represented by our monomials. x_i is connected to operators[i] - operators = get_operators(type, rank, operators, lie_algebra, chevalley_basis) - weights = weights_for_operators(lie_algebra, chevalley_basis[3], operators) # weights of the operators + operators = get_operators(lie_algebra, operators, chevalley_basis) + weights = weights_for_operators(lie_algebra.lie_algebra_gap, chevalley_basis[3], operators) # weights of the operators weights = (weight->Int.(weight)).(weights) weights_eps = [w_to_eps(type, rank, w) for w in weights] # other root system + + asVec(v) = fromGap(GAP.Globals.ExtRepOfObj(v)) # TODO + birational_sequence = BirationalSequence(operators, [asVec(v) for v in operators], weights, weights_eps) + ZZx, _ = PolynomialRing(ZZ, length(operators)) # for our monomials - monomial_order_lt = get_monomial_order_lt(monomial_order, ZZx, gens(ZZx)) # less than function to sort monomials by order + monomial_order_lt = get_monomial_order_lt(monomial_order, ZZx) # less than function to sort monomials by order # save computations from recursions calc_highest_weight = Dict{Vector{Int}, Set{ZZMPolyRingElem}}([0 for i in 1:rank] => Set([ZZx(1)])) @@ -165,32 +198,33 @@ function basis_lie_highest_weight(type::String, rank::Int, highest_weight::Vecto no_minkowski = Set{Vector{Int}}() # start recursion over highest_weight - monomial_basis = compute_monomials(type, rank, lie_algebra, ZZx, highest_weight, operators, weights, - weights_eps, monomial_order_lt, calc_highest_weight, cache_size, parallel, no_minkowski) + set_mon = compute_monomials(lie_algebra, birational_sequence, ZZx, highest_weight, monomial_order_lt, calc_highest_weight, cache_size, no_minkowski) # output - if return_no_minkowski && return_operators - return monomial_basis, no_minkowski, operators - elseif return_no_minkowski - return monomial_basis, no_minkowski - elseif return_operators - return monomial_basis, operators - else - return monomial_basis - end + return BasisLieHighestWeightStructure( + lie_algebra, + birational_sequence, + highest_weight, + monomial_order, + MonomialBasis( + set_mon, + length(set_mon), + no_minkowski + ) + ) end -function sub_simple_refl(word::Vector{Int}, lie_algebra::GAP.Obj)::GAP.Obj +function sub_simple_refl(word::Vector{Int}, lie_algebra_gap::GAP.Obj)::GAP.Obj """ substitute simple reflections (i,i+1), saved in dec by i, with E_{i,i+1} """ - root_system = GAP.Globals.RootSystem(lie_algebra) + root_system = GAP.Globals.RootSystem(lie_algebra_gap) canonical_generators = fromGap(GAP.Globals.CanonicalGenerators(root_system)[1], recursive = false) operators = GAP.Obj([canonical_generators[i] for i in word], recursive = false) return operators end -function get_operators(type::String, rank::Int, operators::Union{String, Vector{Int}}, lie_algebra::GAP.Obj, +function get_operators(lie_algebra::LieAlgebra, operators::Union{String, Vector{Int}}, chevalley_basis::GAP.Obj)::GAP.Obj """ handles user input for operators @@ -218,22 +252,26 @@ function get_operators(type::String, rank::Int, operators::Union{String, Vector{ println("operators needs to be of type Vector{Int}") return -1 end - if !(all([(1 <= i <= rank) for i in operators])) + if !(all([(1 <= i <= lie_algebra.rank) for i in operators])) println("all values of operators need to between 1 and the rank of the lie algebra.") end # If one of the conditions is met, the algorithms works. Otherwise a warning is printed (and can be ignored). #if !(is_longest_weyl_word(type, rank, operators)) && !(Set(operators) == [i for i=1:n]) # println("WARNING: operators may be incorrect input.") #end - operators = sub_simple_refl(operators, lie_algebra) + operators = sub_simple_refl(operators, lie_algebra.lie_algebra_gap) return operators end -function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, - highest_weight::Vector{Int}, operators::GAP.Obj, weights::Vector{Vector{Int64}}, - weights_eps::Vector{Vector{Int64}}, monomial_order_lt::Function, - calc_highest_weight::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, cache_size::Int, - parallel::Bool, no_minkowski::Set{Vector{Int}})::Set{ZZMPolyRingElem} +function compute_monomials( + lie_algebra::LieAlgebra, + birational_sequence::BirationalSequence, + ZZx::ZZMPolyRing, + highest_weight::Vector{Int}, + monomial_order_lt::Function, + calc_highest_weight::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, + cache_size::Int, + no_minkowski::Set{Vector{Int}})::Set{ZZMPolyRingElem} """ This function calculates the monomial basis M_{highest_weight} recursively. The recursion saves all computed results in calc_highest_weight and we first check, if we already encountered this highest weight in a prior step. @@ -249,7 +287,7 @@ function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::Z # we already computed the highest_weight result in a prior recursion step if haskey(calc_highest_weight, highest_weight) return calc_highest_weight[highest_weight] - elseif highest_weight == [0 for i in 1:rank] # we mathematically know the solution + elseif highest_weight == [0 for i in 1:lie_algebra.rank] # we mathematically know the solution return Set(ZZx(1)) end @@ -257,11 +295,11 @@ function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::Z # gap_dim is number of monomials that we need to find, i.e. |M_{highest_weight}|. # if highest_weight is a fundamental weight, partition into smaller summands is possible. This is the basecase of # the recursion. - gap_dim = GAP.Globals.DimensionOfHighestWeightModule(lie_algebra, GAP.Obj(highest_weight)) # fundamental weights + gap_dim = GAP.Globals.DimensionOfHighestWeightModule(lie_algebra.lie_algebra_gap, GAP.Obj(highest_weight)) # fundamental weights if is_fundamental(highest_weight) || sum(abs.(highest_weight)) == 0 push!(no_minkowski, highest_weight) - set_mon = add_by_hand(type, rank, lie_algebra, ZZx, highest_weight, operators, weights, weights_eps, - monomial_order_lt, gap_dim, Set{ZZMPolyRingElem}(), cache_size, parallel) + set_mon = add_by_hand(lie_algebra, birational_sequence, ZZx, highest_weight, + monomial_order_lt, gap_dim, Set{ZZMPolyRingElem}(), cache_size) push!(calc_highest_weight, highest_weight => set_mon) return set_mon else @@ -272,15 +310,15 @@ function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::Z l = length(sub_weights) # go through all partitions lambda_1 + lambda_2 = highest_weight until we have enough monomials or used all # partitions - while length(set_mon) < gap_dim && i < l + while length(set_mon) < gap_dim && i < l i += 1 lambda_1 = sub_weights[i] lambda_2 = highest_weight .- lambda_1 - mon_lambda_1 = compute_monomials(type, rank, lie_algebra, ZZx, lambda_1, operators, weights, weights_eps, - monomial_order_lt, calc_highest_weight, cache_size, parallel, + mon_lambda_1 = compute_monomials(lie_algebra, birational_sequence, ZZx, lambda_1, + monomial_order_lt, calc_highest_weight, cache_size, no_minkowski) - mon_lambda_2 = compute_monomials(type, rank, lie_algebra, ZZx, lambda_2, operators, weights, weights_eps, - monomial_order_lt, calc_highest_weight, cache_size, parallel, + mon_lambda_2 = compute_monomials(lie_algebra, birational_sequence, ZZx, lambda_2, + monomial_order_lt, calc_highest_weight, cache_size, no_minkowski) # Minkowski-sum: M_{lambda_1} + M_{lambda_2} \subseteq M_{highest_weight}, if monomials get identified with # points in ZZ^n @@ -290,8 +328,8 @@ function compute_monomials(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::Z # check if we found enough monomials if length(set_mon) < gap_dim push!(no_minkowski, highest_weight) - set_mon = add_by_hand(type, rank, lie_algebra, ZZx, highest_weight, operators, weights, weights_eps, - monomial_order_lt, gap_dim, set_mon, cache_size, parallel) + set_mon = add_by_hand(lie_algebra, birational_sequence, ZZx, highest_weight, + monomial_order_lt, gap_dim, set_mon, cache_size) end push!(calc_highest_weight, highest_weight => set_mon) return set_mon @@ -328,7 +366,7 @@ end function compute_sub_weights(highest_weight::Vector{Int})::Vector{Vector{Int}} """ - returns list of weights w != 0 with 0 <= w <= highest_weight elementwise + returns list of weights w != 0 with 0 <= w <= highest_weight elementwise, ordered by l_2-norm """ sub_weights = [] foreach(Iterators.product((0:x for x in highest_weight)...)) do i @@ -340,12 +378,17 @@ function compute_sub_weights(highest_weight::Vector{Int})::Vector{Vector{Int}} return sub_weights end -function add_known_monomials!(weight::Vector{Int}, set_mon_in_weightspace::Dict{Vector{Int64}, - Set{ZZMPolyRingElem}}, number_of_operators::Int, weights::Vector{Vector{Int64}}, - matrices_of_operators::Vector{SMat{ZZRingElem}}, calc_monomials::Dict{ZZMPolyRingElem, - Tuple{TVec, Vector{Int}}}, - space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}}, - v0::SRow{ZZRingElem}, cache_size::Int) +function add_known_monomials!( + birational_sequence::BirationalSequence, + ZZx::ZZMPolyRing, + weight::Vector{Int}, + set_mon_in_weightspace::Dict{Vector{Int64}, + Set{ZZMPolyRingElem}}, + matrices_of_operators::Vector{SMat{ZZRingElem}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, + v0::SRow{ZZRingElem}, + cache_size::Int) """ By using the Minkowski-sum, we know that all monomials in set_mon_in_weightspace are in our basis. Since we want to extend the weightspace with missing monomials, we need to calculate and add the vector of each monomial to our @@ -357,7 +400,7 @@ function add_known_monomials!(weight::Vector{Int}, set_mon_in_weightspace::Dict{ d = sz(matrices_of_operators[1]) vec = calc_vec(v0, mon, matrices_of_operators) else - vec = calc_new_mon!(gens(ZZx) , mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, + vec = calc_new_mon!(gens(ZZx) , mon, birational_sequence.weights, matrices_of_operators, calc_monomials, space, cache_size) end @@ -369,14 +412,21 @@ function add_known_monomials!(weight::Vector{Int}, set_mon_in_weightspace::Dict{ end end -function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, - matrices_of_operators::Vector{SMat{ZZRingElem}}, number_of_operators::Int, - weights::Vector{Vector{Int}}, monomial_order_lt::Function, weight::Vector{Int}, - dim_weightspace::Int, weights_eps::Vector{Vector{Int}}, - set_mon_in_weightspace::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, - calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, space::Dict{Vector{Int64}, - Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}}, v0::SRow{ZZRingElem}, - cache_size::Int, set_mon::Set{ZZMPolyRingElem}) +function add_new_monomials!( + lie_algebra::LieAlgebra, + birational_sequence::BirationalSequence, + ZZx::ZZMPolyRing, + matrices_of_operators::Vector{SMat{ZZRingElem}}, + monomial_order_lt::Function, + dim_weightspace::Int, + weight::Vector{Int}, + set_mon_in_weightspace::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + space::Dict{Vector{Int64}, + Oscar.BasisLieHighestWeight.VSBasis}, + v0::SRow{ZZRingElem}, + cache_size::Int, + set_mon::Set{ZZMPolyRingElem}) """ If a weightspace is missing monomials, we need to calculate them by trial and error. We would like to go through all monomials in the order monomial_order_lt and calculate the corresponding vector. If it extends the basis, we add it @@ -386,8 +436,8 @@ function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, """ # get monomials that are in the weightspace, sorted by monomial_order_lt - poss_mon_in_weightspace = convert_lattice_points_to_monomials(ZZx, - get_lattice_points_of_weightspace(weights_eps, w_to_eps(type, rank, weight), type)) + poss_mon_in_weightspace = convert_lattice_points_to_monomials(ZZx, get_lattice_points_of_weightspace( + birational_sequence.weights_eps, w_to_eps(lie_algebra.lie_type, lie_algebra.rank, weight), lie_algebra.lie_type)) poss_mon_in_weightspace = sort(poss_mon_in_weightspace, lt=monomial_order_lt) # check which monomials should get added to the basis @@ -410,7 +460,7 @@ function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, d = sz(matrices_of_operators[1]) vec = calc_vec(v0, mon, matrices_of_operators) else - vec = calc_new_mon!(gens(ZZx), mon, weights, matrices_of_operators, number_of_operators, calc_monomials, space, e, + vec = calc_new_mon!(gens(ZZx), mon, birational_sequence.weights, matrices_of_operators, calc_monomials, space, cache_size) end @@ -430,26 +480,30 @@ function add_new_monomials!(type::String, rank::Int, ZZx::ZZMPolyRing, end -function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPolyRing, - highest_weight::Vector{Int}, operators::GAP.Obj, weights::Vector{Vector{Int64}}, - weights_eps::Vector{Vector{Int64}}, monomial_order_lt::Function, gap_dim::Int, - set_mon::Set{ZZMPolyRingElem}, cache_size::Int, parallel::Bool)::Set{ZZMPolyRingElem} +function add_by_hand( + lie_algebra::LieAlgebra, + birational_sequence::BirationalSequence, + ZZx::ZZMPolyRing, + highest_weight::Vector{Int}, + monomial_order_lt::Function, + gap_dim::Int, + set_mon::Set{ZZMPolyRingElem}, + cache_size::Int, + )::Set{ZZMPolyRingElem} """ This function calculates the missing monomials by going through each non full weightspace and adding possible monomials manually by computing their corresponding vectors and checking if they enlargen the basis. """ # initialization # matrices g_i for (g_1^a_1 * ... * g_k^a_k)*v - matrices_of_operators = tensorMatricesForOperators(lie_algebra, highest_weight, operators) - number_of_operators = length(matrices_of_operators) - e = [1*(1:number_of_operators .== i) for i in 1:number_of_operators] # e_i - space = Dict(0*weights[1] => nullSpace()) # span of basis vectors to keep track of the basis + matrices_of_operators = tensorMatricesForOperators(lie_algebra.lie_algebra_gap, highest_weight, birational_sequence.operators) + space = Dict(0*birational_sequence.weights[1] => nullSpace()) # span of basis vectors to keep track of the basis v0 = sparse_row(ZZ, [(1,1)]) # starting vector v # saves the calculated vectors to decrease necessary matrix multiplicatons - calc_monomials = Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}(ZZx(1) => (v0, 0 * weights[1])) + calc_monomials = Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}(ZZx(1) => (v0, 0 * birational_sequence.weights[1])) push!(set_mon, ZZx(1)) # required monomials of each weightspace - weightspaces = get_dim_weightspace(type, rank, lie_algebra, highest_weight) + weightspaces = get_dim_weightspace(lie_algebra, highest_weight) # sort the monomials from the minkowski-sum by their weightspaces set_mon_in_weightspace = Dict{Vector{Int}, Set{ZZMPolyRingElem}}() @@ -457,7 +511,7 @@ function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPoly set_mon_in_weightspace[weight] = Set{ZZMPolyRingElem}() end for mon in set_mon - weight = calc_weight(mon, weights) + weight = calc_weight(mon, birational_sequence.weights) push!(set_mon_in_weightspace[weight], mon) end @@ -469,34 +523,39 @@ function add_by_hand(type::String, rank::Int, lie_algebra::GAP.Obj, ZZx::ZZMPoly end end delete!(weightspaces, weights_with_full_weightspace) - - # use parallel computations if parallel = true. The weightspaces could be calculated completely indepent (except for + + # The weightspaces could be calculated completely indepent (except for # the caching). This is not implemented, since I used the package Distributed.jl for this, which is not in the # Oscar dependencies. But I plan to reimplement this. # insert known monomials into basis for (weight, _) in weightspaces - add_known_monomials!(weight, set_mon_in_weightspace, number_of_operators, weights, matrices_of_operators, - calc_monomials, space, e, v0, cache_size) + add_known_monomials!(birational_sequence, ZZx, weight, set_mon_in_weightspace, matrices_of_operators, + calc_monomials, space, v0, cache_size) end # calculate new monomials for (weight, dim_weightspace) in weightspaces - add_new_monomials!(type, rank, ZZx, matrices_of_operators, number_of_operators, weights, monomial_order_lt, - weight, dim_weightspace, weights_eps, set_mon_in_weightspace, calc_monomials, space, e, v0, - cache_size, set_mon) + add_new_monomials!(lie_algebra, birational_sequence, ZZx, matrices_of_operators, + monomial_order_lt, + dim_weightspace, weight, + set_mon_in_weightspace, + calc_monomials, space, v0, + cache_size, set_mon) end return set_mon end -function get_dim_weightspace(type::String, rank::Int, lie_algebra::GAP.Obj, - highest_weight::Vector{Int})::Dict{Vector{Int}, Int} +function get_dim_weightspace( + lie_algebra::LieAlgebra, + highest_weight::Vector{Int} + )::Dict{Vector{Int}, Int} """ Calculates dictionary with weights as keys and dimension of corresponding weightspace as value. GAP computes the dimension for all positive weights. The dimension is constant on orbits of the weylgroup, and we can therefore calculate the dimension of each weightspace. """ # calculate dimension for dominant weights with GAP - root_system = GAP.Globals.RootSystem(lie_algebra) + root_system = GAP.Globals.RootSystem(lie_algebra.lie_algebra_gap) result = GAP.Globals.DominantCharacter(root_system, GAP.Obj(highest_weight)) dominant_weights = [map(Int, item) for item in result[1]] dominant_weights_dim = map(Int, result[2]) @@ -505,7 +564,7 @@ function get_dim_weightspace(type::String, rank::Int, lie_algebra::GAP.Obj, # calculate dimension for the rest by checking which positive weights lies in the orbit. for i in 1:length(dominant_weights) - orbit_weights = orbit_weylgroup(type, rank, lie_algebra, dominant_weights[i]) + orbit_weights = orbit_weylgroup(lie_algebra, dominant_weights[i]) dim_weightspace = dominant_weights_dim[i] for weight in orbit_weights weightspaces[highest_weight - weight] = dim_weightspace diff --git a/experimental/BasisLieHighestWeight/src/MonomialOrder.jl b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl index 5185fd8f3ebf..de08134fabd2 100644 --- a/experimental/BasisLieHighestWeight/src/MonomialOrder.jl +++ b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl @@ -1,10 +1,10 @@ -function get_monomial_order_lt(monomial_order::Union{String, Function}, ZZx::ZZMPolyRing, - x::Vector{ZZMPolyRingElem})::Function +function get_monomial_order_lt(monomial_order::Union{String, Function}, ZZx::ZZMPolyRing)::Function """ Returns the desired monomial_order function less than """ #if isa(monomial_order, Function) # choosen_monomial_order = monomial_order + x = gens(ZZx) if monomial_order == "GRevLex" choosen_monomial_order = degrevlex(x) elseif monomial_order == "RevLex" diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl index 6d0e9ce21885..e4d1128c96d8 100644 --- a/experimental/BasisLieHighestWeight/src/NewMonomial.jl +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -10,8 +10,11 @@ function calc_weight(mon::ZZMPolyRingElem, weights::Vector{Vector{Int}})::Vector return weight end -function calc_vec(v0::SRow{ZZRingElem}, mon::ZZMPolyRingElem, - matrices_of_operators::Vector{SMat{ZZRingElem}})::SRow{ZZRingElem} +function calc_vec( + v0::SRow{ZZRingElem}, + mon::ZZMPolyRingElem, + matrices_of_operators::Vector{SMat{ZZRingElem}} + )::SRow{ZZRingElem} """ calculates vector associated with monomial mon """ @@ -25,14 +28,16 @@ function calc_vec(v0::SRow{ZZRingElem}, mon::ZZMPolyRingElem, return vec end -function highest_calc_sub_monomial(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, - calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, - number_of_operators::Int)::ZZMPolyRingElem +function highest_calc_sub_monomial( + x::Vector{ZZMPolyRingElem}, + mon::ZZMPolyRingElem, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + )::ZZMPolyRingElem """ returns the key in calc_monomials that can be extended by the least amount of left-operations to mon """ sub_mon = copy(mon) - + number_of_operators = length(mon) for i in 1:number_of_operators while is_divisible_by(sub_mon, x[i]) if haskey(calc_monomials, sub_mon) @@ -45,16 +50,20 @@ function highest_calc_sub_monomial(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingE return sub_mon end -function calc_new_mon!(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, weights::Vector{Vector{Int}}, - matrices_of_operators::Vector{SMat{ZZRingElem}}, number_of_operators::Int, - calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, - space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}}, - cache_size::Int)::SRow{ZZRingElem} +function calc_new_mon!(x::Vector{ZZMPolyRingElem}, + mon::ZZMPolyRingElem, + weights::Vector{Vector{Int}}, + matrices_of_operators::Vector{SMat{ZZRingElem}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, + cache_size::Int + )::SRow{ZZRingElem} # calculate vector of mon by extending a previous calculated vector to a # monom that differs only by left-multiplication, save results in calc_monomials - sub_mon = highest_calc_sub_monomial(x, mon, calc_monomials, number_of_operators) + sub_mon = highest_calc_sub_monomial(x, mon, calc_monomials) #println("sub_mon: ", sub_mon) sub_mon_cur = copy(sub_mon) + number_of_operators = length(mon) (vec, weight) = calc_monomials[sub_mon] for i in number_of_operators:-1:1 for k in degrees(sub_mon)[i]:(degrees(mon)[i]-1) diff --git a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl index 49f7c485b5fd..5a48ffa3565d 100644 --- a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl +++ b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl @@ -1,10 +1,10 @@ -function orbit_weylgroup(type::String, rank::Int, lie_algebra::GAP.Obj, weight_vector::Vector{Int}) +function orbit_weylgroup(lie_algebra::LieAlgebra, weight_vector::Vector{Int}) """ operates weyl-group of type type and rank rank on vector weight_vector and returns list of vectors in orbit input and output weights in terms of w_i """ # initialization - weyl_group = GAP.Globals.WeylGroup(GAP.Globals.RootSystem(lie_algebra)) + weyl_group = GAP.Globals.WeylGroup(GAP.Globals.RootSystem(lie_algebra.lie_algebra_gap)) orbit_iterator = GAP.Globals.WeylOrbitIterator(weyl_group, GAP.Obj(weight_vector)) vertices = [] @@ -36,7 +36,7 @@ function convert_lattice_points_to_monomials(ZZx, lattice_points_weightspace) for lattice_point in lattice_points_weightspace] end -function get_lattice_points_of_weightspace(weights, weight, type) +function get_lattice_points_of_weightspace(weights, weight, lie_type) """ calculates all lattice points in a given weightspace for a lie algebra of type type input: @@ -45,7 +45,7 @@ function get_lattice_points_of_weightspace(weights, weight, type) output: all lattice points with weight weight """ - if type in ["A", "G"] + if lie_type in ["A", "G"] return get_lattice_points_of_weightspace_A_G_n(weights, weight) else return get_lattice_points_of_weightspace_Xn(weights, weight) @@ -71,7 +71,7 @@ function get_lattice_points_of_weightspace_A_G_n(weights, weight) example: weights = [[1, 0, 2], [-1, 1, 1], [0, -1, 0]] (i.e. a_1 = eps_1 - eps_2, a_2 = eps_2 - eps_3, a_12 = eps_1 - eps_3) weight = [2, 1, 0] - -> poly = polytope.Polytope(INEQUALITIES=[0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1], + -> poly = polytope.polytope(INEQUALITIES=[0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1], EQUATIONS=[-2 1 1 0 2; -1 1 -1 1 1; 0 1 0 -1 0]) => returns [[1 0 0], [1 1 0]] """ @@ -97,7 +97,7 @@ end function get_lattice_points_of_weightspace_Xn(weights, weight) """ - calculates all lattice points in a given weightspace for lie algebras that don't have type A + calculates all lattice points in a given weightspace for lie algebras that don't have type A or G input: weights: the operator weights in eps_i weight: lambda - mu diff --git a/experimental/BasisLieHighestWeight/test/runtests.jl b/experimental/BasisLieHighestWeight/test/runtests.jl index 2e1a353586d2..798878d552d6 100644 --- a/experimental/BasisLieHighestWeight/test/runtests.jl +++ b/experimental/BasisLieHighestWeight/test/runtests.jl @@ -18,7 +18,8 @@ function compare_algorithms(dynkin::Char, n::Int64, lambda::Vector{Int64}) mons_old = MBOld.basisLieHighestWeight(string(dynkin), n, lambda) # basic algorithm # new algorithm - mons_new = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda) + base = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda) + mons_new = base.monomial_basis.set_mon L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension @@ -30,17 +31,20 @@ function compare_algorithms(dynkin::Char, n::Int64, lambda::Vector{Int64}) end function check_dimension(dynkin::Char, n::Int64, lambda::Vector{Int64}, monomial_order::String) - w = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda, monomial_order=monomial_order) + base = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda, monomial_order=monomial_order) + mons_new = base.monomial_basis.set_mon L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension - @test gap_dim == length(w) # check if dimension is correct + @test gap_dim == length(mons_new) # check if dimension is correct end @testset "Test basisLieHighestWeight" begin @testset "Known examples" begin - mons = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0]) + base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0]) + mons = base.monomial_basis.set_mon @test issetequal(string.(mons), Set(["1", "x3", "x1"])) - mons = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0], operators=[1,2,1]) + base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0], operators=[1,2,1]) + mons = base.monomial_basis.set_mon @test issetequal(string.(mons), Set(["1", "x2*x3", "x3"])) end @testset "Compare with simple algorithm and check dimension" begin From 01ce341c300a50d427b9c8ef9b943a5c991883e2 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 13:15:52 +0200 Subject: [PATCH 04/15] TODOs --- .../src/BasisLieHighestWeight.jl | 45 +++++++++++++++++-- .../BasisLieHighestWeight/src/LieAlgebras.jl | 2 - .../BasisLieHighestWeight/src/WeylPolytope.jl | 11 ----- 3 files changed, 41 insertions(+), 17 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 2952bd316561..779137be550f 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -7,9 +7,42 @@ using ..Oscar: GAPWrap using Polymake # TODO Change operators from GAP.Obj to Oscar objects -# TODO Doctests # TODO BirationalSequence needs better names for weights and operators +# TODO Fix weightnames in general # TODO Adapt arguments in function outside of BasisLieHighestWeight.jl +# TODO Summarize function +# TODO Groundup-structure for the special functions +# TODO Bugfix cache_size != 0 + +# TODO Export and docstring: +# basis_lie_highest_weight +# get_dim_weightspace +# orbit_weylgroup +# get_lattice_points_of_weightspace +# convert_lattice_points_to_monomials +# convert_monomials_to_lattice_points + +# tensorMatricesForOperators +# weights_for_operators + +# w_to_eps +# eps_to_w +# alpha_to_eps +# eps_to_alpha +# w_to_eps +# eps_to_w + +# Unittests for +# VectorSpaceBases.jl: +# reduce_col +# normalize +# add_and_reduce! + +# NewMonomial.jl +# calc_weight +# calc_vec +# highest_calc_sub_monomial +# calc_new_mon! struct LieAlgebra lie_type::String @@ -50,10 +83,14 @@ include("./WeylPolytope.jl") fromGap = Oscar.GAP.gap_to_julia @doc """ - basisLieHighestWeight(type::String, rank::Int, highest_weight::Vector{Int}; +basis_lie_highest_weight( + type::String, + rank::Int, + highest_weight::Vector{Int}; operators::Union{String, Vector{Int}} = "regular", - monomial_order::Union{String, Function} = "GRevLex", cache_size::Int = 0, return_no_minkowski::Bool = false, - return_operators::Bool = false) + monomial_order::Union{String, Function} = "GRevLex", + cache_size::Int = 0, +)::BasisLieHighestWeightStructure Computes a monomial basis for the highest weight module with highest weight ``highest_weight`` (in terms of the fundamental weights), for a simple Lie algebra of type diff --git a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl index a25386409c0c..3818dd6b47e3 100644 --- a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl +++ b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl @@ -49,8 +49,6 @@ function weights_for_operators(lie_algebra::GAP.Obj, cartan::GAP.Obj, operators: return [ [(dot(h, v))[findfirst(v .!= 0)] / (v)[findfirst(v .!= 0)] for h in cartan] for v in operators ] - - """ # TODO delete fromGap. Multiplication of cartan and operators is not regular matrix multiplication cartan = fromGap(cartan, recursive=false) diff --git a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl index 5a48ffa3565d..42427711d1af 100644 --- a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl +++ b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl @@ -20,17 +20,6 @@ function orbit_weylgroup(lie_algebra::LieAlgebra, weight_vector::Vector{Int}) return vertices end -function get_points_polytope(polytope) - """ - returns all points (interior and vertices) of a polytope in regular (i.e. not homogenoues coordinates). - """ - interior_points = convert(Matrix{Int64}, polytope.INTERIOR_LATTICE_POINTS) - vertices_points = convert(Matrix{Int64}, polytope.VERTICES) - points = [interior_points; vertices_points][:, 2:end] - return points -end - - function convert_lattice_points_to_monomials(ZZx, lattice_points_weightspace) return [finish(push_term!(MPolyBuildCtx(ZZx), ZZ(1), convert(Vector{Int}, convert(Vector{Int64}, lattice_point)))) for lattice_point in lattice_points_weightspace] From addd9854090bd734d6ce54d3d0a90d2128170f89 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 14:13:42 +0200 Subject: [PATCH 05/15] Added test-files, removed TVec --- .../docs/src/basisLieHighestWeight.md | 2 +- .../src/BasisLieHighestWeight.jl | 10 +- .../BasisLieHighestWeight/src/NewMonomial.jl | 4 +- .../src/VectorSpaceBases.jl | 10 +- .../test/BasisLieHighestWeight-test.jl | 88 ++++++++++++++++++ .../test/NewMonomial-test.jl | 10 ++ .../test/VectorSpaceBases-test.jl | 10 ++ .../BasisLieHighestWeight/test/runtests.jl | 91 +------------------ 8 files changed, 125 insertions(+), 100 deletions(-) create mode 100644 experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl create mode 100644 experimental/BasisLieHighestWeight/test/NewMonomial-test.jl create mode 100644 experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl diff --git a/experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md b/experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md index 7fc64f8ce43c..434c63f185c6 100644 --- a/experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md +++ b/experimental/BasisLieHighestWeight/docs/src/basisLieHighestWeight.md @@ -12,5 +12,5 @@ Pages = ["basisLieHighestWeight.md"] # Monomial bases for Lie algebras ```@docs -basisLieHighestWeight2 +BasisLieHighestWeight ``` diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 779137be550f..e7e83fbb60cb 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -13,6 +13,7 @@ using Polymake # TODO Summarize function # TODO Groundup-structure for the special functions # TODO Bugfix cache_size != 0 +# TODO Remove TVec # TODO Export and docstring: # basis_lie_highest_weight @@ -33,6 +34,9 @@ using Polymake # eps_to_w # Unittests for +# BasisLieHighestWeight.jl +# compute_sub_weights + # VectorSpaceBases.jl: # reduce_col # normalize @@ -422,7 +426,7 @@ function add_known_monomials!( set_mon_in_weightspace::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, matrices_of_operators::Vector{SMat{ZZRingElem}}, - calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, v0::SRow{ZZRingElem}, cache_size::Int) @@ -458,7 +462,7 @@ function add_new_monomials!( dim_weightspace::Int, weight::Vector{Int}, set_mon_in_weightspace::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, - calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, v0::SRow{ZZRingElem}, @@ -537,7 +541,7 @@ function add_by_hand( space = Dict(0*birational_sequence.weights[1] => nullSpace()) # span of basis vectors to keep track of the basis v0 = sparse_row(ZZ, [(1,1)]) # starting vector v # saves the calculated vectors to decrease necessary matrix multiplicatons - calc_monomials = Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}(ZZx(1) => (v0, 0 * birational_sequence.weights[1])) + calc_monomials = Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}(ZZx(1) => (v0, 0 * birational_sequence.weights[1])) push!(set_mon, ZZx(1)) # required monomials of each weightspace weightspaces = get_dim_weightspace(lie_algebra, highest_weight) diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl index e4d1128c96d8..0e717597c167 100644 --- a/experimental/BasisLieHighestWeight/src/NewMonomial.jl +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -31,7 +31,7 @@ end function highest_calc_sub_monomial( x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, - calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, )::ZZMPolyRingElem """ returns the key in calc_monomials that can be extended by the least amount of left-operations to mon @@ -54,7 +54,7 @@ function calc_new_mon!(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, weights::Vector{Vector{Int}}, matrices_of_operators::Vector{SMat{ZZRingElem}}, - calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}}, + calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, cache_size::Int )::SRow{ZZRingElem} diff --git a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl index a5f3bdb417a8..24d3a9c82385 100644 --- a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl +++ b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl @@ -1,18 +1,16 @@ # manages the linear (in-)dependence of integer vectors # this file is only of use to basis_lie_highest_weight for the basis vectors -TVec = SRow{ZZRingElem} # TVec is datatype of basisvectors in basisLieHighestWeight -Short = UInt8 # for exponents of monomials; max. 255 struct VSBasis - basis_vectors::Vector{TVec} # vector of basisvectors + basis_vectors::Vector{SRow{ZZRingElem}} # vector of basisvectors pivot::Vector{Int} # vector of pivotelements, i.e. pivot[i] is first nonzero element of basis_vectors[i] end nullSpace() = VSBasis([], []) # empty Vektorraum -reduce_col(a::TVec, b::TVec, i::Int) = (b[i]*a - a[i]*b)::TVec # create zero entry in a +reduce_col(a::SRow{ZZRingElem}, b::SRow{ZZRingElem}, i::Int) = (b[i]*a - a[i]*b)::SRow{ZZRingElem} # create zero entry in a -function normalize(v::TVec)::Tuple{TVec, Int64} +function normalize(v::SRow{ZZRingElem})::Tuple{SRow{ZZRingElem}, Int64} """ divides vector by gcd of nonzero entries, returns vector and first nonzero index used: add_and_reduce! @@ -24,7 +22,7 @@ function normalize(v::TVec)::Tuple{TVec, Int64} return divexact(v, gcd(map(y->y[2], union(v)))), pivot end -function add_and_reduce!(sp::VSBasis, v::TVec)::TVec +function add_and_reduce!(sp::VSBasis, v::SRow{ZZRingElem})::SRow{ZZRingElem} """ for each pivot of sp.basis_vectors we make entry of v zero and return the result and insert it into sp 0 => linear dependent diff --git a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl new file mode 100644 index 000000000000..798878d552d6 --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl @@ -0,0 +1,88 @@ +using Oscar +using Test +# using TestSetExtensions + +include("MBOld.jl") +forGap = Oscar.GAP.julia_to_gap + +""" +We are testing our code in multiple ways. First, we calculated two small examples per hand and compare those. Then we +check basic properties of the result. For example we know the size of our monomial basis. These properties get partially +used in the algorithm and could therefore be true for false results. We have another basic algorithm that solves the +problem without the recursion, weightspaces and saving of computations. The third test compares the results we can +compute with the weaker version. +""" + +function compare_algorithms(dynkin::Char, n::Int64, lambda::Vector{Int64}) + # old algorithm + mons_old = MBOld.basisLieHighestWeight(string(dynkin), n, lambda) # basic algorithm + + # new algorithm + base = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda) + mons_new = base.monomial_basis.set_mon + L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) + gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension + + # comparison + # convert set of monomials over different ring objects to string representation to compare for equality + @test issetequal(string.(mons_old), string.(mons_new)) # compare if result of old and new algorithm match + @test gap_dim == length(mons_new) # check if dimension is correct + print(".") +end + +function check_dimension(dynkin::Char, n::Int64, lambda::Vector{Int64}, monomial_order::String) + base = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda, monomial_order=monomial_order) + mons_new = base.monomial_basis.set_mon + L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) + gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension + @test gap_dim == length(mons_new) # check if dimension is correct +end + +@testset "Test basisLieHighestWeight" begin + @testset "Known examples" begin + base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0]) + mons = base.monomial_basis.set_mon + @test issetequal(string.(mons), Set(["1", "x3", "x1"])) + base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0], operators=[1,2,1]) + mons = base.monomial_basis.set_mon + @test issetequal(string.(mons), Set(["1", "x2*x3", "x3"])) + end + @testset "Compare with simple algorithm and check dimension" begin + @testset "Dynkin type $dynkin" for dynkin in ('A', 'B', 'C', 'D') + @testset "n = $n" for n in 1:4 + if (!(dynkin == 'B' && n < 2) && !(dynkin == 'C' && n < 2) && !(dynkin == 'D' && n < 4)) + for i in 1:n # w_i + lambda = zeros(Int64,n) + lambda[i] = 1 + compare_algorithms(dynkin, n, lambda) + end + + if (n > 1) + lambda = [1, (0 for i in 1:n-2)..., 1] # w_1 + w_n + compare_algorithms(dynkin, n, lambda) + end + + if (n < 4) + lambda = ones(Int64,n) # w_1 + ... + w_n + compare_algorithms(dynkin, n, lambda) + end + end + end + end + end + @testset "Check dimension" begin + @testset "Monomial order $monomial_order" for monomial_order in ("Lex", "RevLex", "GRevLex") + # the functionality longest-word was temporarily removed because it required coxeter groups from + # https://github.com/jmichel7/Gapjm.jl + #@testset "Operators $ops" for ops in ("regular", "longest-word") + check_dimension('A', 3, [1,1,1], monomial_order) + #check_dimension('B', 3, [2,1,0], monomial_order, ops) + #check_dimension('C', 3, [1,1,1], monomial_order, ops) + #check_dimension('D', 4, [3,0,1,1], monomial_order, ops) + #check_dimension('F', 4, [2,0,1,0], monomial_order, ops) + #check_dimension('G', 2, [1,0], monomial_order, ops) + #check_dimension('G', 2, [2,2], monomial_order, ops) + #end + end + end +end diff --git a/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl new file mode 100644 index 000000000000..0da9693decc7 --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl @@ -0,0 +1,10 @@ +using Oscar +using Test + +include("../src/NewMonomial.jl") + +@testset "Test NewMonomial" begin + @testset "2" begin + @test isequal(1, 1) + end +end diff --git a/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl b/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl new file mode 100644 index 000000000000..b7d4953d7e6b --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl @@ -0,0 +1,10 @@ +using Oscar +using Test + +include("../src/VectorSpaceBases.jl") + +@testset "Test VectorSpaceBases" begin + @testset "1" begin + @test isequal(1, 1) + end +end diff --git a/experimental/BasisLieHighestWeight/test/runtests.jl b/experimental/BasisLieHighestWeight/test/runtests.jl index 798878d552d6..60548b8cbb35 100644 --- a/experimental/BasisLieHighestWeight/test/runtests.jl +++ b/experimental/BasisLieHighestWeight/test/runtests.jl @@ -1,88 +1,3 @@ -using Oscar -using Test -# using TestSetExtensions - -include("MBOld.jl") -forGap = Oscar.GAP.julia_to_gap - -""" -We are testing our code in multiple ways. First, we calculated two small examples per hand and compare those. Then we -check basic properties of the result. For example we know the size of our monomial basis. These properties get partially -used in the algorithm and could therefore be true for false results. We have another basic algorithm that solves the -problem without the recursion, weightspaces and saving of computations. The third test compares the results we can -compute with the weaker version. -""" - -function compare_algorithms(dynkin::Char, n::Int64, lambda::Vector{Int64}) - # old algorithm - mons_old = MBOld.basisLieHighestWeight(string(dynkin), n, lambda) # basic algorithm - - # new algorithm - base = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda) - mons_new = base.monomial_basis.set_mon - L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) - gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension - - # comparison - # convert set of monomials over different ring objects to string representation to compare for equality - @test issetequal(string.(mons_old), string.(mons_new)) # compare if result of old and new algorithm match - @test gap_dim == length(mons_new) # check if dimension is correct - print(".") -end - -function check_dimension(dynkin::Char, n::Int64, lambda::Vector{Int64}, monomial_order::String) - base = BasisLieHighestWeight.basis_lie_highest_weight(string(dynkin), n, lambda, monomial_order=monomial_order) - mons_new = base.monomial_basis.set_mon - L = Oscar.GAP.Globals.SimpleLieAlgebra(forGap(string(dynkin)), n, Oscar.GAP.Globals.Rationals) - gap_dim = Oscar.GAP.Globals.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension - @test gap_dim == length(mons_new) # check if dimension is correct -end - -@testset "Test basisLieHighestWeight" begin - @testset "Known examples" begin - base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0]) - mons = base.monomial_basis.set_mon - @test issetequal(string.(mons), Set(["1", "x3", "x1"])) - base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0], operators=[1,2,1]) - mons = base.monomial_basis.set_mon - @test issetequal(string.(mons), Set(["1", "x2*x3", "x3"])) - end - @testset "Compare with simple algorithm and check dimension" begin - @testset "Dynkin type $dynkin" for dynkin in ('A', 'B', 'C', 'D') - @testset "n = $n" for n in 1:4 - if (!(dynkin == 'B' && n < 2) && !(dynkin == 'C' && n < 2) && !(dynkin == 'D' && n < 4)) - for i in 1:n # w_i - lambda = zeros(Int64,n) - lambda[i] = 1 - compare_algorithms(dynkin, n, lambda) - end - - if (n > 1) - lambda = [1, (0 for i in 1:n-2)..., 1] # w_1 + w_n - compare_algorithms(dynkin, n, lambda) - end - - if (n < 4) - lambda = ones(Int64,n) # w_1 + ... + w_n - compare_algorithms(dynkin, n, lambda) - end - end - end - end - end - @testset "Check dimension" begin - @testset "Monomial order $monomial_order" for monomial_order in ("Lex", "RevLex", "GRevLex") - # the functionality longest-word was temporarily removed because it required coxeter groups from - # https://github.com/jmichel7/Gapjm.jl - #@testset "Operators $ops" for ops in ("regular", "longest-word") - check_dimension('A', 3, [1,1,1], monomial_order) - #check_dimension('B', 3, [2,1,0], monomial_order, ops) - #check_dimension('C', 3, [1,1,1], monomial_order, ops) - #check_dimension('D', 4, [3,0,1,1], monomial_order, ops) - #check_dimension('F', 4, [2,0,1,0], monomial_order, ops) - #check_dimension('G', 2, [1,0], monomial_order, ops) - #check_dimension('G', 2, [2,2], monomial_order, ops) - #end - end - end -end +include("VectorSpaceBases-test.jl") +include("NewMonomial-test.jl") +include("BasisLieHighestWeight-test.jl") \ No newline at end of file From ee8fd05d095edc7c471759d9f032403873668594 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 14:18:07 +0200 Subject: [PATCH 06/15] Removed nullSpace() --- .../BasisLieHighestWeight/src/BasisLieHighestWeight.jl | 7 +++---- experimental/BasisLieHighestWeight/src/NewMonomial.jl | 2 +- .../BasisLieHighestWeight/src/VectorSpaceBases.jl | 2 -- experimental/BasisLieHighestWeight/test/MBOld.jl | 8 ++------ 4 files changed, 6 insertions(+), 13 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index e7e83fbb60cb..c1bfc8f4cfec 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -13,7 +13,6 @@ using Polymake # TODO Summarize function # TODO Groundup-structure for the special functions # TODO Bugfix cache_size != 0 -# TODO Remove TVec # TODO Export and docstring: # basis_lie_highest_weight @@ -447,7 +446,7 @@ function add_known_monomials!( # check if vec extends the basis if !haskey(space, weight) - space[weight] = nullSpace() + space[weight] = VSBasis([], []) end add_and_reduce!(space[weight], vec) end @@ -507,7 +506,7 @@ function add_new_monomials!( # check if vec extends the basis if !haskey(space, weight) - space[weight] = nullSpace() + space[weight] = VSBasis([], []) end vec_red = add_and_reduce!(space[weight], vec) if isempty(vec_red) # v0 == 0 @@ -538,7 +537,7 @@ function add_by_hand( # initialization # matrices g_i for (g_1^a_1 * ... * g_k^a_k)*v matrices_of_operators = tensorMatricesForOperators(lie_algebra.lie_algebra_gap, highest_weight, birational_sequence.operators) - space = Dict(0*birational_sequence.weights[1] => nullSpace()) # span of basis vectors to keep track of the basis + space = Dict(0*birational_sequence.weights[1] => VSBasis([], [])) # span of basis vectors to keep track of the basis v0 = sparse_row(ZZ, [(1,1)]) # starting vector v # saves the calculated vectors to decrease necessary matrix multiplicatons calc_monomials = Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}(ZZx(1) => (v0, 0 * birational_sequence.weights[1])) diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl index 0e717597c167..c6ec34632ed0 100644 --- a/experimental/BasisLieHighestWeight/src/NewMonomial.jl +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -70,7 +70,7 @@ function calc_new_mon!(x::Vector{ZZMPolyRingElem}, sub_mon_cur *= x[i] weight += weights[i] if !haskey(space, weight) - space[weight] = nullSpace() + space[weight] = VSBasis([], []) end vec = mul(vec, transpose(matrices_of_operators[i])) # currently there is no sparse matrix * vector mult diff --git a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl index 24d3a9c82385..c5ee66373796 100644 --- a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl +++ b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl @@ -6,8 +6,6 @@ struct VSBasis pivot::Vector{Int} # vector of pivotelements, i.e. pivot[i] is first nonzero element of basis_vectors[i] end -nullSpace() = VSBasis([], []) # empty Vektorraum - reduce_col(a::SRow{ZZRingElem}, b::SRow{ZZRingElem}, i::Int) = (b[i]*a - a[i]*b)::SRow{ZZRingElem} # create zero entry in a function normalize(v::SRow{ZZRingElem})::Tuple{SRow{ZZRingElem}, Int64} diff --git a/experimental/BasisLieHighestWeight/test/MBOld.jl b/experimental/BasisLieHighestWeight/test/MBOld.jl index b990e3d39127..c8fe936d22ee 100644 --- a/experimental/BasisLieHighestWeight/test/MBOld.jl +++ b/experimental/BasisLieHighestWeight/test/MBOld.jl @@ -24,10 +24,6 @@ struct VSBasis pivot::Vector{Int} end - -nullSpace() = VSBasis([], []) - - function normalize(v::TVec) """ divides vector by gcd of nonzero entries, returns vector and first nonzero index @@ -244,7 +240,7 @@ function compute(v0, mats, wts::Vector{Vector{Int}}) vectors = [v0] weights = [0 * wts[1]] - space = Dict(weights[1] => nullSpace()) + space = Dict(weights[1] => VSBasis([], [])) addAndReduce!(space[weights[1]], v0) deg = 0 @@ -277,7 +273,7 @@ function compute(v0, mats, wts::Vector{Vector{Int}}) wt = weights[p] + wts[i] if !haskey(space, wt) - space[wt] = nullSpace() + space[wt] = VSBasis([], []) end vec = mul(vectors[p], transpose(mats[i])) From 1f112ecc8180b210b30a53863bfa6ac8995a6c67 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 17:25:22 +0200 Subject: [PATCH 07/15] VectorSpaceBases-test.jl --- .../src/BasisLieHighestWeight.jl | 15 ++++------ .../BasisLieHighestWeight/src/NewMonomial.jl | 4 +-- .../src/VectorSpaceBases.jl | 12 ++++---- .../test/BasisLieHighestWeight-test.jl | 2 +- .../BasisLieHighestWeight/test/MBOld.jl | 17 +++++------ .../test/VectorSpaceBases-test.jl | 30 +++++++++++++++++-- 6 files changed, 50 insertions(+), 30 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index c1bfc8f4cfec..6b572b86d183 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -36,11 +36,6 @@ using Polymake # BasisLieHighestWeight.jl # compute_sub_weights -# VectorSpaceBases.jl: -# reduce_col -# normalize -# add_and_reduce! - # NewMonomial.jl # calc_weight # calc_vec @@ -426,7 +421,7 @@ function add_known_monomials!( Set{ZZMPolyRingElem}}, matrices_of_operators::Vector{SMat{ZZRingElem}}, calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, - space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, + space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.SparseVectorSpaceBasis}, v0::SRow{ZZRingElem}, cache_size::Int) """ @@ -446,7 +441,7 @@ function add_known_monomials!( # check if vec extends the basis if !haskey(space, weight) - space[weight] = VSBasis([], []) + space[weight] = SparseVectorSpaceBasis([], []) end add_and_reduce!(space[weight], vec) end @@ -463,7 +458,7 @@ function add_new_monomials!( set_mon_in_weightspace::Dict{Vector{Int64}, Set{ZZMPolyRingElem}}, calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, space::Dict{Vector{Int64}, - Oscar.BasisLieHighestWeight.VSBasis}, + Oscar.BasisLieHighestWeight.SparseVectorSpaceBasis}, v0::SRow{ZZRingElem}, cache_size::Int, set_mon::Set{ZZMPolyRingElem}) @@ -506,7 +501,7 @@ function add_new_monomials!( # check if vec extends the basis if !haskey(space, weight) - space[weight] = VSBasis([], []) + space[weight] = SparseVectorSpaceBasis([], []) end vec_red = add_and_reduce!(space[weight], vec) if isempty(vec_red) # v0 == 0 @@ -537,7 +532,7 @@ function add_by_hand( # initialization # matrices g_i for (g_1^a_1 * ... * g_k^a_k)*v matrices_of_operators = tensorMatricesForOperators(lie_algebra.lie_algebra_gap, highest_weight, birational_sequence.operators) - space = Dict(0*birational_sequence.weights[1] => VSBasis([], [])) # span of basis vectors to keep track of the basis + space = Dict(0*birational_sequence.weights[1] => SparseVectorSpaceBasis([], [])) # span of basis vectors to keep track of the basis v0 = sparse_row(ZZ, [(1,1)]) # starting vector v # saves the calculated vectors to decrease necessary matrix multiplicatons calc_monomials = Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}(ZZx(1) => (v0, 0 * birational_sequence.weights[1])) diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl index c6ec34632ed0..fd9b02eecda6 100644 --- a/experimental/BasisLieHighestWeight/src/NewMonomial.jl +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -55,7 +55,7 @@ function calc_new_mon!(x::Vector{ZZMPolyRingElem}, weights::Vector{Vector{Int}}, matrices_of_operators::Vector{SMat{ZZRingElem}}, calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, - space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, + space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.SparseVectorSpaceBasis}, cache_size::Int )::SRow{ZZRingElem} # calculate vector of mon by extending a previous calculated vector to a @@ -70,7 +70,7 @@ function calc_new_mon!(x::Vector{ZZMPolyRingElem}, sub_mon_cur *= x[i] weight += weights[i] if !haskey(space, weight) - space[weight] = VSBasis([], []) + space[weight] = SparseVectorSpaceBasis([], []) end vec = mul(vec, transpose(matrices_of_operators[i])) # currently there is no sparse matrix * vector mult diff --git a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl index c5ee66373796..663fa6cd9095 100644 --- a/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl +++ b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl @@ -1,12 +1,13 @@ # manages the linear (in-)dependence of integer vectors # this file is only of use to basis_lie_highest_weight for the basis vectors -struct VSBasis +struct SparseVectorSpaceBasis basis_vectors::Vector{SRow{ZZRingElem}} # vector of basisvectors pivot::Vector{Int} # vector of pivotelements, i.e. pivot[i] is first nonzero element of basis_vectors[i] end -reduce_col(a::SRow{ZZRingElem}, b::SRow{ZZRingElem}, i::Int) = (b[i]*a - a[i]*b)::SRow{ZZRingElem} # create zero entry in a + # create zero entry in i-th entry +reduce_col(a::SRow{ZZRingElem}, b::SRow{ZZRingElem}, i::Int) = (b[i]*a - a[i]*b)::SRow{ZZRingElem} function normalize(v::SRow{ZZRingElem})::Tuple{SRow{ZZRingElem}, Int64} """ @@ -20,7 +21,7 @@ function normalize(v::SRow{ZZRingElem})::Tuple{SRow{ZZRingElem}, Int64} return divexact(v, gcd(map(y->y[2], union(v)))), pivot end -function add_and_reduce!(sp::VSBasis, v::SRow{ZZRingElem})::SRow{ZZRingElem} +function add_and_reduce!(sp::SparseVectorSpaceBasis, v::SRow{ZZRingElem})::SRow{ZZRingElem} """ for each pivot of sp.basis_vectors we make entry of v zero and return the result and insert it into sp 0 => linear dependent @@ -40,11 +41,10 @@ function add_and_reduce!(sp::VSBasis, v::SRow{ZZRingElem})::SRow{ZZRingElem} # use pivots of basis basis_vectors to create zeros in v for j in 1:length(basis_vectors) - i = pivot[j] - if i != newPivot + if pivot[j] != newPivot continue end - v = reduce_col(v, basis_vectors[j], i) + v = reduce_col(v, basis_vectors[j], newPivot) v, newPivot = normalize(v) if newPivot == 0 #return 0 diff --git a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl index 798878d552d6..945a9e3d1308 100644 --- a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl +++ b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl @@ -49,7 +49,7 @@ end end @testset "Compare with simple algorithm and check dimension" begin @testset "Dynkin type $dynkin" for dynkin in ('A', 'B', 'C', 'D') - @testset "n = $n" for n in 1:4 + @testset "n = $n" for n in 1:5 if (!(dynkin == 'B' && n < 2) && !(dynkin == 'C' && n < 2) && !(dynkin == 'D' && n < 4)) for i in 1:n # w_i lambda = zeros(Int64,n) diff --git a/experimental/BasisLieHighestWeight/test/MBOld.jl b/experimental/BasisLieHighestWeight/test/MBOld.jl index c8fe936d22ee..9d660e8796db 100644 --- a/experimental/BasisLieHighestWeight/test/MBOld.jl +++ b/experimental/BasisLieHighestWeight/test/MBOld.jl @@ -16,15 +16,14 @@ G = Oscar.GAP.Globals forGap = Oscar.GAP.julia_to_gap fromGap = Oscar.GAP.gap_to_julia -TVec = SRow{ZZRingElem} Short = UInt8 -struct VSBasis - A::Vector{TVec} +struct SparseVectorSpaceBasis + A::Vector{SRow{ZZRingElem}} pivot::Vector{Int} end -function normalize(v::TVec) +function normalize(v::SRow{ZZRingElem}) """ divides vector by gcd of nonzero entries, returns vector and first nonzero index used: addAndReduce! @@ -42,7 +41,7 @@ end reduceCol(a, b, i::Int) = b[i]*a - a[i]*b -function addAndReduce!(sp::VSBasis, v::TVec) +function addAndReduce!(sp::SparseVectorSpaceBasis, v::SRow{ZZRingElem}) """ for each pivot of sp.A we make entry of v zero and return the result 0 => linear dependent @@ -192,7 +191,7 @@ end """ - basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; parallel::Bool = true) :: Tuple{Vector{Vector{Short}},Vector{TVec}} + basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; parallel::Bool = true) :: Tuple{Vector{Vector{Short}},Vector{SRow{ZZRingElem}}} Compute a monomial basis for the highest weight module with highest weight ``hw`` (in terms of the fundamental weights), for a simple Lie algebra of type ``t`` and rank ``n``. @@ -204,7 +203,7 @@ julia> dim, monomials, vectors = PolyBases.MB.basisLieHighestWeight("A", 2, [1,0 ``` """ -function basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; roots = []) #--- :: Tuple{Int64,Vector{Vector{Short}},Vector{TVec}} +function basisLieHighestWeight(t::String, n::Int, hw::Vector{Int}; roots = []) #--- :: Tuple{Int64,Vector{Vector{Short}},Vector{SRow{ZZRingElem}}} L, CH = lieAlgebra(t, n) ops = CH[1] # positive root vectors # .. reorder.. @@ -240,7 +239,7 @@ function compute(v0, mats, wts::Vector{Vector{Int}}) vectors = [v0] weights = [0 * wts[1]] - space = Dict(weights[1] => VSBasis([], [])) + space = Dict(weights[1] => SparseVectorSpaceBasis([], [])) addAndReduce!(space[weights[1]], v0) deg = 0 @@ -273,7 +272,7 @@ function compute(v0, mats, wts::Vector{Vector{Int}}) wt = weights[p] + wts[i] if !haskey(space, wt) - space[wt] = VSBasis([], []) + space[wt] = SparseVectorSpaceBasis([], []) end vec = mul(vectors[p], transpose(mats[i])) diff --git a/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl b/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl index b7d4953d7e6b..f2f9abd8eed8 100644 --- a/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl +++ b/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl @@ -4,7 +4,33 @@ using Test include("../src/VectorSpaceBases.jl") @testset "Test VectorSpaceBases" begin - @testset "1" begin - @test isequal(1, 1) + a = sparse_row(ZZ, [1, 3, 6], [3, 2, 4])::SRow{ZZRingElem} # [3, 0, 2, 0, 0, 4] + b = sparse_row(ZZ, [3], [2])::SRow{ZZRingElem} # [0, 0, 2, 0, 0, 0] + c = sparse_row(ZZ, [1, 6], [4, 3])::SRow{ZZRingElem} # [4, 0, 0, 0, 0, 3] + d = sparse_row(ZZ, [1, 3], [4, 3])::SRow{ZZRingElem} # [6, 0, 4, 0, 0, 0] + sparse_vector_space_basis = SparseVectorSpaceBasis([a, b], [1, 3]) + + @testset "reduce_col" begin + a_reduced_b_3 = sparse_row(ZZ, [1, 6], [6, 8])::SRow{ZZRingElem} # [6, 0, 0, 0, 0, 8] + a_reduced_c_1 = sparse_row(ZZ, [3, 6], [8, 7])::SRow{ZZRingElem} # [0, 0, 8, 0, 0, 7] + @test isequal(reduce_col(a, b, 3), a_reduced_b_3) + @test isequal(reduce_col(a, c, 1), a_reduced_c_1) + end + + @testset "normalize" begin + a_normalized = sparse_row(ZZ, [1, 3, 6], [3, 2, 4])::SRow{ZZRingElem} # [3, 0, 2, 0, 0, 4] + b_normalized = sparse_row(ZZ, [3], [1])::SRow{ZZRingElem} # [0, 0, 1, 0, 0, 0] + @test isequal(normalize(a), (a_normalized, 1)) + @test isequal(normalize(b), (b_normalized, 3)) + end + + @testset "add_and_reduce!" begin + add_and_reduce!(sparse_vector_space_basis, c) + c_reduced = sparse_row(ZZ, [6], [1])::SRow{ZZRingElem} # [0, 0, 0, 0, 0, 1] + @test isequal(sparse_vector_space_basis.basis_vectors, [a, b, c_reduced]) + @test isequal(sparse_vector_space_basis.pivot, [1, 3, 6]) + add_and_reduce!(sparse_vector_space_basis, d) # d = 2*a, therefore basis should not change + @test isequal(sparse_vector_space_basis.basis_vectors, [a, b, c_reduced]) + @test isequal(sparse_vector_space_basis.pivot, [1, 3, 6]) end end From d74516e50c17fe452c90fbb6dc612884f9dc9a80 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 19:08:11 +0200 Subject: [PATCH 08/15] NewMonomial-test.jl --- .../src/BasisLieHighestWeight.jl | 10 ++---- .../BasisLieHighestWeight/src/LieAlgebras.jl | 2 ++ .../BasisLieHighestWeight/src/NewMonomial.jl | 10 +++--- .../test/BasisLieHighestWeight-test.jl | 2 +- .../test/NewMonomial-test.jl | 31 +++++++++++++++++-- 5 files changed, 40 insertions(+), 15 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 6b572b86d183..6576e05fe234 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -36,12 +36,6 @@ using Polymake # BasisLieHighestWeight.jl # compute_sub_weights -# NewMonomial.jl -# calc_weight -# calc_vec -# highest_calc_sub_monomial -# calc_new_mon! - struct LieAlgebra lie_type::String rank::Int @@ -324,8 +318,8 @@ function compute_monomials( return calc_highest_weight[highest_weight] elseif highest_weight == [0 for i in 1:lie_algebra.rank] # we mathematically know the solution return Set(ZZx(1)) - end - + end + # calculation required # gap_dim is number of monomials that we need to find, i.e. |M_{highest_weight}|. # if highest_weight is a fundamental weight, partition into smaller summands is possible. This is the basecase of diff --git a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl index 3818dd6b47e3..a25386409c0c 100644 --- a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl +++ b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl @@ -49,6 +49,8 @@ function weights_for_operators(lie_algebra::GAP.Obj, cartan::GAP.Obj, operators: return [ [(dot(h, v))[findfirst(v .!= 0)] / (v)[findfirst(v .!= 0)] for h in cartan] for v in operators ] + + """ # TODO delete fromGap. Multiplication of cartan and operators is not regular matrix multiplication cartan = fromGap(cartan, recursive=false) diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl index fd9b02eecda6..31d7e9ddb86d 100644 --- a/experimental/BasisLieHighestWeight/src/NewMonomial.jl +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -13,7 +13,7 @@ end function calc_vec( v0::SRow{ZZRingElem}, mon::ZZMPolyRingElem, - matrices_of_operators::Vector{SMat{ZZRingElem}} + matrices_of_operators::Union{Vector{SMat{ZZRingElem, Hecke.ZZRingElem_Array_Mod.ZZRingElem_Array}}, Vector{SMat{ZZRingElem}}} )::SRow{ZZRingElem} """ calculates vector associated with monomial mon @@ -22,7 +22,9 @@ function calc_vec( degree_mon = degrees(mon) for i in length(degree_mon):-1:1 for j in 1:degree_mon[i] - vec = mul(vec, transpose(matrices_of_operators[i])) # currently there is no sparse matrix * vector mult + # currently there is no sparse matrix * vector mult + # this is also the line that takes up almost all the computation time for big examples + vec = mul(vec, transpose(matrices_of_operators[i])) end end return vec @@ -37,7 +39,7 @@ function highest_calc_sub_monomial( returns the key in calc_monomials that can be extended by the least amount of left-operations to mon """ sub_mon = copy(mon) - number_of_operators = length(mon) + number_of_operators = length(x) for i in 1:number_of_operators while is_divisible_by(sub_mon, x[i]) if haskey(calc_monomials, sub_mon) @@ -53,7 +55,7 @@ end function calc_new_mon!(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, weights::Vector{Vector{Int}}, - matrices_of_operators::Vector{SMat{ZZRingElem}}, + matrices_of_operators::Union{Vector{SMat{ZZRingElem, Hecke.ZZRingElem_Array_Mod.ZZRingElem_Array}}, Vector{SMat{ZZRingElem}}}, calc_monomials::Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}, space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.SparseVectorSpaceBasis}, cache_size::Int diff --git a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl index 945a9e3d1308..798878d552d6 100644 --- a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl +++ b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl @@ -49,7 +49,7 @@ end end @testset "Compare with simple algorithm and check dimension" begin @testset "Dynkin type $dynkin" for dynkin in ('A', 'B', 'C', 'D') - @testset "n = $n" for n in 1:5 + @testset "n = $n" for n in 1:4 if (!(dynkin == 'B' && n < 2) && !(dynkin == 'C' && n < 2) && !(dynkin == 'D' && n < 4)) for i in 1:n # w_i lambda = zeros(Int64,n) diff --git a/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl index 0da9693decc7..81d7d680042a 100644 --- a/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl +++ b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl @@ -4,7 +4,34 @@ using Test include("../src/NewMonomial.jl") @testset "Test NewMonomial" begin - @testset "2" begin - @test isequal(1, 1) + ZZx, _ = PolynomialRing(ZZ, 2) + x = gens(ZZx) + mon1 = ZZx(1) + mon2 = x[1]^2 * x[2] + weights = [[1, 1], [2, 1]] + A = sparse_matrix(ZZ, 2, 2) # [0, 1; 2, 1] + setindex!(A, sparse_row(ZZ, [2], [ZZ(1)]), 1) + setindex!(A, sparse_row(ZZ, [1, 2], [ZZ(2), ZZ(1)]), 2) + B = sparse_matrix(ZZ, 2, 2) # [1, 0; 2, 0] + setindex!(B, sparse_row(ZZ, [1], [ZZ(1)]), 1) + setindex!(B, sparse_row(ZZ, [2], [ZZ(2)]), 2) + matrices_of_operators = [A, B] + v0 = sparse_row(ZZ, [1], [1])::SRow{ZZRingElem} # [1, 0] + calc_monomials = Dict{ZZMPolyRingElem, Tuple{SRow{ZZRingElem}, Vector{Int}}}(ZZx(1) => (v0, [0, 0])) + + mon2_vec = sparse_row(ZZ, [1, 2], [2, 2])::SRow{ZZRingElem} + + @testset "calc_weight" begin + @test isequal(calc_weight(mon1, weights), [0, 0]) + @test isequal(calc_weight(mon2, weights), [4, 3]) + end + + @testset "calc_vec" begin + @test isequal(calc_vec(v0, mon1, matrices_of_operators), v0) + @test isequal(calc_vec(v0, mon2, matrices_of_operators), mon2_vec) + end + + @testset "highest_calc_sub_monomial" begin + @test isequal(highest_calc_sub_monomial(x, mon2, calc_monomials), ZZx(1)) end end From a74c531d3e3aee2df4c355b7cc7f2bcb3a34023c Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 19:24:11 +0200 Subject: [PATCH 09/15] compute_sub_weights --- .../src/BasisLieHighestWeight.jl | 18 ++++++++--------- .../test/BasisLieHighestWeight-test.jl | 20 +++++++++++++++---- 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 6576e05fe234..0a749692dbbd 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -32,10 +32,6 @@ using Polymake # w_to_eps # eps_to_w -# Unittests for -# BasisLieHighestWeight.jl -# compute_sub_weights - struct LieAlgebra lie_type::String rank::Int @@ -395,16 +391,20 @@ end function compute_sub_weights(highest_weight::Vector{Int})::Vector{Vector{Int}} """ - returns list of weights w != 0 with 0 <= w <= highest_weight elementwise, ordered by l_2-norm + returns list of weights w != 0, highest_weight with 0 <= w <= highest_weight elementwise, ordered by l_2-norm """ sub_weights = [] foreach(Iterators.product((0:x for x in highest_weight)...)) do i push!(sub_weights, [i...]) end - popfirst!(sub_weights) # [0, ..., 0] - pop!(sub_weights) # highest_weight - sort!(sub_weights, by=x->sum((x).^2)) - return sub_weights + if isempty(sub_weights) || length(sub_weights) == 1 # case [] or [[0, ..., 0]] + return [] + else + popfirst!(sub_weights) # [0, ..., 0] + pop!(sub_weights) # highest_weight + sort!(sub_weights, by=x->sum((x).^2)) + return sub_weights + end end function add_known_monomials!( diff --git a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl index 798878d552d6..5d07421c5b6c 100644 --- a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl +++ b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl @@ -1,7 +1,6 @@ using Oscar using Test # using TestSetExtensions - include("MBOld.jl") forGap = Oscar.GAP.julia_to_gap @@ -38,8 +37,21 @@ function check_dimension(dynkin::Char, n::Int64, lambda::Vector{Int64}, monomial @test gap_dim == length(mons_new) # check if dimension is correct end -@testset "Test basisLieHighestWeight" begin - @testset "Known examples" begin +@testset "Test BasisLieHighestWeight" begin + @testset "is_fundamental" begin + @test BasisLieHighestWeight.is_fundamental([0, 1, 0]) + @test !BasisLieHighestWeight.is_fundamental([0, 1, 1]) + end + + @testset "compute_sub_weights" begin + @test isequal(BasisLieHighestWeight.compute_sub_weights([0, 0, 0]), []) + sub_weights = [[1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 0, 1], [0, 1, 1], [1, 1, 1], [2, 0, 0], + [0, 2, 0], [2, 1, 0], [1, 2, 0], [2, 0, 1], [0, 2, 1], [2, 1, 1], [1, 2, 1], [2, 2, 0], + [0, 3, 0], [2, 2, 1], [1, 3, 0], [0, 3, 1], [1, 3, 1], [2, 3, 0]] + @test isequal(BasisLieHighestWeight.compute_sub_weights([2,3,1]), sub_weights) + end + + @testset "Known examples basis_lie_highest_weight" begin base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1,0]) mons = base.monomial_basis.set_mon @test issetequal(string.(mons), Set(["1", "x3", "x1"])) @@ -47,7 +59,7 @@ end mons = base.monomial_basis.set_mon @test issetequal(string.(mons), Set(["1", "x2*x3", "x3"])) end - @testset "Compare with simple algorithm and check dimension" begin + @testset "Compare basis_lie_highest_weight with algorithm of Johannes and check dimension" begin @testset "Dynkin type $dynkin" for dynkin in ('A', 'B', 'C', 'D') @testset "n = $n" for n in 1:4 if (!(dynkin == 'B' && n < 2) && !(dynkin == 'C' && n < 2) && !(dynkin == 'D' && n < 4)) From beaedafbb0131e3afaa5abd87722ffccae73879f Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 22:18:20 +0200 Subject: [PATCH 10/15] Added GAPWrap instead of GAP.Globals where possible --- Project.toml | 1 + .../src/BasisLieHighestWeight.jl | 79 +++++++++---------- .../BasisLieHighestWeight/src/LieAlgebras.jl | 12 +-- .../BasisLieHighestWeight/src/WeylPolytope.jl | 41 +++++++++- 4 files changed, 81 insertions(+), 52 deletions(-) diff --git a/Project.toml b/Project.toml index b480eba139ef..3f792cd1155c 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.12.1-DEV" AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" AlgebraicSolving = "66b61cbe-0446-4d5d-9090-1ff510639f9d" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GAP = "c863536a-3901-11e9-33e7-d5cd0df7b904" Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 0a749692dbbd..5e02dab632c7 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -1,6 +1,9 @@ module BasisLieHighestWeight +export LieAlgebra export basis_lie_highest_weight export is_fundamental +export get_dim_weightspace +export orbit_weylgroup using ..Oscar using ..Oscar: GAPWrap @@ -12,19 +15,16 @@ using Polymake # TODO Adapt arguments in function outside of BasisLieHighestWeight.jl # TODO Summarize function # TODO Groundup-structure for the special functions -# TODO Bugfix cache_size != 0 -# TODO Export and docstring: +# TODO Export and docstring (?): # basis_lie_highest_weight # get_dim_weightspace # orbit_weylgroup # get_lattice_points_of_weightspace # convert_lattice_points_to_monomials # convert_monomials_to_lattice_points - # tensorMatricesForOperators # weights_for_operators - # w_to_eps # eps_to_w # alpha_to_eps @@ -32,12 +32,33 @@ using Polymake # w_to_eps # eps_to_w -struct LieAlgebra +# TODO GAPWrap-wrappers are missing for +# ChevalleyBasis +# DimensionOfHighestWeightModule +# SimpleLieAlgebra +# Rationals +# HighestWeightModule +# List +# MatrixOfAction +# RootSystem +# CartanMatrix +# WeylGroup +# DominantCharacter +# DimensionOfHighestWeightModule +# CanonicalGenerators + + +struct LieAlgebraStructure lie_type::String rank::Int lie_algebra_gap::GAP.Obj end +function LieAlgebraStructure(lie_type::String, rank::Int) + return LieAlgebraStructure(lie_type, rank, create_lie_algebra(lie_type, rank)) +end + + struct BirationalSequence operators::GAP.Obj # TODO Integer operators_vectors::Vector{Vector{Any}} @@ -53,7 +74,7 @@ struct MonomialBasis end struct BasisLieHighestWeightStructure - lie_algebra::LieAlgebra + lie_algebra::LieAlgebraStructure birational_sequence::BirationalSequence highest_weight::Vector{Int} monomial_order::Union{String, Function} @@ -203,15 +224,17 @@ function basis_lie_highest_weight( # steps. Then it starts the recursion and returns the result. # initialization of objects that can be precomputed - lie_algebra_gap, chevalley_basis = create_lie_algebra(type, rank) # lie_algebra of type, rank and its chevalley_basis - lie_algebra = LieAlgebra(type, rank, lie_algebra_gap) + # lie_algebra of type, rank and its chevalley_basis + lie_algebra = LieAlgebraStructure(type, rank) + chevalley_basis = GAP.Globals.ChevalleyBasis(lie_algebra.lie_algebra_gap) + # operators that are represented by our monomials. x_i is connected to operators[i] operators = get_operators(lie_algebra, operators, chevalley_basis) weights = weights_for_operators(lie_algebra.lie_algebra_gap, chevalley_basis[3], operators) # weights of the operators weights = (weight->Int.(weight)).(weights) weights_eps = [w_to_eps(type, rank, w) for w in weights] # other root system - asVec(v) = fromGap(GAP.Globals.ExtRepOfObj(v)) # TODO + asVec(v) = fromGap(GAPWrap.ExtRepOfObj(v)) # TODO birational_sequence = BirationalSequence(operators, [asVec(v) for v in operators], weights, weights_eps) ZZx, _ = PolynomialRing(ZZ, length(operators)) # for our monomials @@ -249,7 +272,7 @@ function sub_simple_refl(word::Vector{Int}, lie_algebra_gap::GAP.Obj)::GAP.Obj return operators end -function get_operators(lie_algebra::LieAlgebra, operators::Union{String, Vector{Int}}, +function get_operators(lie_algebra::LieAlgebraStructure, operators::Union{String, Vector{Int}}, chevalley_basis::GAP.Obj)::GAP.Obj """ handles user input for operators @@ -289,7 +312,7 @@ function get_operators(lie_algebra::LieAlgebra, operators::Union{String, Vector{ end function compute_monomials( - lie_algebra::LieAlgebra, + lie_algebra::LieAlgebraStructure, birational_sequence::BirationalSequence, ZZx::ZZMPolyRing, highest_weight::Vector{Int}, @@ -442,7 +465,7 @@ function add_known_monomials!( end function add_new_monomials!( - lie_algebra::LieAlgebra, + lie_algebra::LieAlgebraStructure, birational_sequence::BirationalSequence, ZZx::ZZMPolyRing, matrices_of_operators::Vector{SMat{ZZRingElem}}, @@ -510,7 +533,7 @@ end function add_by_hand( - lie_algebra::LieAlgebra, + lie_algebra::LieAlgebraStructure, birational_sequence::BirationalSequence, ZZx::ZZMPolyRing, highest_weight::Vector{Int}, @@ -574,33 +597,5 @@ function add_by_hand( return set_mon end -function get_dim_weightspace( - lie_algebra::LieAlgebra, - highest_weight::Vector{Int} - )::Dict{Vector{Int}, Int} - """ - Calculates dictionary with weights as keys and dimension of corresponding weightspace as value. GAP computes the - dimension for all positive weights. The dimension is constant on orbits of the weylgroup, and we can therefore - calculate the dimension of each weightspace. - """ - # calculate dimension for dominant weights with GAP - root_system = GAP.Globals.RootSystem(lie_algebra.lie_algebra_gap) - result = GAP.Globals.DominantCharacter(root_system, GAP.Obj(highest_weight)) - dominant_weights = [map(Int, item) for item in result[1]] - dominant_weights_dim = map(Int, result[2]) - dominant_weights = convert(Vector{Vector{Int}}, dominant_weights) - weightspaces = Dict{Vector{Int}, Int}() - - # calculate dimension for the rest by checking which positive weights lies in the orbit. - for i in 1:length(dominant_weights) - orbit_weights = orbit_weylgroup(lie_algebra, dominant_weights[i]) - dim_weightspace = dominant_weights_dim[i] - for weight in orbit_weights - weightspaces[highest_weight - weight] = dim_weightspace - end - end - return weightspaces -end - end -export BasisLieHighestWeight +export BasisLieHighestWeight \ No newline at end of file diff --git a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl index a25386409c0c..5e6aa83f9ad4 100644 --- a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl +++ b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl @@ -1,12 +1,12 @@ fromGap = Oscar.GAP.gap_to_julia -function create_lie_algebra(type::String, rank::Int)::Tuple{GAP.Obj, GAP.Obj} +function create_lie_algebra(type::String, rank::Int)::GAP.Obj """ Creates the Lie-algebra as a GAP object that gets used for a lot of other computations with GAP """ lie_algebra = GAP.Globals.SimpleLieAlgebra(GAP.Obj(type), rank, GAP.Globals.Rationals) - return lie_algebra, GAP.Globals.ChevalleyBasis(lie_algebra) + return lie_algebra end @@ -25,9 +25,9 @@ function matricesForOperators(lie_algebra::GAP.Obj, highest_weight::Vector{Int}, """ used to create tensorMatricesForOperators """ - M = Oscar.GAP.Globals.HighestWeightModule(lie_algebra, Oscar.GAP.julia_to_gap(highest_weight)) - matrices_of_operators = Oscar.GAP.Globals.List(operators, o -> Oscar.GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o)) - matrices_of_operators = gapReshape.( Oscar.GAP.gap_to_julia(matrices_of_operators)) + M = GAP.Globals.HighestWeightModule(lie_algebra, GAP.julia_to_gap(highest_weight)) + matrices_of_operators = GAP.Globals.List(operators, o -> GAP.Globals.MatrixOfAction(GAPWrap.Basis(M), o)) + matrices_of_operators = gapReshape.(GAP.gap_to_julia(matrices_of_operators)) denominators = map(y->denominator(y[2]), union(union(matrices_of_operators...)...)) common_denominator = lcm(denominators)# // 1 matrices_of_operators = (A->change_base_ring(ZZ, multiply_scalar(A, common_denominator))).(matrices_of_operators) @@ -55,7 +55,7 @@ function weights_for_operators(lie_algebra::GAP.Obj, cartan::GAP.Obj, operators: # TODO delete fromGap. Multiplication of cartan and operators is not regular matrix multiplication cartan = fromGap(cartan, recursive=false) operators = fromGap(operators, recursive=false) - asVec(v) = fromGap(GAP.Globals.ExtRepOfObj(v)) + asVec(v) = fromGap(GAPWrap.ExtRepOfObj(v)) #println(cartan) #println(operators) if any(iszero.(asVec.(operators))) diff --git a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl index 42427711d1af..a9ce7cdddd18 100644 --- a/experimental/BasisLieHighestWeight/src/WeylPolytope.jl +++ b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl @@ -1,4 +1,6 @@ -function orbit_weylgroup(lie_algebra::LieAlgebra, weight_vector::Vector{Int}) + + +function orbit_weylgroup(lie_algebra::LieAlgebraStructure, weight_vector::Vector{Int}) """ operates weyl-group of type type and rank rank on vector weight_vector and returns list of vectors in orbit input and output weights in terms of w_i @@ -9,9 +11,9 @@ function orbit_weylgroup(lie_algebra::LieAlgebra, weight_vector::Vector{Int}) vertices = [] # operate with the weylgroup on weight_vector - GAP.Globals.IsDoneIterator(orbit_iterator) - while !(GAP.Globals.IsDoneIterator(orbit_iterator)) - w = GAP.Globals.NextIterator(orbit_iterator) + GAPWrap.IsDoneIterator(orbit_iterator) + while !(GAPWrap.IsDoneIterator(orbit_iterator)) + w = GAPWrap.NextIterator(orbit_iterator) push!(vertices, Vector{Int}(w)) end @@ -20,6 +22,37 @@ function orbit_weylgroup(lie_algebra::LieAlgebra, weight_vector::Vector{Int}) return vertices end +function get_dim_weightspace( + lie_algebra::LieAlgebraStructure, + highest_weight::Vector{Int} + )::Dict{Vector{Int}, Int} + """ + Calculates dictionary with weights as keys and dimension of corresponding weightspace as value. GAP computes the + dimension for all positive weights. The dimension is constant on orbits of the weylgroup, and we can therefore + calculate the dimension of each weightspace. + """ + # calculate dimension for dominant weights with GAP + root_system = GAP.Globals.RootSystem(lie_algebra.lie_algebra_gap) + result = GAP.Globals.DominantCharacter(root_system, GAP.Obj(highest_weight)) + dominant_weights = [map(Int, item) for item in result[1]] + dominant_weights_dim = map(Int, result[2]) + dominant_weights = convert(Vector{Vector{Int}}, dominant_weights) + weightspaces = Dict{Vector{Int}, Int}() + + # calculate dimension for the rest by checking which positive weights lies in the orbit. + for i in 1:length(dominant_weights) + orbit_weights = orbit_weylgroup(lie_algebra, dominant_weights[i]) + dim_weightspace = dominant_weights_dim[i] + for weight in orbit_weights + weightspaces[highest_weight - weight] = dim_weightspace + end + end + return weightspaces +end + + + + function convert_lattice_points_to_monomials(ZZx, lattice_points_weightspace) return [finish(push_term!(MPolyBuildCtx(ZZx), ZZ(1), convert(Vector{Int}, convert(Vector{Int64}, lattice_point)))) for lattice_point in lattice_points_weightspace] From 5bea68baa21a8986c4a8078944c22d02dd8be170 Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 22:58:35 +0200 Subject: [PATCH 11/15] Custom print functions for the structures --- .../src/BasisLieHighestWeight.jl | 29 +++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 5e02dab632c7..5e90298be0c9 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -15,6 +15,7 @@ using Polymake # TODO Adapt arguments in function outside of BasisLieHighestWeight.jl # TODO Summarize function # TODO Groundup-structure for the special functions +# TODO write methods small when using Polymake, f.e. polyhedron instead of Polyhedron # TODO Export and docstring (?): # basis_lie_highest_weight @@ -58,6 +59,9 @@ function LieAlgebraStructure(lie_type::String, rank::Int) return LieAlgebraStructure(lie_type, rank, create_lie_algebra(lie_type, rank)) end +function Base.show(io::IO, lie_algebra::LieAlgebraStructure) + print(io, "Lie-Algebra of type ", lie_algebra.lie_type, " and rank ", lie_algebra.rank) +end struct BirationalSequence operators::GAP.Obj # TODO Integer @@ -66,6 +70,12 @@ struct BirationalSequence weights_eps::Vector{Vector{Int}} end +function Base.show(io::IO, birational_sequence::BirationalSequence) + println(io, "BirationalSequence") + println(io, "Operators: ", birational_sequence.operators) + print(io, "Weights in w_i:", birational_sequence.weights) +end + struct MonomialBasis set_mon::Set{ZZMPolyRingElem} dimension::Int @@ -73,6 +83,13 @@ struct MonomialBasis # polytope::polytope end +function Base.show(io::IO, monomial_basis::MonomialBasis) + println(io, "MonomialBasis") + println(io, "Dimension: ", monomial_basis.dimension) + println(io, "Generators within semi-group: ", monomial_basis.no_minkowski) + print(io, "Monomials: ", monomial_basis.set_mon) +end + struct BasisLieHighestWeightStructure lie_algebra::LieAlgebraStructure birational_sequence::BirationalSequence @@ -81,6 +98,18 @@ struct BasisLieHighestWeightStructure monomial_basis::MonomialBasis end +function Base.show(io::IO, base::BasisLieHighestWeightStructure) + println(io, base.lie_algebra) + println("") + println(io, base.birational_sequence) + println("") + println(io, "Highest-weight: ", base.highest_weight) + println(io, "Monomial-order: ", base.monomial_order) + println("") + print(io, base.monomial_basis) +end + + include("./VectorSpaceBases.jl") include("./NewMonomial.jl") include("./TensorModels.jl") From 374c86118f8e001ab87ac5ea09fcb633a8740a6f Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 15 Jul 2023 23:25:01 +0200 Subject: [PATCH 12/15] tabs to indents --- experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl | 2 +- experimental/BasisLieHighestWeight/test/runtests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 5e90298be0c9..09bafe9ab3b4 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -627,4 +627,4 @@ function add_by_hand( end end -export BasisLieHighestWeight \ No newline at end of file +export BasisLieHighestWeight diff --git a/experimental/BasisLieHighestWeight/test/runtests.jl b/experimental/BasisLieHighestWeight/test/runtests.jl index 60548b8cbb35..9c3fe7fa09c0 100644 --- a/experimental/BasisLieHighestWeight/test/runtests.jl +++ b/experimental/BasisLieHighestWeight/test/runtests.jl @@ -1,3 +1,3 @@ include("VectorSpaceBases-test.jl") include("NewMonomial-test.jl") -include("BasisLieHighestWeight-test.jl") \ No newline at end of file +include("BasisLieHighestWeight-test.jl") From 0e5a82cc0a430a17738d145af3635274bcc7947a Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sun, 16 Jul 2023 00:03:06 +0200 Subject: [PATCH 13/15] Doctest basis_lie_highest_weight --- .../src/BasisLieHighestWeight.jl | 129 ++++++++++-------- 1 file changed, 69 insertions(+), 60 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 09bafe9ab3b4..3128c0bd071f 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -77,6 +77,7 @@ function Base.show(io::IO, birational_sequence::BirationalSequence) end struct MonomialBasis + ZZx::ZZMPolyRing set_mon::Set{ZZMPolyRingElem} dimension::Int no_minkowski::Set{Vector{Int}} @@ -87,7 +88,8 @@ function Base.show(io::IO, monomial_basis::MonomialBasis) println(io, "MonomialBasis") println(io, "Dimension: ", monomial_basis.dimension) println(io, "Generators within semi-group: ", monomial_basis.no_minkowski) - print(io, "Monomials: ", monomial_basis.set_mon) + print(io, "First 10 Monomials in GRevLex: ", sort(collect(monomial_basis.set_mon), + lt = get_monomial_order_lt("GRevLex", monomial_basis.ZZx))[1:min(end, 10)]) end struct BasisLieHighestWeightStructure @@ -146,65 +148,71 @@ Computes a monomial basis for the highest weight module with highest weight # Examples ```jldoctest -julia> BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1, 1], return_no_minkowski = true, return_operators = true) -(Set(ZZMPolyRingElem[x1*x2, x2, 1, x1*x3, x3^2, x1, x3, x2*x3]), Set([[1, 0], [0, 1]]), GAP: [ v.1, v.2, v.3 ]) - - -julia> BasisLieHighestWeight.basis_lie_highest_weight("A", 3, [2, 2, 3], monomial_order = "Lex") -Set{ZZMPolyRingElem} with 1260 elements: - x3*x5^2*x6^2 - x2*x3*x5^2*x6 - x4^2*x5^2*x6 - x1^2*x3^3*x5 - x2*x3*x4*x5^3*x6^2 - x1*x3*x4*x5^3*x6^2 - x1^2*x3*x4*x6 - x1*x3*x4^3 - x4^2*x5^3*x6^4 - x1*x2*x3*x5^2 - x3^2*x4^4*x5^2*x6 - x2^2*x3*x6^2 - x1*x2^2*x3^2*x5 - x1*x3*x4*x5^2 - x1^2*x2*x6 - x1*x3^2*x4*x5*x6^3 - x1^2*x2*x4*x5^2*x6^2 - x4^4*x5 - x1^2*x2*x3^2*x6 - x1*x3^2*x5^2 - x2*x3*x4*x5^3 - ⋮ - -julia> BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1, 0], operators = [1,2,1]) -Set{ZZMPolyRingElem} with 3 elements: - 1 - x3 - x2*x3 - -julia> BasisLieHighestWeight.basis_lie_highest_weight("C", 3, [1, 1, 1], monomial_order = "Lex") -Set{ZZMPolyRingElem} with 512 elements: - x1*x5*x6*x8 - x6^4 - x3*x4^2*x8 - x3*x4*x6*x7 - x8^3 - x3*x6^2 - x2*x3 - x5*x6^2*x9 - x6*x8^2*x9 - x1*x6*x7 - x5*x6*x9^2 - x6^2*x7^2*x8 - x5*x7*x8 - x4*x6^2*x7*x8^2 - x4^2*x5*x7 - x1*x5^2*x6 - x1*x6*x8 - x3*x4*x5 - x2*x4*x6^2*x7 - x4*x6*x7 - x1*x4*x7*x8^2 - ⋮ +julia> base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1, 1]) +Lie-Algebra of type A and rank 2 + +BirationalSequence +Operators: GAP: [ v.1, v.2, v.3 ] +Weights in w_i:[[2, -1], [-1, 2], [1, 1]] + +Highest-weight: [1, 1] +Monomial-order: GRevLex + +MonomialBasis +Dimension: 8 +Generators within semi-group: Set([[1, 0], [0, 1]]) +First 10 Monomials in GRevLex: ZZMPolyRingElem[1, x3, x2, x1, x3^2, x2*x3, x1*x3, x1*x2] + + + +julia> base = BasisLieHighestWeight.basis_lie_highest_weight("A", 3, [2, 2, 3], monomial_order = "Lex") +Lie-Algebra of type A and rank 3 + +BirationalSequence +Operators: GAP: [ v.1, v.2, v.3, v.4, v.5, v.6 ] +Weights in w_i:[[2, -1, 0], [-1, 2, -1], [0, -1, 2], [1, 1, -1], [-1, 1, 1], [1, 0, 1]] + +Highest-weight: [2, 2, 3] +Monomial-order: Lex + +MonomialBasis +Dimension: 1260 +Generators within semi-group: Set([[0, 1, 0], [1, 0, 0], [0, 0, 1]]) +First 10 Monomials in GRevLex: ZZMPolyRingElem[1, x6, x5, x4, x3, x2, x1, x6^2, x5*x6, x4*x6] + + + +julia> base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1, 0], operators = [1,2,1]) +Lie-Algebra of type A and rank 2 + +BirationalSequence +Operators: GAP: [ v.1, v.2, v.1 ] +Weights in w_i:[[2, -1], [-1, 2], [2, -1]] + +Highest-weight: [1, 0] +Monomial-order: GRevLex + +MonomialBasis +Dimension: 3 +Generators within semi-group: Set([[1, 0]]) +First 10 Monomials in GRevLex: ZZMPolyRingElem[1, x3, x2*x3] + + + +julia> base = BasisLieHighestWeight.basis_lie_highest_weight("C", 3, [1, 1, 1], monomial_order = "Lex") +Lie-Algebra of type C and rank 3 + +BirationalSequence +Operators: GAP: [ v.1, v.2, v.3, v.4, v.5, v.6, v.7, v.8, v.9 ] +Weights in w_i:[[2, -1, 0], [-1, 2, -1], [0, -2, 2], [1, 1, -1], [-1, 0, 1], [1, -1, 1], [-2, 2, 0], [0, 1, 0], [2, 0, 0]] + +Highest-weight: [1, 1, 1] +Monomial-order: Lex + +MonomialBasis +Dimension: 512 +Generators within semi-group: Set([[0, 1, 0], [1, 0, 0], [0, 0, 1], [1, 1, 1], [0, 1, 1]]) +First 10 Monomials in GRevLex: ZZMPolyRingElem[1, x9, x8, x7, x6, x5, x4, x3, x2, x1] ``` """ function basis_lie_highest_weight( @@ -284,6 +292,7 @@ function basis_lie_highest_weight( highest_weight, monomial_order, MonomialBasis( + ZZx, set_mon, length(set_mon), no_minkowski From b6bac29512c2b52cc4d49064f4dc0c2f70ba303a Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 22 Jul 2023 12:31:47 +0200 Subject: [PATCH 14/15] Changed get_monomial_order_lt to accept name of monomial-order as defined in Oscar.jl --- .../src/BasisLieHighestWeight.jl | 35 +++++++++---------- .../src/MonomialOrder.jl | 21 ++++++----- .../test/BasisLieHighestWeight-test.jl | 2 +- 3 files changed, 28 insertions(+), 30 deletions(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 3128c0bd071f..1c284a73851f 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -12,12 +12,10 @@ using Polymake # TODO Change operators from GAP.Obj to Oscar objects # TODO BirationalSequence needs better names for weights and operators # TODO Fix weightnames in general -# TODO Adapt arguments in function outside of BasisLieHighestWeight.jl -# TODO Summarize function # TODO Groundup-structure for the special functions -# TODO write methods small when using Polymake, f.e. polyhedron instead of Polyhedron +# TODO (?) write methods small when using Polymake, f.e. polyhedron instead of Polyhedron -# TODO Export and docstring (?): +# TODO (?) Export and docstring: # basis_lie_highest_weight # get_dim_weightspace # orbit_weylgroup @@ -33,7 +31,7 @@ using Polymake # w_to_eps # eps_to_w -# TODO GAPWrap-wrappers are missing for +# TODO (?) GAPWrap-wrappers are missing for # ChevalleyBasis # DimensionOfHighestWeightModule # SimpleLieAlgebra @@ -81,15 +79,23 @@ struct MonomialBasis set_mon::Set{ZZMPolyRingElem} dimension::Int no_minkowski::Set{Vector{Int}} - # polytope::polytope + polytope::Oscar.Polymake.BigObjectAllocated +end + +function MonomialBasis(ZZx::ZZMPolyRing, set_mon::Set{ZZMPolyRingElem}, no_minkowski::Set{Vector{Int}}) + vertices = degrees.(collect(set_mon)) + vertices_hom = transpose(reduce(hcat, [prepend!(vec, 1) for vec in vertices])) # homogenoues coordinate system + poly = Oscar.Polymake.polytope.Polytope(POINTS=vertices_hom) + return MonomialBasis(ZZx, set_mon, length(set_mon), no_minkowski, poly) end function Base.show(io::IO, monomial_basis::MonomialBasis) println(io, "MonomialBasis") println(io, "Dimension: ", monomial_basis.dimension) println(io, "Generators within semi-group: ", monomial_basis.no_minkowski) - print(io, "First 10 Monomials in GRevLex: ", sort(collect(monomial_basis.set_mon), - lt = get_monomial_order_lt("GRevLex", monomial_basis.ZZx))[1:min(end, 10)]) + println(io, "First 10 Monomials in degrevlex: ", sort(collect(monomial_basis.set_mon), + lt = get_monomial_order_lt("degrevlex", monomial_basis.ZZx))[1:min(end, 10)]) + print(io, "Volume polytope: ", monomial_basis.polytope.VOLUME) end struct BasisLieHighestWeightStructure @@ -163,8 +169,6 @@ Dimension: 8 Generators within semi-group: Set([[1, 0], [0, 1]]) First 10 Monomials in GRevLex: ZZMPolyRingElem[1, x3, x2, x1, x3^2, x2*x3, x1*x3, x1*x2] - - julia> base = BasisLieHighestWeight.basis_lie_highest_weight("A", 3, [2, 2, 3], monomial_order = "Lex") Lie-Algebra of type A and rank 3 @@ -180,8 +184,6 @@ Dimension: 1260 Generators within semi-group: Set([[0, 1, 0], [1, 0, 0], [0, 0, 1]]) First 10 Monomials in GRevLex: ZZMPolyRingElem[1, x6, x5, x4, x3, x2, x1, x6^2, x5*x6, x4*x6] - - julia> base = BasisLieHighestWeight.basis_lie_highest_weight("A", 2, [1, 0], operators = [1,2,1]) Lie-Algebra of type A and rank 2 @@ -197,8 +199,6 @@ Dimension: 3 Generators within semi-group: Set([[1, 0]]) First 10 Monomials in GRevLex: ZZMPolyRingElem[1, x3, x2*x3] - - julia> base = BasisLieHighestWeight.basis_lie_highest_weight("C", 3, [1, 1, 1], monomial_order = "Lex") Lie-Algebra of type C and rank 3 @@ -220,7 +220,7 @@ function basis_lie_highest_weight( rank::Int, highest_weight::Vector{Int}; operators::Union{String, Vector{Int}} = "regular", - monomial_order::Union{String, Function} = "GRevLex", + monomial_order::Union{String, Function} = "degrevlex", cache_size::Int = 0, )::BasisLieHighestWeightStructure @@ -294,8 +294,7 @@ function basis_lie_highest_weight( MonomialBasis( ZZx, set_mon, - length(set_mon), - no_minkowski + no_minkowski, ) ) end @@ -586,7 +585,7 @@ function add_by_hand( """ # initialization # matrices g_i for (g_1^a_1 * ... * g_k^a_k)*v - matrices_of_operators = tensorMatricesForOperators(lie_algebra.lie_algebra_gap, highest_weight, birational_sequence.operators) + matrices_of_operators = tensorMatricesForOperators(lie_algebra.lie_algebra_gap, highest_weight, birational_sequence.operators) space = Dict(0*birational_sequence.weights[1] => SparseVectorSpaceBasis([], [])) # span of basis vectors to keep track of the basis v0 = sparse_row(ZZ, [(1,1)]) # starting vector v # saves the calculated vectors to decrease necessary matrix multiplicatons diff --git a/experimental/BasisLieHighestWeight/src/MonomialOrder.jl b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl index de08134fabd2..8c2998118ef4 100644 --- a/experimental/BasisLieHighestWeight/src/MonomialOrder.jl +++ b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl @@ -2,17 +2,16 @@ function get_monomial_order_lt(monomial_order::Union{String, Function}, ZZx::ZZM """ Returns the desired monomial_order function less than """ - #if isa(monomial_order, Function) - # choosen_monomial_order = monomial_order - x = gens(ZZx) - if monomial_order == "GRevLex" - choosen_monomial_order = degrevlex(x) - elseif monomial_order == "RevLex" - choosen_monomial_order = revlex(x) - elseif monomial_order == "Lex" - choosen_monomial_order = lex(x) + if isa(monomial_order, Function) + choosen_monomial_order = monomial_order else - println("no suitable order picked") + # Ensure that `monomial_order` is a valid function before trying to call it + if isdefined(Main, Symbol(monomial_order)) + x = gens(ZZx) + choosen_monomial_order = eval(Symbol(monomial_order))(x) + else + error("No monomial_order: $monomial_order") + end end - return (mon1, mon2) -> (cmp(choosen_monomial_order, mon1, mon2) == -1) + return (mon1, mon2) -> (cmp(choosen_monomial_order, mon1, mon2) == -1) end diff --git a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl index 5d07421c5b6c..590419b6b3cf 100644 --- a/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl +++ b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl @@ -83,7 +83,7 @@ end end end @testset "Check dimension" begin - @testset "Monomial order $monomial_order" for monomial_order in ("Lex", "RevLex", "GRevLex") + @testset "Monomial order $monomial_order" for monomial_order in ("lex", "revlex", "degrevlex") # the functionality longest-word was temporarily removed because it required coxeter groups from # https://github.com/jmichel7/Gapjm.jl #@testset "Operators $ops" for ops in ("regular", "longest-word") From 35f074a3b22aaeef5973c9f2e566aa173b4575fc Mon Sep 17 00:00:00 2001 From: Ben Wilop Date: Sat, 5 Aug 2023 21:52:06 +0200 Subject: [PATCH 15/15] function body of special basis_lie_highest_weight functions --- .../src/BasisLieHighestWeight.jl | 54 ++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl index 1c284a73851f..dddcfc80aeb7 100644 --- a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -9,13 +9,18 @@ using ..Oscar using ..Oscar: GAPWrap using Polymake +# TODO basis_lie_highest_weight_lustzig +# TODO basis_lie_highest_weight_string +# TODO basis_lie_highest_weight_feigin_fflv +# TODO basis_lie_highest_weight_nZ + # TODO Change operators from GAP.Obj to Oscar objects # TODO BirationalSequence needs better names for weights and operators # TODO Fix weightnames in general # TODO Groundup-structure for the special functions # TODO (?) write methods small when using Polymake, f.e. polyhedron instead of Polyhedron -# TODO (?) Export and docstring: +# TODO (?) Maybe export and docstring: # basis_lie_highest_weight # get_dim_weightspace # orbit_weylgroup @@ -299,6 +304,53 @@ function basis_lie_highest_weight( ) end +function basis_lie_highest_weight_lustzig( + type::String, + rank::Int, + highest_weight::Vector{Int}; + operators::Union{String, Vector{Int}} = "regular", + monomial_order::Union{String, Function} = "oplex", + cache_size::Int = 0, + )::BasisLieHighestWeightStructure + # operators = some sequence of the String / Littelmann-Berenstein-Zelevinsky polytope + return basis_lie_highest_weight(type, rank, highest_weight, monomial_order=monomial_order, cache_size) +end + +function basis_lie_highest_weight_string( + type::String, + rank::Int, + highest_weight::Vector{Int}; + operators::Union{String, Vector{Int}} = "regular", + cache_size::Int = 0, + )::BasisLieHighestWeightStructure + # monomial_order = "oplex" + # operators = some sequence of the String / Littelmann-Berenstein-Zelevinsky polytope + return basis_lie_highest_weight(type, rank, highest_weight, monomial_order=monomial_order, cache_size) +end + +function basis_lie_highest_weight_feigin_fflv( + type::String, + rank::Int, + highest_weight::Vector{Int}; + cache_size::Int = 0, + )::BasisLieHighestWeightStructure + monomial_order = "oplex" + # operators = some sequence of the Feigin-Fourier-Littelmann-Vinberg polytope + return basis_lie_highest_weight(type, rank, highest_weight, monomial_order=monomial_order, cache_size) +end + +function basis_lie_highest_weight_nZ( + type::String, + rank::Int, + highest_weight::Vector{Int}; + operators::Union{String, Vector{Int}} = "regular", + cache_size::Int = 0, + )::BasisLieHighestWeightStructure + # monomial_order = "lex" + # operators = some sequence of the Nakashima-Zelevinsky polytope, same as for _string + return basis_lie_highest_weight(type, rank, highest_weight, monomial_order=monomial_order, cache_size) +end + function sub_simple_refl(word::Vector{Int}, lie_algebra_gap::GAP.Obj)::GAP.Obj """ substitute simple reflections (i,i+1), saved in dec by i, with E_{i,i+1}