substituted Package SparseArray with Oscar SMat

docs/src/Experimental/basisLieHighestWeight.md

@@ -1,17 +1,16 @@
 ```@meta
-CurrentModule = Oscar
+CurrentModule = Oscar.BasisLieHighestWeight
 ```
 
 ```@setup oscar
-using Oscar
+using Oscar.BasisLieHighestWeight
 ```
 
 ```@contents
 Pages = ["basisLieHighestWeight.md"]
 ```
 
 # Monomial bases for Lie algebras
-
-@docs```
-
+```@docs
+basisLieHighestWeight2
 ```

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
+end

experimental/basisLieHighestWeight/LieAlgebras.jl

@@ -3,17 +3,15 @@
 using Oscar
 #using SparseArrays
 
-G = Oscar.GAP.Globals
-forGap = Oscar.GAP.julia_to_gap
 fromGap = Oscar.GAP.gap_to_julia
 
 
 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,12 +51,12 @@ 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
     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
+end

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
 
 

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
+end

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,18 +48,20 @@ 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])
         mats = mats == [] ? _mats : tensorProducts(mats, _mats)
+        #println(spdiagm(0 => [ZZ(1) for _ in 1:5])agm)
         #display(mats)
     end
-
     return mats
 end

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)

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
+end

experimental/basisLieHighestWeight/main.jl

@@ -1,2 +1,2 @@
 include("BasisLieHighestWeight.jl")
-export BasisLieHighestWeight
+export BasisLieHighestWeight

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)
@@ -27,7 +27,7 @@ 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.basisLieHighestWeight2(string(dynkin), n, lambda, monomial_order=monomial_order, ops=ops) # algorithm that needs to be tested
+    w = BasisLieHighestWeight.basisLieHighestWeight2(string(dynkin), n, lambda, monomial_order=monomial_order) # algorithm that needs to be tested
     L = G.SimpleLieAlgebra(forGap(string(dynkin)), n, G.Rationals)
     gapDim = G.DimensionOfHighestWeightModule(L, forGap(lambda)) # dimension
     @test gapDim == length(w) # check if dimension is correct