diff --git a/docs/src/Experimental/basisLieHighestWeight.md b/docs/src/Experimental/basisLieHighestWeight.md index 2b604942f9cb..7fc64f8ce43c 100644 --- a/docs/src/Experimental/basisLieHighestWeight.md +++ b/docs/src/Experimental/basisLieHighestWeight.md @@ -1,9 +1,9 @@ ```@meta -CurrentModule = Oscar +CurrentModule = Oscar.BasisLieHighestWeight ``` ```@setup oscar -using Oscar +using Oscar.BasisLieHighestWeight ``` ```@contents @@ -11,7 +11,6 @@ Pages = ["basisLieHighestWeight.md"] ``` # Monomial bases for Lie algebras - -@docs``` - +```@docs +basisLieHighestWeight2 ``` diff --git a/experimental/basisLieHighestWeight/BasisLieHighestWeight.jl b/experimental/basisLieHighestWeight/BasisLieHighestWeight.jl index df2d0d2850f5..1de6291a3026 100644 --- a/experimental/basisLieHighestWeight/BasisLieHighestWeight.jl +++ b/experimental/basisLieHighestWeight/BasisLieHighestWeight.jl @@ -12,27 +12,27 @@ include("./NewMonomial.jl") #include("./VectorSpaceBases.jl") #--- bekommt gerade noch ZZ, Short und TVEC aus VectorSpaceBases -G = Oscar.GAP.Globals -forGap = Oscar.GAP.julia_to_gap fromGap = Oscar.GAP.gap_to_julia -Markdown.@doc doc""" -basisLieHighestWeight2(t, n, hw; ops, known_monomials, monomial_order, cache_size, parallel:Bool, return_no_minkowski:.Bool, return_ops::Bool) -> Int, Int +@doc Markdown.doc""" + basisLieHighestWeight2(t, n, hw; ops = "regular", known_monomials = [], monomial_order = "GRevLex", cache_size::Int = 1000000, parallel::Bool = false, return_no_minkowski::Bool = false, return_ops::Bool = false) -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``. +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``. + +# Parameters +- `t`: Explain +- `n`: Explain +- `hw`: Explain # Examples ```jldoctest -julia> basisLieHighestWeight2 -output +julia> 1+1 +2 ``` """ - - function basisLieHighestWeight2(t, n, hw; ops = "regular", known_monomials = [], monomial_order = "GRevLex", cache_size::Int = 1000000, parallel::Bool = false, return_no_minkowski::Bool = false, return_ops::Bool = false) - """ - 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``. - """ # The function precomputes objects that are independent of the highest weight and can be used in all recursion steps. Then it starts the recursion and returns the result. # initialization @@ -63,9 +63,9 @@ function sub_simple_refl(word, L, n) """ substitute simple reflections (i,i+1), saved in dec by i, with E_{i,i+1} """ - R = G.RootSystem(L) - CG = fromGap(G.CanonicalGenerators(R)[1], recursive = false) - ops = forGap([CG[i] for i in word], recursive = false) + R = GAP.Globals.RootSystem(L) + CG = fromGap(GAP.Globals.CanonicalGenerators(R)[1], recursive = false) + ops = GAP.Obj([CG[i] for i in word], recursive = false) return ops end @@ -121,7 +121,7 @@ function compute_monomials(t, n, L, hw, ops, wts, wts_eps, monomial_order, calc_ end # calculation required - gapDim = G.DimensionOfHighestWeightModule(L, forGap(hw)) # number of monomials that we need to find, i.e. |M_{hw}|. + gapDim = GAP.Globals.DimensionOfHighestWeightModule(L, GAP.Obj(hw)) # number of monomials that we need to find, i.e. |M_{hw}|. # fundamental weights if is_fundamental(hw) # if hw is a fundamental weight, no partition into smaller summands is possible. This is the basecase of the recursion. push!(no_minkowski, hw) @@ -203,7 +203,7 @@ function add_known_monomials!(weightspace, set_mon_in_weightspace, m, wts, mats, for mon in set_mon_in_weightspace[weight] #vec, wt = calc_new_mon!(mon, m, wts, mats, calc_monomials, space, e, cache_size) d = sz(mats[1]) - v0 = sparse_matrix(ZZ, []) # starting vector v + v0 = sparse_row(ZZ, [(1,1)]) # starting vector v vec = calc_vec(v0, mon, mats) wt = calc_wt(mon, wts) #println("vec:" , vec) @@ -336,9 +336,6 @@ function add_by_hand(t, n, L, hw, ops, wts, wts_eps, monomial_order, gapDim, set #println("size space: ", Int(Base.summarysize(space)) / 2^20) #println("size calc_monomials: ", Int(Base.summarysize(calc_monomials)) / 2^20) add_known_monomials!(weightspace, set_mon_in_weightspace, m, wts, mats, calc_monomials, space, e, cache_size) - if Int(Base.Sys.free_memory()) / 2^20 < 100 - println("-----------------KNOWN------------------------") - end end end @@ -354,9 +351,6 @@ function add_by_hand(t, n, L, hw, ops, wts, wts_eps, monomial_order, gapDim, set #println("size space: ", Int(Base.summarysize(space)) / 2^20) #println("size calc_monomials: ", Int(Base.summarysize(calc_monomials)) / 2^20) add_new_monomials!(t, n, mats, wts, monomial_order, weightspace, wts_eps, set_mon_in_weightspace, calc_monomials, space, e, cache_size, set_mon, m) - if Int(Base.Sys.free_memory()) / 2^20 < 100 - println("-----------------NEW------------------------") - end end end #println("calc_monomials: ", length(calc_monomials)) @@ -370,8 +364,8 @@ function get_dim_weightspace(t, n, L, hw) and we can therefore calculate the dimension of each weightspace """ # calculate dimension for dominant weights with GAP - R = G.RootSystem(L) - W = fromGap(G.DominantCharacter(R, forGap(hw))) + R = GAP.Globals.RootSystem(L) + W = fromGap(GAP.Globals.DominantCharacter(R, GAP.Obj(hw))) dominant_weights = W[1] dominant_weights_dim = W[2] dim_weightspace = [] @@ -389,4 +383,4 @@ end -end \ No newline at end of file +end diff --git a/experimental/basisLieHighestWeight/LieAlgebras.jl b/experimental/basisLieHighestWeight/LieAlgebras.jl index 800e3c85ee29..f02630cce6de 100644 --- a/experimental/basisLieHighestWeight/LieAlgebras.jl +++ b/experimental/basisLieHighestWeight/LieAlgebras.jl @@ -3,8 +3,6 @@ using Oscar #using SparseArrays -G = Oscar.GAP.Globals -forGap = Oscar.GAP.julia_to_gap fromGap = Oscar.GAP.gap_to_julia @@ -12,8 +10,8 @@ function lieAlgebra(t::String, n::Int) """ Creates the Lie-algebra as a GAP object that gets used for a lot other computations with GAP """ - L = G.SimpleLieAlgebra(forGap(t), n, G.Rationals) - return L, G.ChevalleyBasis(L) + L = GAP.Globals.SimpleLieAlgebra(GAP.Obj(t), n, GAP.Globals.Rationals) + return L, GAP.Globals.ChevalleyBasis(L) end @@ -33,9 +31,12 @@ 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)) + #M = GAP.Globals.HighestWeightModule(L, GAP.Obj(hw)) + #mats = GAP.Globals.List(ops, o -> GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o)) + M = Oscar.GAP.Globals.HighestWeightModule(L, Oscar.GAP.julia_to_gap(hw)) + mats = Oscar.GAP.Globals.List(ops, o -> Oscar.GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o)) + #mats = gapReshape.(fromGap(mats)) + mats = gapReshape.( Oscar.GAP.gap_to_julia(mats)) denominators = map(y->denominator(y[2]), union(union(mats...)...)) #d = convert(QQ, lcm(denominators)) d = lcm(denominators)# // 1 @@ -50,7 +51,7 @@ function weightsForOperators(L, cartan, ops) """ cartan = fromGap(cartan, recursive=false) ops = fromGap(ops, recursive=false) - asVec(v) = fromGap(G.ExtRepOfObj(v)) + asVec(v) = fromGap(GAP.Globals.ExtRepOfObj(v)) if any(iszero.(asVec.(ops))) error("ops should be non-zero") end @@ -58,4 +59,4 @@ function weightsForOperators(L, cartan, ops) return [ [asVec(h*v)[nzi(v)] / asVec(v)[nzi(v)] for h in cartan] for v in ops ] -end \ No newline at end of file +end diff --git a/experimental/basisLieHighestWeight/NewMonomial.jl b/experimental/basisLieHighestWeight/NewMonomial.jl index d90d6714833d..6fc26ed1956a 100644 --- a/experimental/basisLieHighestWeight/NewMonomial.jl +++ b/experimental/basisLieHighestWeight/NewMonomial.jl @@ -13,8 +13,6 @@ include("./LieAlgebras.jl") include("./MonomialOrder.jl") include("./WeylPolytope.jl") -G = Oscar.GAP.Globals -forGap = Oscar.GAP.julia_to_gap fromGap = Oscar.GAP.gap_to_julia diff --git a/experimental/basisLieHighestWeight/RootConversion.jl b/experimental/basisLieHighestWeight/RootConversion.jl index c3471cec0a4a..3226ae84e562 100644 --- a/experimental/basisLieHighestWeight/RootConversion.jl +++ b/experimental/basisLieHighestWeight/RootConversion.jl @@ -1,8 +1,6 @@ #using Gapjm using Oscar -G = Oscar.GAP.Globals -forGap = Oscar.GAP.julia_to_gap fromGap = Oscar.GAP.gap_to_julia ############################################ @@ -77,9 +75,9 @@ function alpha_to_w(t, n, weight) end function get_CartanMatrix(t, n) - L = G.SimpleLieAlgebra(forGap(t), n, G.Rationals) - R = G.RootSystem(L) - C_list = fromGap(G.CartanMatrix(R)) + L = GAP.Globals.SimpleLieAlgebra(GAP.Obj(t), n, GAP.Globals.Rationals) + R = GAP.Globals.RootSystem(L) + C_list = fromGap(GAP.Globals.CartanMatrix(R)) C = zeros(n,n) for i in 1:n for j in 1:n @@ -381,4 +379,4 @@ function eps_to_w_A(n, weight) end end return res -end \ No newline at end of file +end diff --git a/experimental/basisLieHighestWeight/TensorModels.jl b/experimental/basisLieHighestWeight/TensorModels.jl index 5199005dc2f4..3665bf9466b8 100644 --- a/experimental/basisLieHighestWeight/TensorModels.jl +++ b/experimental/basisLieHighestWeight/TensorModels.jl @@ -5,9 +5,40 @@ using Oscar # TODO: make the first one a symmetric product, or reduce more generally -spid(n) = spdiagm(0 => [ZZ(1) for _ in 1:n]) -sz(A) = size(A)[1] -tensorProduct(A, B) = kron(A, spid(sz(B))) + kron(spid(sz(A)), B) +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 + #println("ncols(res): ", ncols(res)) + #println("nrows(res): ", nrows(res)) + + 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) #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) @@ -17,9 +48,11 @@ 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) + #println("hw[i]: ", hw[i]) if hw[i] <= 0 continue end @@ -27,8 +60,8 @@ function tensorMatricesForOperators(L, hw, ops) _mats = matricesForOperators(L, wi, ops) _mats = tensorPowers(_mats, hw[i]) mats = mats == [] ? _mats : tensorProducts(mats, _mats) + #println(spdiagm(0 => [ZZ(1) for _ in 1:5])agm) #display(mats) end - return mats end \ No newline at end of file diff --git a/experimental/basisLieHighestWeight/VectorSpaceBases.jl b/experimental/basisLieHighestWeight/VectorSpaceBases.jl index b63c87885960..ccdae590d9a8 100644 --- a/experimental/basisLieHighestWeight/VectorSpaceBases.jl +++ b/experimental/basisLieHighestWeight/VectorSpaceBases.jl @@ -31,7 +31,7 @@ function normalize(v::TVec) end -reduceCol(a, b, i::Int) = a*b[i] - b*a[i] # create zero entry in a +reduceCol(a, b, i::Int) = b[i]*a - a[i]*b # create zero entry in a function addAndReduce!(sp::VSBasis, v::TVec) diff --git a/experimental/basisLieHighestWeight/WeylPolytope.jl b/experimental/basisLieHighestWeight/WeylPolytope.jl index 00d864271374..d9a57d76bd9d 100644 --- a/experimental/basisLieHighestWeight/WeylPolytope.jl +++ b/experimental/basisLieHighestWeight/WeylPolytope.jl @@ -3,8 +3,6 @@ using Oscar include("./RootConversion.jl") -G = Oscar.GAP.Globals -forGap = Oscar.GAP.julia_to_gap fromGap = Oscar.GAP.gap_to_julia ######################### @@ -36,14 +34,14 @@ function orbit_weylgroup(t::String, n::Int, hw) # also possible with polymake orbit_polytope(Vector input_point, Group g), root_system(String type) to save the equations, constant summand missing # initialization L, CH = lieAlgebra(t, n) - W = G.WeylGroup(G.RootSystem(L)) - orb = G.WeylOrbitIterator(W, forGap(hw)) + W = GAP.Globals.WeylGroup(GAP.Globals.RootSystem(L)) + orb = GAP.Globals.WeylOrbitIterator(W, GAP.Obj(hw)) vertices = [] # operate with the weylgroup on hw - G.IsDoneIterator(orb) - while !(G.IsDoneIterator(orb)) - w = G.NextIterator(orb) + GAP.Globals.IsDoneIterator(orb) + while !(GAP.Globals.IsDoneIterator(orb)) + w = GAP.Globals.NextIterator(orb) push!(vertices, fromGap(w)) end @@ -180,4 +178,4 @@ function get_monomials_of_weightspace_Xn(wts, weight) poly = polytope.Polytope(INEQUALITIES=ineq, EQUATIONS=equ) monomials = lattice_points(Polyhedron(poly)) return monomials -end \ No newline at end of file +end diff --git a/experimental/basisLieHighestWeight/main.jl b/experimental/basisLieHighestWeight/main.jl index 9f49cb5b266f..9a1de065dda7 100644 --- a/experimental/basisLieHighestWeight/main.jl +++ b/experimental/basisLieHighestWeight/main.jl @@ -1,2 +1,2 @@ include("BasisLieHighestWeight.jl") -export BasisLieHighestWeight \ No newline at end of file +export BasisLieHighestWeight diff --git a/test/Experimental/basisLieHighestWeight/MB-test.jl b/test/Experimental/basisLieHighestWeight/MB-test.jl index 37af6cc81643..a4e64b720be9 100644 --- a/test/Experimental/basisLieHighestWeight/MB-test.jl +++ b/test/Experimental/basisLieHighestWeight/MB-test.jl @@ -1,7 +1,7 @@ using Oscar using Test using TestSetExtensions -using SparseArrays +#using SparseArrays include("MBOld.jl") @@ -17,7 +17,7 @@ can compute with the weaker version. """ function compare_algorithms(dynkin::Char, n::Int64, lambda::Vector{Int64}) - print("TESTETSETSET", dynkin, n, lambda) + #print("TESTETSETSET", dynkin, n, lambda) dim, m, v = MBOld.basisLieHighestWeight(string(dynkin), n, lambda) # basic algorithm w = BasisLieHighestWeight.basisLieHighestWeight2(string(dynkin), n, lambda) # algorithm that needs to be tested L = G.SimpleLieAlgebra(forGap(string(dynkin)), n, G.Rationals) diff --git a/test/Experimental/basisLieHighestWeight/MBOld.jl b/test/Experimental/basisLieHighestWeight/MBOld.jl index b5f37300bcb0..fd25caaeafa5 100644 --- a/test/Experimental/basisLieHighestWeight/MBOld.jl +++ b/test/Experimental/basisLieHighestWeight/MBOld.jl @@ -7,6 +7,7 @@ export basisLieHighestWeight # use parallel = true for parallel computing using Oscar + G = Oscar.GAP.Globals forGap = Oscar.GAP.julia_to_gap fromGap = Oscar.GAP.gap_to_julia @@ -40,7 +41,7 @@ function normalize(v::TVec) end -reduceCol(a, b, i::Int) = a*b[i] - a[i]*b +reduceCol(a, b, i::Int) = b[i]*a - a[i]*b function addAndReduce!(sp::VSBasis, v::TVec) @@ -109,13 +110,9 @@ function matricesForOperators(L, hw, ops) """ used to create tensorMatricesForOperators """ - println("L: ", L) - println("hw: ", hw) - println("ops: ", ops) M = G.HighestWeightModule(L, forGap(hw)) mats = G.List(ops, o -> G.MatrixOfAction(G.Basis(M), o)) mats = gapReshape.(fromGap(mats)) - println("mats: ", mats) denominators = map(y->denominator(y[2]), union(union(mats...)...)) #d = convert(QQ, lcm(denominators)) d = lcm(denominators)# // 1 @@ -138,30 +135,72 @@ 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 + #println("ncols(res): ", ncols(res)) + #println("nrows(res): ", nrows(res)) + + 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 + # TODO: make the first one a symmetric product, or reduce more generally -spid(n) = spdiagm(0 => [ZZ(1) for _ in 1:n]) -sz(A) = size(A)[1] -tensorProduct(A, B) = kron(A, spid(sz(B))) + kron(spid(sz(A)), B) +#spid(n) = spdiagm(0 => [ZZ(1) for _ in 1:n]) +spid(n) = identity_matrix(SMat, ZZ, n) +sz(A) = nrows(A) #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(L, hw, ops) + """ + Calculates the matrices g_i corresponding to the operator ops[i]. + """ + #println("hw: ", hw) mats = [] for i in 1:length(hw) + #println("hw[i]: ", hw[i]) 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] + #println("Addition cols ", ncols(kron(A, spid(sz(B))) + kron(spid(sz(A)), B))) + #println("Addition rows ",nrows(kron(A, spid(sz(B))) + kron(spid(sz(A)), B))) + end mats = mats == [] ? _mats : tensorProducts(mats, _mats) + #println(spdiagm(0 => [ZZ(1) for _ in 1:5])agm) #display(mats) end - return mats end @@ -259,9 +298,11 @@ function compute(v0, mats, wts::Vector{Vector{Int}}) #--- jedes Monom erhaelt id #--- id(mon) returnt die id des entsprechenden mon. id(mon) = sum((1 << (sum(mon[1:i])+i-1) for i in 1:m-1) , init = 1) # init = 1 for case m = 1 - e = sparse_row(ZZ, [(1,1)]) #--- Einheitsvektoren z.B.: e = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] für m = 3 + #e = [sparse_row(ZZ, [(i,1)]) for i=1:m] #--- Einheitsvektoren z.B.: e = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] für m = 3 + e = [Short.(1:m .== i) for i in 1:m] #maxid(deg) = id(deg.*e[1]) #--- returnt maximale id für Monom vom Grad deg (worst case ist [1, 0, 0], da dann erste Potenz in jedem Summanden ist) - maxid(deg) = id(deg.*[1]) + maxid(deg) = id(deg.*e[1]) + #maxid(deg) = id(sparse_row(ZZ, [(1,deg)])) # TODO #--- Sperrung von Monome @@ -325,7 +366,7 @@ function compute(v0, mats, wts::Vector{Vector{Int}}) space[wt] = nullSpace() end - vec = mul(vector[p], transpose(mats[i])) #--- vec ist neuer potenzieller Kandidat für Basis. + vec = mul(vectors[p], transpose(mats[i])) #--- vec ist neuer potenzieller Kandidat für Basis. vec = addAndReduce!(space[wt], vec) #--- vec ist altes vec ohne Anteil des bereits vorhandenen Erzeugnisses. Wenn 0, dann linear abhängig zu altem Unterraum if vec == 0 #--- wenn vec linear abhängig war, fügen wir vec nicht hinzu continue diff --git a/test/runtests.jl b/test/runtests.jl index d682d42322aa..3bc7b7266608 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -125,7 +125,7 @@ include("Serialization/runtests.jl") include("StraightLinePrograms/runtests.jl") -include("experimental/basisLieHighestWeight/MB-test.jl") +include("Experimental/basisLieHighestWeight/MB-test.jl") @static if compiletimes Base.cumulative_compile_timing(false);