add graded-by-root-system Adj

src/PropertyT.jl

@@ -19,6 +19,10 @@ include("certify.jl")
 
 include("sqadjop.jl")
 
+include("roots.jl")
+import .Roots
+include("gradings.jl")
+
 include("1712.07167.jl")
 include("1812.03456.jl")
 

src/gradings.jl

@@ -0,0 +1,80 @@
+## something about roots
+
+Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
+    Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
+
+function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
+    if s.symbol === :A
+        return Roots.𝕖(N ÷ 2, s.i) - Roots.𝕖(N ÷ 2, s.j)
+    else#if s.symbol === :B
+        n = N ÷ 2
+        i, j = ifelse(s.i <= n, s.i, s.i - n), ifelse(s.j <= n, s.j, s.j - n)
+        return (-1)^(s.i > s.j) * (Roots.𝕖(n, i) + Roots.𝕖(n, j))
+    end
+end
+
+function Roots.positive(
+    generating_set::AbstractVector{<:MatrixGroups.ElementarySymplectic},
+)
+    r = Roots._positive_direction(Roots.Root(first(generating_set)))
+    pos_gens = [
+        s for s in generating_set if s.val > 0.0 && dot(Roots.Root(s), r) ≥ 0.0
+    ]
+    return pos_gens
+end
+
+grading(s::MatrixGroups.ElementarySymplectic) = Roots.Root(s)
+grading(e::MatrixGroups.ElementaryMatrix) = Roots.Root(e)
+grading(t::Groups.Transvection) = grading(Groups._abelianize(t))
+
+function grading(g::FPGroupElement)
+    if length(word(g)) == 1
+        A = alphabet(parent(g))
+        return grading(A[first(word(g))])
+    else
+        throw("Grading is implemented only for generators")
+    end
+end
+
+_groupby(f, iter::AbstractVector) = _groupby(f.(iter), iter)
+function _groupby(keys::AbstractVector{K}, vals::AbstractVector{V}) where {K,V}
+    @assert length(keys) == length(vals)
+    d = Dict(k => V[] for k in keys)
+    for (k, v) in zip(keys, vals)
+        push!(d[k], v)
+    end
+    return d
+end
+
+function laplacians(RG::StarAlgebra, S, grading)
+    d = _groupby(grading, S)
+    Δs = Dict(α => RG(length(Sα)) - sum(RG(s) for s in Sα) for (α, Sα) in d)
+    return Δs
+end
+
+function Adj(rootsystem::AbstractDict, subtype::Symbol)
+    roots = let W = mapreduce(collect, union, keys(rootsystem))
+        W = union!(W, -1 .* W)
+    end
+
+    return reduce(
+        +,
+        (
+            Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if
+            PropertyT_new.Roots.classify_sub_root_system(
+                roots,
+                first(α),
+                first(β),
+            ) == subtype
+        ),
+        init=zero(first(values(rootsystem))),
+    )
+end
+
+function level(rootsystem, level::Integer)
+    1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}")
+    level == 1 && return Adj(rootsystem, :C₁) # always positive
+    level == 2 && return Adj(rootsystem, :A₁) + Adj(rootsystem, Symbol("C₁×C₁")) + Adj(rootsystem, :C₂) # C₂ is not positive
+    level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁"))
+    level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive
+end

src/roots.jl

@@ -0,0 +1,183 @@
+module Roots
+
+using StaticArrays
+using LinearAlgebra
+
+export Root, isproportional, isorthogonal, ~, ⟂
+
+abstract type AbstractRoot{N,T} end
+
+struct Root{N,T} <: AbstractRoot{N,T}
+    coord::SVector{N,T}
+end
+
+Root(a) = Root(SVector(a...))
+
+function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
+    M == N || return false
+    r.coord == s.coord || return false
+    return true
+end
+
+Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
+
+Base.:+(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord + s.coord)
+Base.:-(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord - s.coord)
+Base.:-(r::Root{N}) where {N} = Root(-r.coord)
+
+Base.:*(a::Number, r::Root) = Root(a * r.coord)
+Base.:*(r::Root, a::Number) = a * r
+
+Base.length(r::AbstractRoot) = norm(r, 2)
+
+LinearAlgebra.norm(r::Root, p::Real=2) = norm(r.coord, p)
+LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
+
+cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
+
+function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
+    N == M || return false
+    val = abs(cos_angle(α, β))
+    return isapprox(val, one(val), atol=eps(one(val)))
+end
+
+function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
+    N == M || return false
+    val = cos_angle(α, β)
+    return isapprox(val, zero(val), atol=eps(one(val)))
+end
+
+function _positive_direction(α::Root{N}) where {N}
+    last = -1 / √2^(N - 1)
+    return Root{N,Float64}(
+        SVector(ntuple(i -> ifelse(i == N, last, (√2)^-i), N)),
+    )
+end
+
+function positive(roots::AbstractVector{<:Root{N}}) where {N}
+    # return those roots for which dot(α, Root([½, ¼, …])) > 0.0
+    pd = _positive_direction(first(roots))
+    return filter(α -> dot(α, pd) > 0.0, roots)
+end
+
+Base.:~(α::AbstractRoot, β::AbstractRoot) = isproportional(α, β)
+⟂(α::AbstractRoot, β::AbstractRoot) = isorthogonal(α, β)
+
+function Base.show(io::IO, r::Root{N}) where {N}
+    print(io, "Root$(r.coord)")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
+    lngth² = sum(x -> x^2, r.coord)
+    l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
+    print(io, "Root in ℝ^$N of length $l\n", r.coord)
+end
+
+E(N, i::Integer) = Root(ntuple(k -> k == i ? 1 : 0, N))
+𝕖(N, i) = E(N, i)
+𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
+
+"""
+    classify_root_system(α, β)
+Return the symbol of smallest system generated by roots `α` and `β`.
+
+The classification is based only on roots length and
+proportionality/orthogonality.
+"""
+function classify_root_system(α::AbstractRoot, β::AbstractRoot)
+    lα, lβ = length(α), length(β)
+    if isproportional(α, β)
+        if lα ≈ lβ ≈ √2
+            return :A₁
+        elseif lα ≈ lβ ≈ 2.0
+            return :C₁
+        else
+            error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
+        end
+    elseif isorthogonal(α, β)
+        if lα ≈ lβ ≈ √2
+            return Symbol("A₁×A₁")
+        elseif lα ≈ lβ ≈ 2.0
+            return Symbol("C₁×C₁")
+        elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
+            return Symbol("A₁×C₁")
+        else
+            error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
+        end
+    else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
+        if lα ≈ lβ ≈ √2
+            return :A₂
+        elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
+            return :C₂
+        else
+            error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
+        end
+    end
+end
+
+function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
+    k = findfirst(v -> isproportional(α, v), Ω)
+    if isnothing(k)
+        error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
+    end
+    return Ω[k]
+end
+
+struct Plane{R<:Root}
+    v1::R
+    v2::R
+    vectors::Vector{R}
+end
+
+Plane(α::R, β::R) where {R<:Root} =
+    Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
+
+function Base.in(r::R, plane::Plane{R}) where {R}
+    return any(isproportional(r, v) for v in plane.vectors)
+end
+
+function classify_sub_root_system(
+    Ω::AbstractVector{<:Root{N}},
+    α::Root{N},
+    β::Root{N},
+) where {N}
+
+    v = proportional_root_from_system(Ω, α)
+    w = proportional_root_from_system(Ω, β)
+
+    subsystem = filter(ω -> ω in Plane(v, w), Ω)
+    @assert length(subsystem) > 0
+    subsystem = positive(union(subsystem, -1 .* subsystem))
+
+    l = length(subsystem)
+    if l == 1
+        x = first(subsystem)
+        return classify_root_system(x, x)
+    elseif l == 2
+        return classify_root_system(subsystem...)
+    elseif l == 3
+        a = classify_root_system(subsystem[1], subsystem[2])
+        b = classify_root_system(subsystem[2], subsystem[3])
+        c = classify_root_system(subsystem[1], subsystem[3])
+
+        if a == b == c # it's only A₂
+            return a
+        end
+
+        C = (:C₂, Symbol("C₁×C₁"))
+        if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c))
+            return :C₂
+        end
+    elseif l == 4
+        for i = 1:l
+            for j = (i+1):l
+                T = classify_root_system(subsystem[i], subsystem[j])
+                T == :C₂ && return :C₂
+            end
+        end
+    end
+    @error "Unknown root subsystem generated by" α β
+    throw("Unknown root system: $subsystem")
+end
+
+end # of module Roots