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/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..434c63f185c6 --- /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 +BasisLieHighestWeight +``` diff --git a/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl new file mode 100644 index 000000000000..dddcfc80aeb7 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl @@ -0,0 +1,690 @@ +module BasisLieHighestWeight +export LieAlgebra +export basis_lie_highest_weight +export is_fundamental +export get_dim_weightspace +export orbit_weylgroup + +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 (?) Maybe 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 + +# 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 + +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 + operators_vectors::Vector{Vector{Any}} + weights::Vector{Vector{Int}} + 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 + ZZx::ZZMPolyRing + set_mon::Set{ZZMPolyRingElem} + dimension::Int + no_minkowski::Set{Vector{Int}} + 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) + 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 + lie_algebra::LieAlgebraStructure + birational_sequence::BirationalSequence + highest_weight::Vector{Int} + monomial_order::Union{String, Function} + 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") +include("./LieAlgebras.jl") +include("./MonomialOrder.jl") +include("./RootConversion.jl") +include("./WeylPolytope.jl") + +fromGap = Oscar.GAP.gap_to_julia + +@doc """ +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 + +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 + +# Examples +```jldoctest +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( + type::String, + rank::Int, + highest_weight::Vector{Int}; + operators::Union{String, Vector{Int}} = "regular", + monomial_order::Union{String, Function} = "degrevlex", + cache_size::Int = 0, + )::BasisLieHighestWeightStructure + + """ + 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 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(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 + 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)])) + # 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 + set_mon = compute_monomials(lie_algebra, birational_sequence, ZZx, highest_weight, monomial_order_lt, calc_highest_weight, cache_size, no_minkowski) + + # output + return BasisLieHighestWeightStructure( + lie_algebra, + birational_sequence, + highest_weight, + monomial_order, + MonomialBasis( + ZZx, + set_mon, + no_minkowski, + ) + ) +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} + """ + 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(lie_algebra::LieAlgebraStructure, operators::Union{String, Vector{Int}}, + 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 <= 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.lie_algebra_gap) + return operators +end + +function compute_monomials( + lie_algebra::LieAlgebraStructure, + 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. + 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:lie_algebra.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.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(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 + # 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(lie_algebra, birational_sequence, ZZx, lambda_1, + monomial_order_lt, calc_highest_weight, cache_size, + no_minkowski) + 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 + 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(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 + 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 one +end + +function compute_sub_weights(highest_weight::Vector{Int})::Vector{Vector{Int}} + """ + 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 + 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!( + 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{SRow{ZZRingElem}, Vector{Int}}}, + space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.SparseVectorSpaceBasis}, + 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!(gens(ZZx) , mon, birational_sequence.weights, matrices_of_operators, calc_monomials, space, + cache_size) + end + + # check if vec extends the basis + if !haskey(space, weight) + space[weight] = SparseVectorSpaceBasis([], []) + end + add_and_reduce!(space[weight], vec) + end +end + +function add_new_monomials!( + lie_algebra::LieAlgebraStructure, + 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{SRow{ZZRingElem}, Vector{Int}}}, + space::Dict{Vector{Int64}, + Oscar.BasisLieHighestWeight.SparseVectorSpaceBasis}, + 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 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( + 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 + 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!(gens(ZZx), mon, birational_sequence.weights, matrices_of_operators, calc_monomials, space, + cache_size) + end + + # check if vec extends the basis + if !haskey(space, weight) + space[weight] = SparseVectorSpaceBasis([], []) + 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( + lie_algebra::LieAlgebraStructure, + 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.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 + 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) + + # 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, birational_sequence.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) + + # 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!(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!(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 + +end +export BasisLieHighestWeight diff --git a/experimental/BasisLieHighestWeight/src/LieAlgebras.jl b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl new file mode 100644 index 000000000000..5e6aa83f9ad4 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/LieAlgebras.jl @@ -0,0 +1,78 @@ +fromGap = Oscar.GAP.gap_to_julia + + +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 +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}, + operators::GAP.Obj)::Vector{SMat{ZZRingElem}} + """ + used to create tensorMatricesForOperators + """ + 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) + 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 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(GAPWrap.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 + ] +end diff --git a/experimental/BasisLieHighestWeight/src/MonomialOrder.jl b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl new file mode 100644 index 000000000000..8c2998118ef4 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/MonomialOrder.jl @@ -0,0 +1,17 @@ +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 + else + # 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) +end diff --git a/experimental/BasisLieHighestWeight/src/NewMonomial.jl b/experimental/BasisLieHighestWeight/src/NewMonomial.jl new file mode 100644 index 000000000000..31d7e9ddb86d --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/NewMonomial.jl @@ -0,0 +1,103 @@ +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::Union{Vector{SMat{ZZRingElem, Hecke.ZZRingElem_Array_Mod.ZZRingElem_Array}}, 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] + # 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 +end + +function highest_calc_sub_monomial( + x::Vector{ZZMPolyRingElem}, + mon::ZZMPolyRingElem, + 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 + """ + sub_mon = copy(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) + 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::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 + )::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) + #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) + sub_mon_cur *= x[i] + weight += weights[i] + if !haskey(space, weight) + space[weight] = SparseVectorSpaceBasis([], []) + 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..34d2be40a1f8 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/RootConversion.jl @@ -0,0 +1,355 @@ +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 = Matrix{Int}(GAP.Globals.CartanMatrix(R)) + 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..663fa6cd9095 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/VectorSpaceBases.jl @@ -0,0 +1,68 @@ +# manages the linear (in-)dependence of integer vectors +# this file is only of use to basis_lie_highest_weight for the basis vectors + +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 + + # 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} + """ + 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::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 + * => 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) + if pivot[j] != newPivot + continue + end + v = reduce_col(v, basis_vectors[j], newPivot) + 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..a9ce7cdddd18 --- /dev/null +++ b/experimental/BasisLieHighestWeight/src/WeylPolytope.jl @@ -0,0 +1,158 @@ + + +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 + """ + # initialization + 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 = [] + + # operate with the weylgroup on weight_vector + GAPWrap.IsDoneIterator(orbit_iterator) + while !(GAPWrap.IsDoneIterator(orbit_iterator)) + w = GAPWrap.NextIterator(orbit_iterator) + push!(vertices, Vector{Int}(w)) + end + + # return + vertices = convert(Vector{Vector{Int}}, vertices) + 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] +end + +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: + weights: the operator weights in eps_i + weight: lambda - mu + + output: all lattice points with weight weight + """ + 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) + 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 or G + 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/BasisLieHighestWeight-test.jl b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl new file mode 100644 index 000000000000..590419b6b3cf --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/BasisLieHighestWeight-test.jl @@ -0,0 +1,100 @@ +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 "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"])) + 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 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)) + 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", "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") + 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/MBOld.jl b/experimental/BasisLieHighestWeight/test/MBOld.jl new file mode 100644 index 000000000000..9d660e8796db --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/MBOld.jl @@ -0,0 +1,296 @@ + # 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 + +Short = UInt8 + +struct SparseVectorSpaceBasis + A::Vector{SRow{ZZRingElem}} + pivot::Vector{Int} +end + +function normalize(v::SRow{ZZRingElem}) + """ + 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::SparseVectorSpaceBasis, v::SRow{ZZRingElem}) + """ + 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{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``. + +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{SRow{ZZRingElem}}} + 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] => SparseVectorSpaceBasis([], [])) + 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] = SparseVectorSpaceBasis([], []) + 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/NewMonomial-test.jl b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl new file mode 100644 index 000000000000..81d7d680042a --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/NewMonomial-test.jl @@ -0,0 +1,37 @@ +using Oscar +using Test + +include("../src/NewMonomial.jl") + +@testset "Test NewMonomial" begin + 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 diff --git a/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl b/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl new file mode 100644 index 000000000000..f2f9abd8eed8 --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/VectorSpaceBases-test.jl @@ -0,0 +1,36 @@ +using Oscar +using Test + +include("../src/VectorSpaceBases.jl") + +@testset "Test VectorSpaceBases" begin + 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 diff --git a/experimental/BasisLieHighestWeight/test/runtests.jl b/experimental/BasisLieHighestWeight/test/runtests.jl new file mode 100644 index 000000000000..9c3fe7fa09c0 --- /dev/null +++ b/experimental/BasisLieHighestWeight/test/runtests.jl @@ -0,0 +1,3 @@ +include("VectorSpaceBases-test.jl") +include("NewMonomial-test.jl") +include("BasisLieHighestWeight-test.jl")