Skip to content

Commit

Permalink
Code from BasisLieHighestWeight, Pull request oscar-system#2115
Browse files Browse the repository at this point in the history
  • Loading branch information
BenWilop authored and lgoettgens committed Oct 18, 2023
1 parent 006f4c7 commit 56bed3e
Show file tree
Hide file tree
Showing 12 changed files with 1,743 additions and 0 deletions.
3 changes: 3 additions & 0 deletions experimental/BasisLieHighestWeight/docs/doc.main
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
[
"basisLieHighestWeight.md"
]
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
```@meta
CurrentModule = Oscar.BasisLieHighestWeight
```

```@setup oscar
using Oscar.BasisLieHighestWeight
```

```@contents
Pages = ["basisLieHighestWeight.md"]
```

# Monomial bases for Lie algebras
```@docs
basisLieHighestWeight2
```
507 changes: 507 additions & 0 deletions experimental/BasisLieHighestWeight/src/BasisLieHighestWeight.jl

Large diffs are not rendered by default.

54 changes: 54 additions & 0 deletions experimental/BasisLieHighestWeight/src/LieAlgebras.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
using Oscar

fromGap = Oscar.GAP.gap_to_julia


function create_lie_lgebra(type::String, rank::Int)::Tuple{GAP.Obj, GAP.Obj}
"""
Creates the Lie-algebra as a GAP object that gets used for a lot of other computations with GAP
"""
lie_algebra = GAP.Globals.SimpleLieAlgebra(GAP.Obj(type), rank, GAP.Globals.Rationals)
return lie_algebra, GAP.Globals.ChevalleyBasis(lie_algebra)
end


gapReshape(A) = sparse_matrix(QQ, hcat(A...))

# temporary workaround for issue 2128
function multiply_scalar(A::SMat{T}, d) where T
for i in 1:nrows(A)
scale_row!(A, i, T(d))
end
return A
end

function matricesForOperators(lie_algebra::GAP.Obj, highest_weight::Vector{Int},
ops::GAP.Obj)::Vector{SMat{ZZRingElem}}
"""
used to create tensorMatricesForOperators
"""
M = Oscar.GAP.Globals.HighestWeightModule(lie_algebra, Oscar.GAP.julia_to_gap(highest_weight))
matrices_of_operators = Oscar.GAP.Globals.List(ops, o -> Oscar.GAP.Globals.MatrixOfAction(GAP.Globals.Basis(M), o))
matrices_of_operators = gapReshape.( Oscar.GAP.gap_to_julia(matrices_of_operators))
denominators = map(y->denominator(y[2]), union(union(matrices_of_operators...)...))
common_denominator = lcm(denominators)# // 1
matrices_of_operators = (A->change_base_ring(ZZ, multiply_scalar(A, common_denominator))).(matrices_of_operators)
return matrices_of_operators
end


function weights_for_operators(lie_algebra::GAP.Obj, cartan::GAP.Obj, operators::GAP.Obj)::Vector{Vector{Int}}
"""
Calculates the weight wts[i] for each operator ops[i]
"""
cartan = fromGap(cartan, recursive=false)
operators = fromGap(operators, recursive=false)
asVec(v) = fromGap(GAP.Globals.ExtRepOfObj(v))
if any(iszero.(asVec.(operators)))
error("ops should be non-zero")
end
nzi(v) = findfirst(asVec(v) .!= 0)
return [
[asVec(h*v)[nzi(v)] / asVec(v)[nzi(v)] for h in cartan] for v in operators
]
end
18 changes: 18 additions & 0 deletions experimental/BasisLieHighestWeight/src/MonomialOrder.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
function get_monomial_order_lt(monomial_order::Union{String, Function}, ZZx::ZZMPolyRing,
x::Vector{ZZMPolyRingElem})::Function
"""
Returns the desired monomial_order function less than
"""
#if isa(monomial_order, Function)
# choosen_monomial_order = monomial_order
if monomial_order == "GRevLex"
choosen_monomial_order = degrevlex(x)
elseif monomial_order == "RevLex"
choosen_monomial_order = revlex(x)
elseif monomial_order == "Lex"
choosen_monomial_order = lex(x)
else
println("no suitable order picked")
end
return (mon1, mon2) -> (cmp(choosen_monomial_order, mon1, mon2) == -1)
end
101 changes: 101 additions & 0 deletions experimental/BasisLieHighestWeight/src/NewMonomial.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
include("./VectorSpaceBases.jl")
include("./TensorModels.jl")
include("./LieAlgebras.jl")
include("./MonomialOrder.jl")
include("./WeylPolytope.jl")

fromGap = Oscar.GAP.gap_to_julia


function calc_weight(mon::ZZMPolyRingElem, weights::Vector{Vector{Int}})::Vector{Int}
"""
calculates weight associated with monomial mon
"""
degree_mon = degrees(mon)
weight = [0 for i in 1:length(weights[1])]
for i in 1:length(degree_mon)
weight .+= degree_mon[i] * weights[i]
end
return weight
end

function calc_vec(v0::SRow{ZZRingElem}, mon::ZZMPolyRingElem,
matrices_of_operators::Vector{SMat{ZZRingElem}})::SRow{ZZRingElem}
"""
calculates vector associated with monomial mon
"""
vec = v0
degree_mon = degrees(mon)
for i in length(degree_mon):-1:1
for j in 1:degree_mon[i]
vec = mul(vec, transpose(matrices_of_operators[i])) # currently there is no sparse matrix * vector mult
end
end
return vec
end

function highest_calc_sub_monomial(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem,
calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}},
number_of_operators::Int)::ZZMPolyRingElem
"""
returns the key in calc_monomials that can be extended by the least amount of left-operations to mon
"""
sub_mon = copy(mon)

for i in 1:number_of_operators
while is_divisible_by(sub_mon, x[i])
if haskey(calc_monomials, sub_mon)
return sub_mon
else
sub_mon /= x[i]
end
end
end
return sub_mon
end

function calc_new_mon!(x::Vector{ZZMPolyRingElem}, mon::ZZMPolyRingElem, weights::Vector{Vector{Int}},
matrices_of_operators::Vector{SMat{ZZRingElem}}, number_of_operators::Int,
calc_monomials::Dict{ZZMPolyRingElem, Tuple{TVec, Vector{Int}}},
space::Dict{Vector{Int64}, Oscar.BasisLieHighestWeight.VSBasis}, e::Vector{Vector{Int}},
cache_size::Int)::SRow{ZZRingElem}
# calculate vector of mon by extending a previous calculated vector to a
# monom that differs only by left-multiplication, save results in calc_monomials
sub_mon = highest_calc_sub_monomial(x, mon, calc_monomials, number_of_operators)
#println("sub_mon: ", sub_mon)
sub_mon_cur = copy(sub_mon)
(vec, weight) = calc_monomials[sub_mon]
for i in number_of_operators:-1:1
for k in degrees(sub_mon)[i]:(degrees(mon)[i]-1)
sub_mon_cur *= x[i]
weight += weights[i]
if !haskey(space, weight)
space[weight] = nullSpace()
end

vec = mul(vec, transpose(matrices_of_operators[i])) # currently there is no sparse matrix * vector mult
if length(calc_monomials) < cache_size
calc_monomials[sub_mon_cur] = (vec, weight)
end

# check if the extended monomial can be deleted from calculated_monomials, i.e. the other possible
# extensions by left multiplication with some x[i] are already contained
can_be_deleted = true
k = number_of_operators
for l in 1:number_of_operators
if degrees(sub_mon_cur - x[i])[l] != 0
k = l
end
end
for l in 1:k
can_be_deleted = can_be_deleted && haskey(calc_monomials, sub_mon_cur - x[i] + x[l])
end
if can_be_deleted && sub_mon_cur != x[i]
delete!(calc_monomials, sub_mon_cur - x[i])
end
end
end
#println(length(calc_monomials))
return vec
end

Loading

0 comments on commit 56bed3e

Please sign in to comment.