add deprecation warning

everything is handled by run.jl

.gitignore

@@ -10,4 +10,6 @@ SL*_*
 *ipynb*
 *.gws
 .*
-*test*
+tests*
+*.py
+*.pyc

AutFn.jl

@@ -1,69 +0,0 @@
-using ArgParse
-
-###############################################################################
-#
-#  Parsing command line
-#
-###############################################################################
-
-function parse_commandline()
-    s = ArgParseSettings()
-
-    @add_arg_table s begin
-        "--tol"
-            help = "set numerical tolerance for the SDP solver"
-            arg_type = Float64
-            default = 1e-6
-        "--iterations"
-            help = "set maximal number of iterations for the SDP solver"
-            arg_type = Int
-            default = 50000
-        "--upper-bound"
-            help = "Set an upper bound for the spectral gap"
-            arg_type = Float64
-            default = Inf
-        "--cpus"
-            help = "Set number of cpus used by solver (default: auto)"
-            arg_type = Int
-            required = false
-        "--radius"
-            help = "Radius of ball B_r(e,S) to find solution over"
-            arg_type = Int
-            default = 2
-        "--warmstart"
-            help = "Use warmstart.jld as the initial guess for SCS"
-            action = :store_true
-        "--nosymmetry"
-            help = "Don't use symmetries of the Laplacian"
-            action = :store_true
-        "N"
-            help = "Compute for the automorphisms group of the free group on N generators"
-            arg_type = Int
-            required = true
-    end
-    return parse_args(s)
-end
-const PARSEDARGS = parse_commandline()
-
-#=
-Note that the element
-    α(i,j,k) = ϱ(i,j)*ϱ(i,k)*inv(ϱ(i,j))*inv(ϱ(i,k)),
-which surely belongs to ball of radius 4 in Aut(Fₙ) becomes trivial under the representation
-    Aut(Fₙ) → GLₙ(ℤ)⋉ℤⁿ → GL_(n+1)(ℂ).
-Moreover, due to work of Potapchik and Rapinchuk [1] every real representation of Aut(Fₙ) into GLₘ(ℂ) (for m ≤ 2n-2) factors through GLₙ(ℤ)⋉ℤⁿ, so will have the same problem.
-
-We need a different approach: Here we actually compute in (S)Aut(𝔽ₙ)
-=#
-
-include("CPUselect.jl")
-set_parallel_mthread(PARSEDARGS, workers=true)
-
-include("main.jl")
-
-G = PropertyTGroups.SpecialAutomorphismGroup(PARSEDARGS)
-
-if PARSEDARGS["nosymmetry"]
-    main(Standard, G)
-else
-    main(Symmetrize, G)
-end

FPgroup.jl

@@ -1,61 +0,0 @@
-using ArgParse
-
-function parse_commandline()
-   args = ArgParseSettings()
-
-    @add_arg_table args begin
-    "--tol"
-        help = "set numerical tolerance for the SDP solver"
-        arg_type = Float64
-        default = 1e-6
-    "--iterations"
-        help = "set maximal number of iterations for the SDP solver"
-        arg_type = Int
-        default = 50000
-    "--upper-bound"
-        help = "Set an upper bound for the spectral gap"
-        arg_type = Float64
-        default = Inf
-    "--cpus"
-        help = "Set number of cpus used by solver (default: auto)"
-        arg_type = Int
-        required = false
-    "--radius"
-        help = "Radius of ball B_r(e,S) to find solution over"
-        arg_type = Int
-        default = 2
-    "--warmstart"
-        help = "Use warmstart.jl as the initial guess for SCS"
-        action = :store_true
-    "--MCG"
-        help = "Compute for mapping class group of surface of genus N"
-        arg_type = Int
-        required = false
-    "--Higman"
-        help = "Compute for Higman Group"
-        action = :store_true
-    "--Caprace"
-        help = "Compute for Higman Group"
-        action = :store_true
-    end
-
-    return parse_args(args)
-end
-const PARSEDARGS = parse_commandline()
-
-include("CPUselect.jl")
-set_parallel_mthread(PARSEDARGS, workers=false)
-
-include("main.jl")
-
-if PARSEDARGS["Caprace"]
-    G = PropertyTGroups.CapraceGroup()
-elseif PARSEDARGS["Higman"]
-    G = PropertyTGroups.HigmanGroup()
-elseif PARSEDARGS["MCG"] != nothing
-    G = PropertyTGroups.MappingClassGroup(parse_args)
-else
-    throw("You need to specify one of the --Higman, --Caprace, --MCG N")
-end
-
-main(G)

README.md

@@ -1,4 +1,10 @@
-This repository contains code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
+# 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).
 
 # Installing
 To run the code You need `julia-v0.5` (should work on `v0.6`, but with warnings).

SLn.jl

@@ -1,66 +0,0 @@
-using ArgParse
-
-###############################################################################
-#
-#  Parsing command line
-#
-###############################################################################
-
-function parse_commandline()
-    settings = ArgParseSettings()
-
-    @add_arg_table settings begin
-        "--tol"
-            help = "set numerical tolerance for the SDP solver"
-            arg_type = Float64
-            default = 1e-6
-        "--iterations"
-            help = "set maximal number of iterations for the SDP solver"
-            arg_type = Int
-            default = 50000
-        "--upper-bound"
-            help = "Set an upper bound for the spectral gap"
-            arg_type = Float64
-            default = Inf
-        "--cpus"
-            help = "Set number of cpus used by solver"
-            arg_type = Int
-            required = false
-        "--radius"
-            help = "Radius of ball B_r(e,S) to find solution over"
-            arg_type = Int
-            default = 2
-        "--warmstart"
-            help = "Use warmstart.jld as the initial guess for SCS"
-            action = :store_true
-        "--nosymmetry"
-            help = "Don't use symmetries of the Laplacian"
-            action = :store_true
-        "-p"
-            help = "Matrices over field of p-elements (p=0 => over ZZ)"
-            arg_type = Int
-            default = 0
-        "-X"
-            help = "Consider EL(N, ZZ⟨X⟩)"
-            action = :store_true
-        "N"
-            help = "Compute with the group generated by elementary matrices of size n×n"
-            arg_type = Int
-            default = 2
-    end
-    return parse_args(settings)
-end
-const PARSEDARGS = parse_commandline()
-
-include("CPUselect.jl")
-set_parallel_mthread(PARSEDARGS, workers=true)
-
-include("main.jl")
-
-G = PropertyTGroups.SpecialLinearGroup(PARSEDARGS)
-
-if PARSEDARGS["nosymmetry"]
-    main(Standard, G)
-else
-    main(Symmetrize, G)
-end

groups/Allgroups.jl

@@ -1,26 +1,56 @@
 module PropertyTGroups
 
+using PropertyT
 using AbstractAlgebra
 using Nemo
 using Groups
+using GroupRings
 
-export PropertyTGroup, SymmetricGroup, GAPGroup
+export PropertyTGroup, SymmetrizedGroup, GAPGroup,
+    SpecialLinearGroup,
+    SpecialAutomorphismGroup,
+    HigmanGroup,
+    CapraceGroup,
+    MappingClassGroup
+
+export PropertyTGroup
 
 abstract type PropertyTGroup end
 
-abstract type SymmetricGroup <: PropertyTGroup end
+abstract type SymmetrizedGroup <: PropertyTGroup end
 
 abstract type GAPGroup <: PropertyTGroup end
 
+function PropertyTGroup(args)
+    if haskey(args, "SL")
+        G = PropertyTGroups.SpecialLinearGroup(args)
+    elseif haskey(args, "SAut")
+        G = PropertyTGroups.SpecialAutomorphismGroup(args)
+    elseif haskey(args, "MCG")
+        G = PropertyTGroups.MappingClassGroup(args)
+    elseif haskey(args, "Higman")
+        G = PropertyTGroups.HigmanGroup(args)
+    elseif haskey(args, "Caprace")
+        G = PropertyTGroups.CapraceGroup(args)
+    else
+        throw("You must provide one of --SL, --SAut, --MCG, --Higman, --Caprace")
+    end
+    return G
+end
+
 include("autfreegroup.jl")
 include("speciallinear.jl")
 
 Comm(x,y) = x*y*x^-1*y^-1
 
-generatingset(G::GAPGroup) = gens(group(G))
+function generatingset(G::GAPGroup)
+    S = gens(group(G))
+    return unique([S; inv.(S)])
+end
 
 include("mappingclassgroup.jl")
 include("higman.jl")
 include("caprace.jl")
+include("actions.jl")
 
 end # of module PropertyTGroups

groups/actions.jl

@@ -0,0 +1,92 @@
+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

groups/autfreegroup.jl

@@ -1,23 +1,14 @@
-struct SpecialAutomorphismGroup <: SymmetricGroup
+struct SpecialAutomorphismGroup{N} <: SymmetrizedGroup
-    args::Dict{String,Any}
     group::AutGroup
-
-    function SpecialAutomorphismGroup(args::Dict)
-        N = args["N"]
-        return new(args, AutGroup(FreeGroup(N), special=true))
-    end
 end
 
-function name(G::SpecialAutomorphismGroup)
+function SpecialAutomorphismGroup(args::Dict)
-    N = G.args["N"]
+    N = args["SAut"]
-
+    return SpecialAutomorphismGroup{N}(AutGroup(FreeGroup(N), special=true))
-    if G.args["nosymmetry"]
-        return "SAutF$(N)"
-    else
-        return "oSAutF$(N)"
-    end
 end
 
+name(G::SpecialAutomorphismGroup{N}) where N = "SAutF$(N)"
+
 group(G::SpecialAutomorphismGroup) = G.group
 
 function generatingset(G::SpecialAutomorphismGroup)
@@ -25,45 +16,6 @@ function generatingset(G::SpecialAutomorphismGroup)
     return unique([S; inv.(S)])
 end
 
-function autS(G::SpecialAutomorphismGroup)
+function autS(G::SpecialAutomorphismGroup{N}) where N
-    N = G.args["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

groups/caprace.jl

@@ -1,6 +1,4 @@
-struct CapraceGroup <: GAPGroup
+struct CapraceGroup <: GAPGroup end
-    args::Dict{String,Any}
-end
 
 name(G::CapraceGroup) = "CapraceGroup"
 

groups/higman.jl

@@ -1,6 +1,4 @@
-struct HigmanGroup <: GAPGroup
+struct HigmanGroup <: GAPGroup end
-    args::Dict{String,Any}
-end
 
 name(G::HigmanGroup) = "HigmanGroup"
 
@@ -11,14 +9,14 @@ function group(G::HigmanGroup)
    a,b,c,d = gens(higman_group)
 
    relations = [
-   b^-1*Comm(b,a),
+   b*Comm(b,a),
-   c^-1*Comm(c,b),
+   c*Comm(c,b),
-   d^-1*Comm(d,c),
+   d*Comm(d,c),
-   a^-1*Comm(a,d)
+   a*Comm(a,d)
    ];
 
    relations = [relations; [inv(rel) for rel in relations]]
 
-   Groups.add_rels!(higman_group, Dict(rel => HigmanGr() for rel in relations))
+   Groups.add_rels!(higman_group, Dict(rel => higman_group() for rel in relations))
    return higman_group
 end

groups/mappingclassgroup.jl

@@ -1,26 +1,23 @@
-struct MappingClassGroup <: GAPGroup
+struct MappingClassGroup{N} <: GAPGroup end
-    args::Dict{String,Any}
-end
 
-function name(G::MappingClassGroup)
+MappingClassGroup(args::Dict) = MappingClassGroup{args["MCG"]}()
-    N = G.args["N"]
+
-    return "MCG($(N))"
+name(G::MappingClassGroup{N}) where N = "MCG(N)"
-end
+
+function group(G::MappingClassGroup{N}) where N
 
-function group(G::MappingClassGroup)
-    N = G.args["N"]
     if N < 2
         throw("Genus must be at least 2!")
     elseif N == 2
         MCGroup = Groups.FPGroup(["a1","a2","a3","a4","a5"]);
         S = gens(MCGroup)
 
-        N = length(S)
+        n = length(S)
         A = prod(reverse(S))*prod(S)
 
         relations = [
-           [Comm(S[i], S[j]) for i in 1:N for j in 1:N if abs(i-j) > 1]...,
+           [Comm(S[i], S[j]) for i in 1:n for j in 1:n if abs(i-j) > 1]...,
-           [S[i]*S[i+1]*S[i]*inv(S[i+1]*S[i]*S[i+1]) for i in 1:N-1]...,
+           [S[i]*S[i+1]*S[i]*inv(S[i+1]*S[i]*S[i+1]) for i in 1:G.n-1]...,
            (S[1]*S[2]*S[3])^4*inv(S[5])^2,
            Comm(A, S[1]),
            A^2

groups/speciallinear.jl

@@ -1,40 +1,55 @@
-struct SpecialLinearGroup <: SymmetricGroup
+struct SpecialLinearGroup{N} <: SymmetrizedGroup
-    args::Dict{String,Any}
     group::AbstractAlgebra.Group
-
+    p::Int
-    function SpecialLinearGroup(args::Dict)
+    X::Bool
-        n = args["N"]
-        p = args["p"]
-        X = args["X"]
-
-        if p == 0
-            G = MatrixSpace(Nemo.ZZ, n, n)
-        else
-            G = MatrixSpace(FiniteField(p), n, n)
-        end
-        return new(args, G)
-    end
 end
 
-function name(G::SpecialLinearGroup)
+function SpecialLinearGroup(args::Dict)
-    N = G.args["N"]
+    N = args["SL"]
-    p = G.args["p"]
+    p = args["p"]
-    X = G.args["X"]
+    X = args["X"]
 
     if p == 0
-       R = (X ? "Z[x]" : "Z")
+        G = MatrixSpace(Nemo.ZZ, N, N)
     else
-       R = "F$p"
+        R = Nemo.NmodRing(UInt(p))
+        G = MatrixSpace(R, N, N)
     end
-    if G.args["nosymmetry"]
+    return SpecialLinearGroup{N}(G, p, X)
-        return "SL($N,$R)"
+end
+
+function name(G::SpecialLinearGroup{N}) where N
+    if G.p == 0
+       R = (G.X ? "Z[x]" : "Z")
     else
-        return "oSL($N,$R)"
+       R = "F$(G.p)"
     end
+    return SL($(G.N),$R)
 end
 
 group(G::SpecialLinearGroup) = G.group
 
+function generatingset(G::SpecialLinearGroup{N}) where N
+    G.p > 0 && G.X && throw("SL(n, F_p[x]) not implemented")
+    SL = group(G)
+    return generatingset(SL, G.X)
+end
+
+# r is the injectivity radius of
+# SL(n, Z[X]) -> SL(n, Z) induced by X -> 100
+
+function generatingset(SL::MatSpace, X::Bool=false, r=5)
+    n = SL.cols
+    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
+
+    if !X
+        S = [E(idx[1],idx[2],SL) for idx in indexing]
+    else
+        S = [E(i,j,SL,v) for (i,j) in indexing for v in [1, 100*r]]
+    end
+    return unique([S; inv.(S)])
+end
+
 function E(i::Int, j::Int, M::MatSpace, val=one(M.base_ring))
     @assert i≠j
     m = one(M)
@@ -42,52 +57,6 @@ function E(i::Int, j::Int, M::MatSpace, val=one(M.base_ring))
     return m
 end
 
-function generatingset(G::SpecialLinearGroup)
+function autS(G::SpecialLinearGroup{N}) where N
-    n = G.args["N"]
-    p = G.args["p"]
-    X = G.args["X"]
-    SL = group(G)
-    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
-
-    p > 0 && X && throw("SL(n, F_p[x]) not implemented")
-
-    if !X
-        S = [E(idx[1],idx[2],SL) for idx in indexing]
-    else
-        r = G.args["radius"]
-        S = [E(i,j,SL,v) for (i,j) in indexing for v in [1, 100*r]]
-    end
-    return unique([S; inv.(S)])
-end
-
-function autS(G::SpecialLinearGroup)
-    N = G.args["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)
-    res = M(G_inv)
-    Nemo.mul!(res, A, res)
-    return Nemo.mul!(res, M(G), res)
-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

logging.jl

@@ -0,0 +1,58 @@
+using Memento
+
+function setup_logging(filename::String, handlername::Symbol=:log)
+    isdir(dirname(filename)) || mkdir(dirname(filename))
+    logger = Memento.config!("info", fmt="{date}| {msg}")
+    handler = DefaultHandler(filename, DefaultFormatter("{date}| {msg}"))
+    logger.handlers[String(handlername)] = handler
+    return logger
+end
+
+macro logtime(logger, ex)
+    quote
+        local stats = Base.gc_num()
+        local elapsedtime = Base.time_ns()
+        local val = $(esc(ex))
+        elapsedtime = Base.time_ns() - elapsedtime
+        local diff = Base.GC_Diff(Base.gc_num(), stats)
+        local ts = time_string(elapsedtime,
+            diff.allocd,
+            diff.total_time,
+            Base.gc_alloc_count(diff)
+            )
+        $(esc(info))($(esc(logger)), ts)
+        val
+    end
+end
+
+function time_string(elapsedtime, bytes, gctime, allocs)
+    str = @sprintf("%10.6f seconds", elapsedtime/1e9)
+    if bytes != 0 || allocs != 0
+        bytes, mb = Base.prettyprint_getunits(bytes, length(Base._mem_units), Int64(1024))
+        allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000))
+        if ma == 1
+            str*= @sprintf(" (%d%s allocation%s: ", allocs, Base._cnt_units[ma], allocs==1 ? "" : "s")
+        else
+            str*= @sprintf(" (%.2f%s allocations: ", allocs, Base._cnt_units[ma])
+        end
+        if mb == 1
+            str*= @sprintf("%d %s%s", bytes, Base._mem_units[mb], bytes==1 ? "" : "s")
+        else
+            str*= @sprintf("%.3f %s", bytes, Base._mem_units[mb])
+        end
+        if gctime > 0
+            str*= @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime)
+        end
+        str*=")"
+    elseif gctime > 0
+        str*= @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime)
+    end
+    return str
+end
+
+import Base: info, @time
+
+Base.info(x) = info(getlogger(), x)
+macro time(x)
+    return :(@logtime(getlogger(Main), $(esc(x))))
+end

main.jl

@@ -1,127 +1,61 @@
-using AbstractAlgebra
-using Nemo
 using PropertyT
-using Groups
 
-using SCS.SCSSolver
+include("FPGroups_GAP.jl")
-# using Mosek
-# using CSDP
-# using SDPA
 
 include("groups/Allgroups.jl")
 using PropertyTGroups
 
-struct Symmetrize end
+import PropertyT.Settings
-struct Standard end
 
-function summarize(logger, groupdir, iterations, tol, upper_bound, radius, G, S)
+function summarize(sett::PropertyT.Settings)
-    info(logger, "Group:       $groupdir")
+    info("Threads:     $(Threads.nthreads())")
-    info(logger, "Iterations:  $iterations")
+    info("Workers:     $(workers())")
-    info(logger, "Precision:   $tol")
+    info("GroupDir:    $(PropertyT.prepath(sett))")
-    info(logger, "Upper bound: $upper_bound")
+    info(string(sett.G))
-    info(logger, "Radius:      $radius")
+    info("with generating set of size $(length(sett.S))")
-    info(logger, "Threads:     $(Threads.nthreads())")
+
-    info(logger, "Workers:     $(workers())")
+    info("Radius:      $(sett.radius)")
-    info(logger, string(G))
+    info("Precision:   $(sett.tol)")
-    info(logger, "with generating set of size $(length(S))")
+    info("Upper bound: $(sett.upper_bound)")
+    info("Solver:      $(sett.solver)")
 end
 
-function params(Gr::SymmetricGroup)
+function Settings(Gr::PropertyTGroup, args, solver)
-    radius = Gr.args["radius"]
+    r = get(args, "radius", 2)
-    tol = Gr.args["tol"]
+    gr_name = PropertyTGroups.name(Gr)*"_r$r"
-    iterations = Gr.args["iterations"]
-    upper_bound = Gr.args["upper-bound"]
-    warm = Gr.args["warmstart"]
-    N = Gr.args["N"]
-    return radius, tol, iterations, upper_bound, warm, N
-end
-
-function params(Gr::PropertyTGroup)
-    radius = Gr.args["radius"]
-    tol = Gr.args["tol"]
-    iterations = Gr.args["iterations"]
-    upper_bound = Gr.args["upper-bound"]
-    warm = Gr.args["warmstart"]
-    return radius, tol, iterations, upper_bound, warm
-end
-
-main(G::SymmetricGroup) = main(Symmetrize, G)
-
-function main(::Type{Symmetrize}, Gr::SymmetricGroup)
-
-    radius, tol, iterations, upper_bound, warm, N = params(Gr)
-
-    groupdir = "$(PropertyTGroups.name(Gr))_r$radius"
-    isdir(groupdir) || mkdir(groupdir)
-    logger = PropertyT.setup_logging(joinpath(groupdir, "$(upper_bound)"), :fulllog)
-
     G = PropertyTGroups.group(Gr)
     S = PropertyTGroups.generatingset(Gr)
 
-    summarize(logger, groupdir, iterations, tol, upper_bound, radius, G, S)
+    sol = solver
+    ub = get(args,"upper-bound", Inf)
+    tol = get(args,"tol", 1e-10)
+    ws = get(args, "warmstart", false)
 
-    autS = PropertyTGroups.autS(Gr)
+    if get(args, "nosymmetry", false)
-    info(logger, "Symmetrising with $(autS)")
+        return PropertyT.Settings(gr_name, G, S, r, sol, ub, tol, ws)
+    else
+        autS = PropertyTGroups.autS(Gr)
+        return PropertyT.Settings(gr_name, G, S, r, sol, ub, tol, ws, autS)
+    end
+end
 
-    solver = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct, alpha=1.95, acceleration_lookback=1)
+function main(::PropertyTGroup, sett::PropertyT.Settings)
+    isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
 
-    # solver = Mosek.MosekSolver(
+    summarize(sett)
-    #    MSK_DPAR_INTPNT_CO_TOL_REL_GAP=tol,
-    #    MSK_IPAR_INTPNT_MAX_ITERATIONS=iterations,
-    #    QUIET=false)
 
-    # solver = CSDP.CSDPSolver(axtol=tol, atytol=tol, objtol=tol, minstepp=tol*10.0^-1, minstepd=tol*10.0^-1)
-
-    # solver = SDPA.SDPASolver(epsilonStar=tol, epsilonDash=tol)
-
-    sett = Settings(groupdir, N, G, S, autS,
-        radius, solver, upper_bound, tol, warm, logger)
     return PropertyT.check_property_T(sett)
 end
 
-function main(::Type{Standard}, Gr::SymmetricGroup)
+function main(::GAPGroup, sett::PropertyT.Settings)
+    isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
 
-    radius, tol, iterations, upper_bound, warm, _ = params(Gr)
+    summarize(sett)
 
-    groupdir = "$(PropertyTGroups.name(Gr))_r$radius"
+    S = [s for s in sett.S if s.symbols[1].pow == 1]
-    isdir(groupdir) || mkdir(groupdir)
+    relations = [k*inv(v) for (k,v) in sett.G.rels]
-    logger = PropertyT.setup_logging(joinpath(groupdir, "$(upper_bound)"), :fulllog)
 
-    G = PropertyTGroups.group(Gr)
+    prepare_pm_delta(PropertyT.prepath(sett), GAP_groupcode(S, relations), sett.radius)
-    S = PropertyTGroups.generatingset(Gr)
-
-    summarize(logger, groupdir, iterations, tol, upper_bound, radius, G, S)
-
-    solver = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct, alpha=1.95, acceleration_lookback=1)
-
-    return PropertyT.check_property_T(groupdir, S, G(),
-        solver, upper_bound, tol, radius, warm)
 
-end
+    return PropertyT.check_property_T(sett)
-
-
-function main(Gr::GAPGroup)
-
-    include("FPGroups_GAP.jl")
-
-    radius, tol, iterations, upper_bound, warm = params(Gr)
-
-    groupdir = groupname(Gr)
-    isdir(groupdir) || mkdir(groupdir)
-    logger = PropertyT.setup_logging(joinpath(groupdir, "$(upper_bound)"), :fulllog)
-
-    G = PropertyTGroups.group(Gr)
-    S = PropertyTGroups.generatingset(Gr)
-
-    relations = [k*inv(v) for (k,v) in G.rels]
-    prepare_pm_delta(groupdir, GAP_groupcode(S, relations), radius)
-
-    S = unique([S; inv.(S)])
-
-    summarize(logger, groupdir, iterations, tol, upper_bound, radius, G, S)
-
-    solver = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct, alpha=1.95, acceleration_lookback=1)
-
-    return PropertyT.check_property_T(groupdir, S, G(),
-        solver, upper_bound, tol, radius, warm)
 end

positivity/check_positivity.jl

@@ -0,0 +1,197 @@
+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)

run.jl

@@ -0,0 +1,108 @@
+using ArgParse
+
+###############################################################################
+#
+#  Parsing command line
+#
+###############################################################################
+
+function parse_commandline()
+    settings = ArgParseSettings()
+
+    @add_arg_table settings begin
+        "--tol"
+            help = "set numerical tolerance for the SDP solver"
+            arg_type = Float64
+            default = 1e-6
+        "--iterations"
+            help = "set maximal number of iterations for the SDP solver"
+            arg_type = Int
+            default = 50000
+        "--upper-bound"
+            help = "Set an upper bound for the spectral gap"
+            arg_type = Float64
+            default = Inf
+        "--cpus"
+            help = "Set number of cpus used by solver"
+            arg_type = Int
+            required = false
+        "--radius"
+            help = "Radius of ball B_r(e,S) to find solution over"
+            arg_type = Int
+            default = 2
+        "--warmstart"
+            help = "Use warmstart.jld as the initial guess for SCS"
+            action = :store_true
+        "--nosymmetry"
+            help = "Don't use symmetries of the Laplacian"
+            action = :store_true
+
+        "--SL "
+            help = "GROUP: the group generated by elementary matrices of size n by n"
+            arg_type = Int
+            required = false
+        "-p"
+            help = "Matrices over field of p-elements (p=0 => over ZZ) [only with --SL]"
+            arg_type = Int
+            default = 0
+        "-X"
+            help = "Consider EL(N, ZZ⟨X⟩) [only with --SL]"
+            action = :store_true
+
+        "--SAut"
+            help = "GROUP: the automorphisms group of the free group on N generators"
+            arg_type = Int
+            required = false
+
+        "--MCG"
+            help = "GROUP: mapping class group of surface of genus N"
+            arg_type = Int
+            required = false
+
+        "--Higman"
+            help = "GROUP: the Higman Group"
+            action = :store_true
+
+        "--Caprace"
+            help = "GROUP: for Caprace Group"
+            action = :store_true
+    end
+    return parse_args(settings)
+end
+
+const PARSEDARGS = parse_commandline()
+
+set_parallel_mthread(PARSEDARGS, workers=false)
+
+include("CPUselect.jl")
+include("logging.jl")
+include("main.jl")
+
+using SCS.SCSSolver
+# using Mosek
+# using CSDP
+# using SDPA
+
+solver(tol, iterations) =
+    SCSSolver(linearsolver=SCS.Direct,
+        eps=tol, max_iters=iterations,
+        alpha=1.95, acceleration_lookback=1)
+
+# Mosek.MosekSolver(
+#    MSK_DPAR_INTPNT_CO_TOL_REL_GAP=tol,
+#    MSK_IPAR_INTPNT_MAX_ITERATIONS=iterations,
+#    QUIET=false)
+
+# CSDP.CSDPSolver(axtol=tol, atytol=tol, objtol=tol, minstepp=tol*10.0^-1, minstepd=tol*10.0^-1)
+
+# SDPA.SDPASolver(epsilonStar=tol, epsilonDash=tol)
+
+const Gr = PropertyTGroups.PropertyTGroup(PARSEDARGS)
+const sett = PropertyT.Settings(Gr, PARSEDARGS,
+    solver(PARSEDARGS["tol"], PARSEDARGS["iterations"]))
+
+fullpath = PropertyT.fullpath(sett)
+isdir(fullpath) || mkpath(fullpath)
+setup_logging(PropertyT.filename(fullpath, :fulllog), :fulllog)
+
+main(Gr, sett)

runtests.jl

@@ -0,0 +1,222 @@
+using Base.Test
+
+include("main.jl")
+
+testdir = "tests_"*string(now())
+mkdir(testdir)
+include("logging.jl")
+logger=setup_logging(joinpath(testdir, "tests.log"))
+info(testdir)
+
+cd(testdir)
+
+# groupname = name(G)
+# ub = PARSEDARGS["upper-bound"]
+#
+# fullpath = joinpath(groupname, string(ub))
+# isdir(fullpath) || mkpath(fullpath)
+
+separator(n=60) = info("\n"*("\n"*"="^n*"\n"^3)*"\n")
+
+
+function SL_tests(args)
+
+
+    args["SL"] = 2
+    args["p"] = 3
+    G = PropertyTGroup(args)
+    @test main(G) == true
+    separator()
+
+    let args = args
+        args["SL"] = 2
+        args["p"] = 5
+        G = PropertyTGroup(args)
+        @test main(G) == false
+        separator()
+
+        args["warmstart"] = true
+        G = PropertyTGroup(args)
+        @test main(G) == false
+        separator()
+
+        args["upper-bound"] = 0.1
+        G = PropertyTGroup(args)
+        @test main(G) == true
+        separator()
+    end
+
+    args["SL"] = 2
+    args["p"] = 7
+    G = PropertyTGroup(args)
+    @test main(G) == false
+    separator()
+
+    args["SL"] = 3
+    args["p"] = 7
+    G = PropertyTGroup(args)
+    @test main(G) == true
+    separator()
+
+    # begin
+    #     args["iterations"] = 25000
+    #     args["N"] = 3
+    #     args["p"] = 0
+    #     args["upper-bound"] = Inf
+    #
+    #     G = PropertyTGroups.SpecialLinearGroup(args)
+    #     @test main(G) == false
+    #     separator()
+    #
+	#     args["warmstart"] = false
+    #     args["upper-bound"] = 0.27
+    #     G = PropertyTGroups.SpecialLinearGroup(args)
+    #     @test main(G) == false
+    #     separator()
+    #
+    #     args["warmstart"] = true
+    #     G = PropertyTGroups.SpecialLinearGroup(args)
+    #     @test main(G) == true
+    #     separator()
+    # end
+
+    return 0
+end
+
+function SAut_tests(args)
+
+    G = PropertyTGroup(args)
+    @test main(G) == false
+    separator()
+
+    args["warmstart"] = true
+    G = PropertyTGroup(args)
+    @test main(G) == false
+    separator()
+
+    args["upper-bound"] = 0.1
+    G = PropertyTGroup(args)
+    @test main(G) == false
+    separator()
+
+    return 0
+end
+
+@testset "Groups with(out) (T)" begin
+
+    @testset "GAPGroups" begin
+        args = Dict(
+            "Higman" => true,
+            "iterations"=>5000,
+            "tol"=>1e-7,
+            "upper-bound"=>Inf,
+            "cpus"=>2,
+            "radius"=>2,
+            "warmstart"=>false,
+            "nosymmetry"=>true,
+        )
+
+        G = PropertyTGroup(args)
+        @test main(G) == false
+
+        args = Dict(
+            "Caprace" => true,
+            "iterations"=>5000,
+            "tol"=>1e-7,
+            "upper-bound"=>Inf,
+            "cpus"=>2,
+            "radius"=>2,
+            "warmstart"=>false,
+            "nosymmetry"=>true,
+        )
+
+        G = PropertyTGroup(args)
+        @test main(G) == false
+
+        args = Dict(
+            "MCG" => 3,
+            "iterations"=>5000,
+            "tol"=>1e-7,
+            "upper-bound"=>Inf,
+            "cpus"=>2,
+            "radius"=>2,
+            "warmstart"=>false,
+            "nosymmetry"=>true,
+        )
+
+        G = PropertyTGroup(args)
+        @test main(G) == false
+    end
+
+    @testset "SLn's" begin
+        @testset "Non-Symmetrized" begin
+
+            args = Dict(
+                "SL" => 2,
+                "p" => 3,
+                "X" => false,
+                "iterations"=>50000,
+                "tol"=>1e-7,
+                "upper-bound"=>Inf,
+                "cpus"=>2,
+                "radius"=>2,
+                "warmstart"=>false,
+                "nosymmetry"=>true,
+            )
+
+            SL_tests(args)
+        end
+
+        @testset "Symmetrized" begin
+
+            args = Dict(
+                "SL" => 2,
+                "p" => 3,
+                "X" => false,
+                "iterations"=>20000,
+                "tol"=>1e-7,
+                "upper-bound"=>Inf,
+                "cpus"=>2,
+                "radius"=>2,
+                "warmstart"=>false,
+                "nosymmetry"=>false,
+            )
+
+            SL_tests(args)
+        end
+    end
+
+    @testset "SAutF_n's" begin
+
+        @testset "Non-Symmetrized" begin
+
+            args = Dict(
+                "SAut" => 2,
+                "iterations"=>5000,
+                "tol"=>1e-7,
+                "upper-bound"=>Inf,
+                "cpus"=>2,
+                "radius"=>2,
+                "warmstart"=>false,
+                "nosymmetry"=>true,
+            )
+            SAut_tests(args)
+        end
+
+        @testset "Symmetrized" begin
+            args = Dict(
+                "SAut" => 3,
+                "iterations"=>500,
+                "tol"=>1e-7,
+                "upper-bound"=>Inf,
+                "cpus"=>2,
+                "radius"=>2,
+                "warmstart"=>false,
+                "nosymmetry"=>false,
+            )
+            SAut_tests(args)
+        end
+    end
+
+
+end