8 changed files with 83 additions and 301 deletions
--- a/.gitignore
+++ b/.gitignore
@ -11,5 +11,3 @@ SL*_*
 *.gws
 .*
 tests*
-*.py
-*.pyc
--- a/README.md
+++ b/README.md
@ -1,10 +1,4 @@
-# DEPRECATED!
-
-This repository has not been updated for a while!
-If You are interested in replicating results for [1712.07167](https://arxiv.org/abs/1712.07167) please check [these instruction](https://kalmar.faculty.wmi.amu.edu.pl/post/1712.07176/)
-Also [this notebook](https://nbviewer.jupyter.org/gist/kalmarek/03510181bc1e7c98615e86e1ec580b2a) could be of some help. If everything else fails the [zenodo dataset](https://zenodo.org/record/1133440) should contain the last-resort instructions.
-
-This repository contains some legacy code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
+This repository contains code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).

 # Installing
 To run the code You need `julia-v0.5` (should work on `v0.6`, but with warnings).
--- a/groups/Allgroups.jl
+++ b/groups/Allgroups.jl
@ -4,7 +4,6 @@ using PropertyT
 using AbstractAlgebra
 using Nemo
 using Groups
-using GroupRings

 export PropertyTGroup, SymmetrizedGroup, GAPGroup,
    SpecialLinearGroup,
@ -51,6 +50,5 @@ end
 include("mappingclassgroup.jl")
 include("higman.jl")
 include("caprace.jl")
-include("actions.jl")

 end # of module PropertyTGroups
--- a/groups/actions.jl
+++ b/groups/actions.jl
@ -1,92 +0,0 @@
-function (p::perm)(A::GroupRingElem)
-    RG = parent(A)
-    result = zero(RG, eltype(A.coeffs))
-    
-    for (idx, c) in enumerate(A.coeffs)
-        if c!= zero(eltype(A.coeffs))
-            result[p(RG.basis[idx])] = c
-        end
-    end
-    return result
-end
-
-###############################################################################
-#
-#  Action of WreathProductElems on Nemo.MatElem
-#
-###############################################################################
-
-function matrix_emb(n::DirectProductGroupElem, p::perm)
-    Id = parent(n.elts[1])()
-    elt = diagm([(-1)^(el == Id ? 0 : 1) for el in n.elts])
-    return elt[:, p.d]
-end
-
-function (g::WreathProductElem)(A::MatElem)
-    g_inv = inv(g)
-    G = matrix_emb(g.n, g_inv.p)
-    G_inv = matrix_emb(g_inv.n, g.p)
-    M = parent(A)
-    return M(G)*A*M(G_inv)
-end
-
-import Base.*
-
-doc"""
-    *(x::AbstractAlgebra.MatElem, P::Generic.perm)
-> Apply the pemutation $P$ to the rows of the matrix $x$ and return the result.
-"""
-function *(x::AbstractAlgebra.MatElem, P::Generic.perm)
-   z = similar(x)
-   m = rows(x)
-   n = cols(x)
-   for i = 1:m
-      for j = 1:n
-         z[i, j] = x[i,P[j]]
-      end
-   end
-   return z
-end
-
-function (p::perm)(A::MatElem)
-    length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
-    return p*A*inv(p)
-end
-
-###############################################################################
-#
-#  Action of WreathProductElems on AutGroupElem
-#
-###############################################################################
-
-function AutFG_emb(A::AutGroup, g::WreathProductElem)
-    isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
-    parent(g).P.n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
-    elt = A()
-    Id = parent(g.n.elts[1])()
-    flips = Groups.AutSymbol[Groups.flip_autsymbol(i) for i in 1:length(g.p.d) if g.n.elts[i] != Id]
-    Groups.r_multiply!(elt, flips, reduced=false)
-    Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])
-    return elt
-end
-
-function AutFG_emb(A::AutGroup, p::perm)
-    isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
-    parent(p).n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(p)) into $A")
-    return A(Groups.perm_autsymbol(p))
-end
-
-function (g::WreathProductElem)(a::Groups.Automorphism)
-    A = parent(a)
-    g = AutFG_emb(A,g)
-    res = A()
-    Groups.r_multiply!(res, g.symbols, reduced=false)
-    Groups.r_multiply!(res, a.symbols, reduced=false)
-    Groups.r_multiply!(res, [inv(s) for s in reverse!(g.symbols)])
-    return res
-end
-
-function (p::perm)(a::Groups.Automorphism)
-    g = AutFG_emb(parent(a),p)
-    return g*a*inv(g)
-end
--- a/groups/autfreegroup.jl
+++ b/groups/autfreegroup.jl
@ -19,3 +19,41 @@ end
 function autS(G::SpecialAutomorphismGroup{N}) where N
    return WreathProduct(PermutationGroup(2), PermutationGroup(N))
 end
+
+###############################################################################
+#
+#  Action of WreathProductElems on AutGroupElem
+#
+###############################################################################
+
+function AutFG_emb(A::AutGroup, g::WreathProductElem)
+    isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
+    parent(g).P.n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
+    elt = A()
+    Id = parent(g.n.elts[1])()
+    flips = Groups.AutSymbol[Groups.flip_autsymbol(i) for i in 1:length(g.p.d) if g.n.elts[i] != Id]
+    Groups.r_multiply!(elt, flips, reduced=false)
+    Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])
+    return elt
+end
+
+function AutFG_emb(A::AutGroup, p::perm)
+    isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
+    parent(p).n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
+    return A(Groups.perm_autsymbol(p))
+end
+
+function (g::WreathProductElem)(a::Groups.Automorphism)
+    A = parent(a)
+    g = AutFG_emb(A,g)
+    res = A()
+    Groups.r_multiply!(res, g.symbols, reduced=false)
+    Groups.r_multiply!(res, a.symbols, reduced=false)
+    Groups.r_multiply!(res, [inv(s) for s in reverse!(g.symbols)])
+    return res
+end
+
+function (p::perm)(a::Groups.Automorphism)
+    g = AutFG_emb(parent(a),p)
+    return g*a*inv(g)
+end
--- a/groups/speciallinear.jl
+++ b/groups/speciallinear.jl
@ -60,3 +60,46 @@ end
 function autS(G::SpecialLinearGroup{N}) where N
    return WreathProduct(PermutationGroup(2), PermutationGroup(N))
 end
+
+###############################################################################
+#
+#  Action of WreathProductElems on Nemo.MatElem
+#
+###############################################################################
+
+function matrix_emb(n::DirectProductGroupElem, p::perm)
+    Id = parent(n.elts[1])()
+    elt = diagm([(-1)^(el == Id ? 0 : 1) for el in n.elts])
+    return elt[:, p.d]
+end
+
+function (g::WreathProductElem)(A::MatElem)
+    g_inv = inv(g)
+    G = matrix_emb(g.n, g_inv.p)
+    G_inv = matrix_emb(g_inv.n, g.p)
+    M = parent(A)
+    return M(G)*A*M(G_inv)
+end
+
+import Base.*
+
+doc"""
+    *(x::AbstractAlgebra.MatElem, P::Generic.perm)
+> Apply the pemutation $P$ to the rows of the matrix $x$ and return the result.
+"""
+function *(x::AbstractAlgebra.MatElem, P::Generic.perm)
+   z = similar(x)
+   m = rows(x)
+   n = cols(x)
+   for i = 1:m
+      for j = 1:n
+         z[i, j] = x[i,P[j]]
+      end
+   end
+   return z
+end
+
+function (p::perm)(A::MatElem)
+    length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
+    return p*A*inv(p)
+end
--- a/main.jl
+++ b/main.jl
@ -27,7 +27,7 @@ function Settings(Gr::PropertyTGroup, args, solver)
    S = PropertyTGroups.generatingset(Gr)

    sol = solver
-    ub = get(args,"upper-bound", Inf)
+    ub = get(args,"upper_bound", Inf)
    tol = get(args,"tol", 1e-10)
    ws = get(args, "warmstart", false)

--- a/positivity/check_positivity.jl
+++ b/positivity/check_positivity.jl
@ -1,197 +0,0 @@
-using AbstractAlgebra
-using Groups
-using GroupRings
-using PropertyT
-
-using SCS
-solver(tol, iterations) =
-    SCSSolver(linearsolver=SCS.Direct,
-        eps=tol, max_iters=iterations,
-        alpha=1.95, acceleration_lookback=1)
-
-include("../main.jl")
-
-using PropertyTGroups
-
-args = Dict("SAut" => 5, "upper-bound" => 50.0, "radius" => 2, "nosymmetry"=>false, "tol"=>1e-12, "iterations" =>200000, "warmstart" => true)
-
-Gr = PropertyTGroups.PropertyTGroup(args)
-sett = PropertyT.Settings(Gr, args,
-    solver(args["tol"], args["iterations"]))
-
-@show sett
-
-fullpath = PropertyT.fullpath(sett)
-isdir(fullpath) || mkpath(fullpath)
-# setup_logging(PropertyT.filename(fullpath, :fulllog), :fulllog)
-
-function small_generating_set(RG::GroupRing{AutGroup{N}}, n) where N
-    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
-
-    rmuls = [Groups.rmul_autsymbol(i,j) for (i,j) in indexing]
-    lmuls = [Groups.lmul_autsymbol(i,j) for (i,j) in indexing]
-    gen_set = RG.group.([rmuls; lmuls])
-
-    return [gen_set; inv.(gen_set)]
-end
-
-function computeX(RG::GroupRing{AutGroup{N}}) where N
-    Tn = small_generating_set(RG, N-1)
-
-    ℤ = Int64
-    Δn = length(Tn)*one(RG, ℤ) - RG(Tn, ℤ);
-
-    Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
-
-    @time X = sum(σ(Δn)*sum(τ(Δn) for τ ∈ Alt_N if τ ≠ σ) for σ in Alt_N);
-    return X
-end
-
-function Sq(RG::GroupRing{AutGroup{N}}) where N
-    T2 = small_generating_set(RG, 2)
-    ℤ = Int64
-    Δ₂ = length(T2)*one(RG, ℤ) - RG(T2, ℤ);
-
-    Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
-    elt = sum(σ(Δ₂)^2 for σ in Alt_N)
-    return elt
-end
-
-function Adj(RG::GroupRing{AutGroup{N}}) where N
-    T2 = small_generating_set(RG, 2)
-
-    ℤ = Int64
-    Δ₂ = length(T2)*one(RG, ℤ) - RG(T2, ℤ);
-
-    Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
-
-    adj(σ::perm, τ::perm, i=1, j=2) = Set([σ[i], σ[j]]) ∩ Set([τ[i], τ[j]])
-    adj(σ::perm) = [τ for τ in Alt_N if length(adj(σ, τ)) == 1]
-
-    @time elt = sum(σ(Δ₂)*sum(τ(Δ₂) for τ in adj(σ)) for σ in Alt_N);
-    # return RG(elt.coeffs÷factorial(N-2)^2)
-    return elt
-end
-
-function Op(RG::GroupRing{AutGroup{N}}) where N
-    T2 = small_generating_set(RG, 2)
-
-    ℤ = Int64
-    Δ₂ = length(T2)*one(RG, ℤ) - RG(T2, ℤ);
-
-    Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
-
-    adj(σ::perm, τ::perm, i=1, j=2) = Set([σ[i], σ[j]]) ∩ Set([τ[i], τ[j]])
-    adj(σ::perm) = [τ for τ in Alt_N if length(adj(σ, τ)) == 0]
-
-    @time elt = sum(σ(Δ₂)*sum(τ(Δ₂) for τ in adj(σ)) for σ in Alt_N);
-    # return RG(elt.coeffs÷factorial(N-2)^2)
-    return elt
-end
-
-const ELT_FILE = joinpath(dirname(PropertyT.filename(sett, :Δ)), "SqAdjOp_coeffs.jld")
-const WARMSTART_FILE = PropertyT.filename(sett, :warmstart)
-
-if isfile(PropertyT.filename(sett,:Δ)) && isfile(ELT_FILE) &&
-    isfile(PropertyT.filename(sett, :OrbitData))
-    # cached
-    Δ = PropertyT.loadGRElem(PropertyT.filename(sett,:Δ), sett.G)
-    RG = parent(Δ)
-    orbit_data = load(PropertyT.filename(sett, :OrbitData), "OrbitData")
-    sq_c, adj_c, op_c = load(ELT_FILE, "Sq", "Adj", "Op")
-    # elt = ELT_FILE, sett.G)
-    sq = GroupRingElem(sq_c, RG)
-    adj = GroupRingElem(adj_c, RG)
-    op = GroupRingElem(op_c, RG);
-else
-    info("Compute Laplacian")
-    Δ = PropertyT.Laplacian(sett.S, sett.radius)
-    RG = parent(Δ)
-
-    info("Compute Sq, Adj, Op")
-    @time sq, adj, op = Sq(RG), Adj(RG), Op(RG)
-
-    PropertyT.saveGRElem(PropertyT.filename(sett, :Δ), Δ)
-    save(ELT_FILE, "Sq", sq.coeffs, "Adj", adj.coeffs, "Op", op.coeffs)
-
-    info("Compute OrbitData")
-    if !isfile(PropertyT.filename(sett, :OrbitData))
-        orbit_data = PropertyT.OrbitData(parent(Y), sett.autS)
-        save(PropertyT.filename(sett, :OrbitData), "OrbitData", orbit_data)
-    else
-        orbit_data = load(PropertyT.filename(sett, :OrbitData), "OrbitData")
-    end
-end;
-
-orbit_data = PropertyT.decimate(orbit_data);
-
-elt = adj+2op;
-
-const SOLUTION_FILE = PropertyT.filename(sett, :solution)
-
-if !isfile(SOLUTION_FILE)
-    
-    SDP_problem, varλ, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=sett.upper_bound)
-
-    begin
-        using SCS
-        scs_solver = SCS.SCSSolver(linearsolver=SCS.Direct,
-            eps=sett.tol,
-            max_iters=args["iterations"],
-            alpha=1.95,
-            acceleration_lookback=1)
-
-        JuMP.setsolver(SDP_problem, scs_solver)
-    end
-
-    λ = Ps  = nothing
-    ws = PropertyT.warmstart(sett)
-
-    # using ProgressMeter
-
-    # @showprogress "Running SDP optimization step... " for i in 1:args["repetitions"]
-    while true
-        λ, Ps, ws = PropertyT.solve(PropertyT.filename(sett, :solverlog),
-        SDP_problem, varλ, varP, ws);
-
-        if all((!isnan).(ws[1]))
-            save(WARMSTART_FILE, "warmstart", ws, "λ", λ, "Ps", Ps)
-            save(WARMSTART_FILE[1:end-4]*"_$(now()).jld", "warmstart", ws, "λ", λ, "Ps", Ps)
-        else
-            warn("No valid solution was saved!")
-        end
-    end
-    
-    info("Reconstructing P...")
-    @time P = PropertyT.reconstruct(Ps, orbit_data);
-    save(SOLUTION_FILE, "λ", λ, "P", P)
-end
-
-P, λ = load(SOLUTION_FILE, "P", "λ")
-@show λ;
-
-@time const Q = real(sqrtm(P));
-
-function SOS_residual(eoi::GroupRingElem, Q::Matrix)
-    RG = parent(eoi)
-    @time sos = PropertyT.compute_SOS(RG, Q);
-    return eoi - sos
-end
-
-info("Floating Point arithmetic:")
-EOI = elt - λ*Δ
-b = SOS_residual(EOI, Q)
-@show norm(b, 1);
-
-info("Interval arithmetic:")
-using IntervalArithmetic
-Qint = PropertyT.augIdproj(Q);
-@assert all([zero(eltype(Q)) in sum(view(Qint, :, i)) for i in 1:size(Qint, 2)])
-
-EOI_int = elt - @interval(λ)*Δ;
-Q_int = PropertyT.augIdproj(Q);
-@assert all([zero(eltype(Q)) in sum(view(Q_int, :, i)) for i in 1:size(Q_int, 2)])
-b_int = SOS_residual(EOI_int, Q_int)
-@show norm(b_int, 1);
-
-info("λ is certified to be > ", (@interval(λ) - 2^2*norm(b_int,1)).lo)