add deprecation warning

everything is handled by run.jl

.gitignore

@@ -1,8 +1,15 @@
+Articles
+Higman
+MCG*
+notebooks
+Oldies
+oSAutF*
+oSL*
+SAutF*
+SL*_*
 *ipynb*
 *.gws
-*/*.jld
-*/*.log
-Old*
 .*
-o*
+tests*
-*test*
+*.py
+*.pyc

AutFN.jl

@@ -1,118 +0,0 @@
-using ArgParse
-
-using Nemo
-import SCS.SCSSolver
-using PropertyT
-
-using Groups
-Nemo.setpermstyle(:cycles)
-
-#=
-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₅(ℂ).
-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 Aut(𝔽₄)
-=#
-
-###############################################################################
-#
-#  Generating set
-#
-###############################################################################
-
-function SOutFN_generating_set(N::Int)
-
-   SOutFN = AutGroup(FreeGroup(N), special=true, outer=true)
-   S = gens(SOutFN);
-   S = [S; [inv(s) for s in S]]
-
-   return SOutFN, unique(S)
-end
-
-###############################################################################
-#
-#  Parsing command line
-#
-###############################################################################
-
-function cpuinfo_physicalcores()
-    maxcore = -1
-    for line in eachline("/proc/cpuinfo")
-        if startswith(line, "core id")
-            maxcore = max(maxcore, parse(Int, split(line, ':')[2]))
-        end
-    end
-    maxcore < 0 && error("failure to read core ids from /proc/cpuinfo")
-    return maxcore + 1
-end
-
-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 = 20000
-        "--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
-        "-N"
-            help = "Consider automorphisms of free group on N generators"
-            arg_type = Int
-            default = 2
-    end
-
-    return parse_args(s)
-end
-
-function main()
-
-    parsed_args = parse_commandline()
-    if parsed_args["cpus"] ≠ nothing
-        if parsed_args["cpus"] > cpuinfo_physicalcores()
-            warn("Number of specified cores exceeds the physical core cound. Performance will suffer.")
-        end
-        Blas.set_num_threads(parsed_args["cpus"])
-    end
-
-    N = parsed_args["N"]
-    radius = 2
-    tol = parsed_args["tol"]
-    iterations = parsed_args["iterations"]
-    upper_bound = parsed_args["upper-bound"]
-
-    dirname = "SOutF$(N)_$(upper_bound)_r=$radius"
-
-    logger = PropertyT.setup_logging(dirname)
-
-    info(logger, "Group:       $dirname")
-    info(logger, "Iterations:  $iterations")
-    info(logger, "Precision:   $tol")
-    info(logger, "Upper bound: $upper_bound")
-
-    G, S = SOutFN_generating_set(N)
-    info(logger, G)
-    info(logger, "Symmetric generating set of size $(length(S))")
-    info(logger, S)
-    Id = G()
-
-    solver = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct)
-
-    @time PropertyT.check_property_T(dirname, S, Id, solver, upper_bound, tol, 2)
-    return 0
-end
-
-main()

AutFN_orbit.jl

@@ -1,54 +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-14
-        "--iterations"
-            help = "set maximal number of iterations for the SDP solver (default: 20000)"
-            arg_type = Int
-            default = 60000
-        "--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
-        "-N"
-            help = "Consider automorphisms of free group on N generators"
-            arg_type = Int
-            default = 2
-        "--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
-    end
-
-    return parse_args(settings)
-end
-
-parsed_args = parse_commandline()
-
-include("CPUselect.jl")
-
-set_parallel_mthread(parsed_args, workers=true)
-
-include("SAutFNs.jl")
-include("Orbit.jl")
-
-main(SAutFNs, parsed_args)

CPUselect.jl

@@ -9,16 +9,10 @@ function cpuinfo_physicalcores()
     return maxcore + 1
 end
 
-function set_parallel_mthread(parsed_args; workers=false)
+function set_parallel_mthread(N::Int, workers::Bool)
-
-    if parsed_args["cpus"] == nothing
-        N = cpuinfo_physicalcores()
-    else
-        N = parsed_args["cpus"]
     if N > cpuinfo_physicalcores()
         warn("Number of specified cores exceeds the physical core count. Performance may suffer.")
     end
-    end
 
     if workers
         addprocs(N)
@@ -28,5 +22,14 @@ function set_parallel_mthread(parsed_args; workers=false)
     BLAS.set_num_threads(N)
     info("Using $N threads in BLAS.")
 
-    return N
+end
+
+function set_parallel_mthread(parsed_args::Dict; workers=false)
+    if parsed_args["cpus"] == nothing
+        N = cpuinfo_physicalcores()
+    else
+        N = parsed_args["cpus"]
+    end
+
+    set_parallel_mthread(N, workers)
 end

FPGroups_GAP.jl

@@ -0,0 +1,157 @@
+using JLD
+
+function GAP_code(group_code, dir, R; maxeqns=10_000, infolevel=2)
+    code = """
+LogTo("$(dir)/GAP.log");
+RequirePackage("kbmag");
+SetInfoLevel(InfoRWS, $infolevel);
+
+
+MetricBalls := function(rws, R)
+    local l, basis, sizes, i;
+    l := EnumerateReducedWords(rws, 0, R);;
+    SortBy(l, Length);
+    sizes := [1..R];
+    Apply(sizes, i -> Number(l, w -> Length(w) <= i));
+    return [l, sizes];
+end;;
+
+ProductMatrix := function(rws, basis, len)
+    local result, dict, g, tmpList, t;
+    result := [];
+    dict := NewDictionary(basis[1], true);
+    t := Runtime();
+    for g in [1..Length(basis)] do;
+        AddDictionary(dict, basis[g], g);
+    od;
+    Print("Creating dictionary: \t\t", StringTime(Runtime()-t), "\\n");
+    for g in basis{[1..len]} do;
+        tmpList := List(Inverse(g)*basis{[1..len]}, w->ReducedForm(rws, w));
+        #t := Runtime();
+        tmpList := List(tmpList, x -> LookupDictionary(dict, x));
+        #Print(Runtime()-t, "\\n");
+        Assert(1, ForAll(tmpList, x -> x <> fail));
+        Add(result, tmpList);
+    od;
+    return result;
+end;;
+
+SaveCSV := function(fname, pm)
+   local file, i, j, k;
+   file := OutputTextFile(fname, false);;
+   for i in pm do;
+      k := 1;
+      for j in i do;
+         if k < Length(i) then
+            AppendTo(file, j, ", ");
+         else
+            AppendTo(file, j, "\\n");
+         fi;
+         k := k+1;
+      od;
+   od;
+   CloseStream(file);
+end;;
+
+$group_code
+
+# G:= SimplifiedFpGroup(G);
+RWS := KBMAGRewritingSystem(G);
+# ResetRewritingSystem(RWS);
+O:=OptionsRecordOfKBMAGRewritingSystem(RWS);;
+O.maxeqns := $maxeqns;
+O.maxstates := 1000*$maxeqns;
+#O.maxstoredlen := [100,100];
+
+before := Runtimes();;
+KnuthBendix(RWS);
+after := Runtimes();;
+delta := after.user_time_children - before.user_time_children;;
+Print("Knuth-Bendix completion: \t", StringTime(delta), "\\n");
+
+t := Runtime();
+res := MetricBalls(RWS,$(2*R));;
+Print("Metric-Balls generation: \t", StringTime(Runtime()-t), "\\n");
+B := res[1];; sizes := res[2];;
+Print("Sizes of generated Balls: \t", sizes, "\\n");
+
+t := Runtime();
+pm := ProductMatrix(RWS, B, sizes[$R]);;
+Print("Computing ProductMatrix: \t", StringTime(Runtime()-t), "\\n");
+
+S := EnumerateReducedWords(RWS, 1, 1);
+S := List(S, s -> Position(B,s));
+
+SaveCSV("$(dir)/pm.csv", pm);
+SaveCSV("$(dir)/S.csv", [S]);
+SaveCSV("$(dir)/sizes.csv", [sizes]);
+
+Print("DONE!\\n");
+
+quit;""";
+    return code
+end
+
+function GAP_groupcode(S, rels)
+    F = "FreeGroup("*join(["\"$v\""for v in S], ", ") *");"
+    m = match(r".*(\[.*\])$", string(rels))
+    rels = replace(m.captures[1], " ", "\n")
+    code = """
+F := $F;
+AssignGeneratorVariables(F);;
+relations := $rels;;
+G := F/relations;
+    """
+    return code
+end
+
+function GAP_execute(gap_code, dir)
+   isdir(dir) || mkdir(dir)
+   GAP_file = joinpath(dir, "GAP_code.g")
+   @show dir
+   @show GAP_file;
+
+   open(GAP_file, "w") do io
+      write(io, gap_code)
+   end
+   run(pipeline(`cat $(GAP_file)`, `gap -q`))
+end
+
+function prepare_pm_delta_csv(name, group_code, R; maxeqns=10_000, infolevel=2)
+   info("Preparing multiplication table using GAP (via kbmag)")
+   gap_code = GAP_code(group_code, name, R, maxeqns=maxeqns, infolevel=infolevel)
+   GAP_execute(gap_code, name)
+end
+
+function prepare_pm_delta(name, group_code, R; maxeqns=100_000, infolevel=2)
+
+   pm_fname = joinpath(name, "pm.csv")
+   S_fname = joinpath(name, "S.csv")
+   sizes_fname = joinpath(name, "sizes.csv")
+   delta_fname = joinpath(name, "delta.jld")
+
+   if !isfile(pm_fname) || !isfile(S_fname) || !isfile(sizes_fname)
+      prepare_pm_delta_csv(name, group_code, R, maxeqns=maxeqns, infolevel=infolevel)
+   end
+
+   if isfile(sizes_fname)
+      sizes = readcsv(sizes_fname, Int)[1,:]
+      if 2R > length(sizes)
+         prepare_pm_delta_csv(name, group_code, R, maxeqns=maxeqns, infolevel=infolevel)
+      end
+   end
+
+   pm = readcsv(pm_fname, Int)
+   S = readcsv(S_fname, Int)[1,:]
+   sizes = readcsv(sizes_fname, Int)[1,:]
+
+   Δ = spzeros(sizes[2R])
+   Δ[S] .= -1
+   Δ[1] = length(S)
+
+   pm = pm[1:sizes[R], 1:sizes[R]]
+
+   save(joinpath(name, "pm.jld"), "pm", pm)
+   save(joinpath(name, "delta.jld"), "Δ", Δ)
+
+end

Orbit.jl

@@ -1,53 +0,0 @@
-using SCS.SCSSolver
-# using Mosek
-# using CSDP
-# using SDPA
-
-using Nemo
-
-using PropertyT
-using Groups
-
-function main(GROUP, parsed_args)
-
-   radius = parsed_args["radius"]
-   tol = parsed_args["tol"]
-   iterations = parsed_args["iterations"]
-   upper_bound = parsed_args["upper-bound"]
-   warm = parsed_args["warmstart"]
-
-   name, N = GROUP.groupname(parsed_args)
-   G, S = GROUP.generatingset(parsed_args)
-   autS = GROUP.autS(parsed_args)
-
-   name = "$(name)_r$radius"
-   isdir(name) || mkdir(name)
-
-   logger = PropertyT.setup_logging(joinpath(name, "$(upper_bound)"))
-
-   info(logger, "Group:       $name")
-   info(logger, "Iterations:  $iterations")
-   info(logger, "Precision:   $tol")
-   info(logger, "Upper bound: $upper_bound")
-   info(logger, "Threads:     $(Threads.nthreads())")
-   info(logger, "Workers:     $(workers())")
-   info(logger, G)
-   info(logger, "Symmetric generating set of size $(length(S))")
-   # info(logger, S)
-
-   solver = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct, alpha=1.9, acceleration_lookback=1)
-
-   # solver = Mosek.MosekSolver(
-   #    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(name, N, G, S, autS, radius, solver, upper_bound, tol, warm)
-
-   PropertyT.check_property_T(sett)
-end

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).

SAutFNs.jl

@@ -1,80 +0,0 @@
-module SAutFNs
-
-using Nemo
-using Groups
-
-if VERSION >= v"0.6.0"
-   import Nemo.Generic.perm
-end
-
-###############################################################################
-#
-#  Generating set
-#
-###############################################################################
-
-function generatingset(N::Int)
-
-   SAutFN = AutGroup(FreeGroup(N), special=true)
-   S = gens(SAutFN);
-   S = [S; [inv(s) for s in S]]
-
-   return SAutFN, unique(S)
-end
-
-function generatingset(parsed_args)
-   N = parsed_args["N"]
-   return generatingset(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")
-    powers = [(elt == parent(elt)() ? 0: 1) for elt in g.n.elts]
-    elt = reduce(*, [A(Groups.flip_autsymbol(i))^pow for (i,pow) in enumerate(powers)])
-    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::AutGroupElem)
-    g = AutFG_emb(parent(a),g)
-    return g*a*g^-1
-end
-
-function (p::perm)(a::AutGroupElem)
-    g = AutFG_emb(parent(a),p)
-    return g*a*g^-1
-end
-
-autS(N::Int) = WreathProduct(PermutationGroup(2), PermutationGroup(N))
-
-function autS(parsed_args)
-   return autS(parsed_args["N"])
-end
-
-###############################################################################
-#
-#  Misc
-#
-###############################################################################
-
-function groupname(parsed_args)
-   N = parsed_args["N"]
-   return groupname(N), N
-end
-
-groupname(N::Int) = "oSAutF$(N)"
-
-end # of module SAutFNs

SL.jl

@@ -1,151 +0,0 @@
-using ArgParse
-
-using Nemo
-import SCS.SCSSolver
-using PropertyT
-
-
-###############################################################################
-#
-#  Generating set
-#
-###############################################################################
-
-function E(i::Int, j::Int, M::MatSpace, val=one(M.base_ring))
-    @assert i≠j
-    m = one(M)
-    m[i,j] = val
-    return m
-end
-
-function SLsize(n,p)
-    result = BigInt(1)
-    for k in 0:n-1
-        result *= p^n - p^k
-    end
-    return div(result, p-1)
-end
-
-function SL_generatingset(n::Int, X::Bool=false)
-   indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
-   G = MatrixSpace(ZZ, n, n)
-   if X
-      S = [E(i,j,G,v) for (i,j) in indexing for v in [1, 100]]
-   else
-      S = [E(i,j,G,v) for (i,j) in indexing for v in [1]]
-   end
-   S = vcat(S, [inv(x) for x in S])
-   return G, unique(S)
-end
-
-function SL_generatingset(n::Int, p::Int, X::Bool=false)
-    p == 0 && return SL_generatingset(n, X)
-    (p > 1 && n > 1) || throw("Both n and p should be positive integers!")
-    info("Size(SL($n,$p)) = $(SLsize(n,p))")
-    F = ResidueRing(ZZ, p)
-    G = MatrixSpace(F, n, n)
-    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
-    S = [E(i, j, G) for (i,j) in indexing]
-    S = vcat(S, [inv(x) for x in S])
-    return G, unique(S)
-end
-
-###############################################################################
-#
-#  Parsing command line
-#
-###############################################################################
-
-function cpuinfo_physicalcores()
-    maxcore = -1
-    for line in eachline("/proc/cpuinfo")
-        if startswith(line, "core id")
-            maxcore = max(maxcore, parse(Int, split(line, ':')[2]))
-        end
-    end
-    maxcore < 0 && error("failure to read core ids from /proc/cpuinfo")
-    return maxcore + 1
-end
-
-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 (default: 20000)"
-            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
-        "-N"
-            help = "Consider elementary matrices EL(N)"
-            arg_type = Int
-            default = 2
-        "-p"
-            help = "Matrices over field of p-elements (p=0 => over ZZ)"
-            arg_type = Int
-            default = 0
-        "--radius"
-            help = "Radius of ball B_r(e,S) to find solution over"
-            arg_type = Int
-            default = 2
-    end
-
-    return parse_args(settings)
-end
-
-function main()
-
-    parsed_args = parse_commandline()
-    if parsed_args["cpus"] ≠ nothing
-       if parsed_args["cpus"] > cpuinfo_physicalcores()
-           warn("Number of specified cores exceeds the physical core cound. Performance will suffer.")
-        end
-        BLAS.set_num_threads(parsed_args["cpus"])
-    end
-
-    N = parsed_args["N"]
-    p = parsed_args["p"]
-
-    if p == 0
-        dirname = "SL$(N)Z"
-    else
-        dirname = "SL$(N)_$p"
-    end
-
-    radius = parsed_args["radius"]
-    tol = parsed_args["tol"]
-    iterations = parsed_args["iterations"]
-    upper_bound = parsed_args["upper-bound"]
-
-    dirname = "$(dirname)_$(upper_bound)_r=$radius"
-
-    logger = PropertyT.setup_logging(dirname)
-
-    info(logger, "Group:       $dirname")
-    info(logger, "Iterations:  $iterations")
-    info(logger, "Precision:   $tol")
-    info(logger, "Upper bound: $upper_bound")
-
-    G, S = SL_generatingset(N, p)
-    info(logger, G)
-    info(logger, "Symmetric generating set of size $(length(S))")
-    Id = one(G)
-
-    solver = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct)
-
-    @time PropertyT.check_property_T(dirname, S, Id, solver, upper_bound, tol, radius)
-    return 0
-end
-
-main()

SLNs.jl

@@ -1,124 +0,0 @@
-module SLNs
-
-using Nemo
-using Groups
-
-if VERSION >= v"0.6.0"
-   import Nemo.Generic.perm
-end
-
-###############################################################################
-#
-#  Generating set
-#
-###############################################################################
-
-function E(i::Int, j::Int, M::MatSpace, val=one(M.base_ring))
-    @assert i≠j
-    m = one(M)
-    m[i,j] = val
-    return m
-end
-
-function SLsize(n,p)
-    result = BigInt(1)
-    for k in 0:n-1
-        result *= p^n - p^k
-    end
-    return div(result, p-1)
-end
-
-function generatingset(n::Int, X::Bool=false)
-   indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
-   G = MatrixSpace(ZZ, n, n)
-   if X
-      S = [E(i,j,G,v) for (i,j) in indexing for v in [1, 100]]
-   else
-      S = [E(i,j,G,v) for (i,j) in indexing for v in [1]]
-   end
-   S = vcat(S, [inv(x) for x in S])
-   return G, unique(S)
-end
-
-function generatingset(n::Int, p::Int, X::Bool=false)
-    (p > 1 && n > 1) || throw("Both n and p should be positive integers!")
-    info("Size(SL($n,$p)) = $(SLsize(n,p))")
-    F = ResidueRing(ZZ, p)
-    G = MatrixSpace(F, n, n)
-    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
-    S = [E(i, j, G) for (i,j) in indexing]
-    S = vcat(S, [inv(x) for x in S])
-    return G, unique(S)
-end
-
-function generatingset(parsed_args)
-   N = parsed_args["N"]
-   p = parsed_args["p"]
-   X = parsed_args["X"]
-
-   if p == 0
-      return generatingset(N, X)
-   else
-      return generatingset(N, p, X)
-   end
-end
-
-###############################################################################
-#
-#  Action of WreathProductElems on Nemo.MatElem
-#
-###############################################################################
-
-function matrix_emb(MM::MatSpace, g::WreathProductElem)
-    parent(g).P.n == MM.cols == MM.rows || throw("No natural embedding of $(parent(g)) in ")
-    powers = [(elt == parent(elt)() ? 0: 1) for elt in g.n.elts]
-    elt = diagm([(-1)^(elt == parent(elt)() ? 0: 1) for elt in g.n.elts])
-    return MM(elt)*MM(Nemo.matrix_repr(g.p)')
-end
-
-function (g::WreathProductElem)(A::MatElem)
-    G = matrix_emb(parent(A), g)
-    inv_G = matrix_emb(parent(A), inv(g))
-    return G*A*inv_G
-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))")
-    R = parent(A)
-    return p*A*inv(p)
-end
-
-autS(N::Int) = WreathProduct(PermutationGroup(2), PermutationGroup(N))
-
-function autS(parsed_args)
-   return autS(parsed_args["N"])
-end
-
-###############################################################################
-#
-#  Misc
-#
-###############################################################################
-
-function groupname(parsed_args)
-   N = parsed_args["N"]
-   p = parsed_args["p"]
-   X = parsed_args["X"]
-
-   return groupname(N, p, X), N
-end
-
-function groupname(N, p=0, X=false)
-    if p == 0
-       if X
-          name = "oSL$(N)Z⟨X⟩"
-       else
-          name = "oSL$(N)Z"
-       end
-    else
-       name = "oSL$(N)_$p"
-    end
-    return name
-end
-
-end # of module SLNs

SL_orbit.jl

@@ -1,61 +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-14
-        "--iterations"
-            help = "set maximal number of iterations for the SDP solver (default: 20000)"
-            arg_type = Int
-            default = 60000
-        "--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
-        "-N"
-            help = "Consider elementary matrices EL(N)"
-            arg_type = Int
-            default = 2
-        "-p"
-            help = "Matrices over field of p-elements (p=0 => over ZZ)"
-            arg_type = Int
-            default = 0
-        "--radius"
-            help = "Radius of ball B_r(e,S) to find solution over"
-            arg_type = Int
-            default = 2
-        "-X"
-            help = "Consider EL(N, ZZ⟨X⟩)"
-            action = :store_true
-        "--warmstart"
-            help = "Use warmstart.jl as the initial guess for SCS"
-            action = :store_true
-    end
-
-    return parse_args(settings)
-end
-
-parsed_args = parse_commandline()
-
-include("CPUselect.jl")
-
-set_parallel_mthread(parsed_args, workers=true)
-
-include("SLNs.jl")
-include("Orbit.jl")
-
-main(SLNs, parsed_args)

groups/Allgroups.jl

@@ -0,0 +1,56 @@
+module PropertyTGroups
+
+using PropertyT
+using AbstractAlgebra
+using Nemo
+using Groups
+using GroupRings
+
+export PropertyTGroup, SymmetrizedGroup, GAPGroup,
+    SpecialLinearGroup,
+    SpecialAutomorphismGroup,
+    HigmanGroup,
+    CapraceGroup,
+    MappingClassGroup
+
+export PropertyTGroup
+
+abstract type 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
+
+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

@@ -0,0 +1,21 @@
+struct SpecialAutomorphismGroup{N} <: SymmetrizedGroup
+    group::AutGroup
+end
+
+function SpecialAutomorphismGroup(args::Dict)
+    N = args["SAut"]
+    return SpecialAutomorphismGroup{N}(AutGroup(FreeGroup(N), special=true))
+end
+
+name(G::SpecialAutomorphismGroup{N}) where N = "SAutF$(N)"
+
+group(G::SpecialAutomorphismGroup) = G.group
+
+function generatingset(G::SpecialAutomorphismGroup)
+    S = gens(group(G));
+    return unique([S; inv.(S)])
+end
+
+function autS(G::SpecialAutomorphismGroup{N}) where N
+    return WreathProduct(PermutationGroup(2), PermutationGroup(N))
+end

groups/caprace.jl

@@ -0,0 +1,35 @@
+struct CapraceGroup <: GAPGroup end
+
+name(G::CapraceGroup) = "CapraceGroup"
+
+function group(G::CapraceGroup)
+
+   caprace_group = Groups.FPGroup(["x","y","z","t","r"])
+
+   x,y,z,t,r = gens(caprace_group)
+
+   relations = [
+   x^7,
+   y^7,
+   t^2,
+   r^73,
+   t*r*t*r,
+   Comm(x,y)*z^-1,
+   Comm(x,z),
+   Comm(y,z),
+   Comm(x^2*y*z^-1, t),
+   Comm(x*y*z^3, t*r),
+   Comm(x^3*y*z^2, t*r^17),
+   Comm(x, t*r^-34),
+   Comm(y, t*r^-32),
+   Comm(z, t*r^-29),
+   Comm(x^-2*y*z, t*r^-25),
+   Comm(x^-1*y*z^-3, t*r^-19),
+   Comm(x^-3*y*z^-2, t*r^-11)
+   ];
+
+   relations = [relations; [inv(rel) for rel in relations]]
+
+   Groups.add_rels!(caprace_group, Dict(rel => caprace_group() for rel in relations))
+   return caprace_group
+end

groups/higman.jl

@@ -0,0 +1,22 @@
+struct HigmanGroup <: GAPGroup end
+
+name(G::HigmanGroup) = "HigmanGroup"
+
+function group(G::HigmanGroup)
+
+   higman_group = Groups.FPGroup(["a","b","c","d"]);
+
+   a,b,c,d = gens(higman_group)
+
+   relations = [
+   b*Comm(b,a),
+   c*Comm(c,b),
+   d*Comm(d,c),
+   a*Comm(a,d)
+   ];
+
+   relations = [relations; [inv(rel) for rel in relations]]
+
+   Groups.add_rels!(higman_group, Dict(rel => higman_group() for rel in relations))
+   return higman_group
+end

groups/mappingclassgroup.jl

@@ -0,0 +1,83 @@
+struct MappingClassGroup{N} <: GAPGroup end
+
+MappingClassGroup(args::Dict) = MappingClassGroup{args["MCG"]}()
+
+name(G::MappingClassGroup{N}) where N = "MCG(N)"
+
+function group(G::MappingClassGroup{N}) where 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)
+        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]...,
+           [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
+           ]
+
+       relations = [relations; [inv(rel) for rel in relations]]
+       Groups.add_rels!(MCGroup, Dict(rel => MCGroup() for rel in relations))
+       return MCGroup
+
+    else
+        MCGroup = Groups.FPGroup(["a$i" for i in 0:2N])
+        S = gens(MCGroup)
+
+        a0 = S[1]
+        A = S[2:end]
+        k = length(A)
+
+        relations = [
+            [Comm(A[i], A[j]) for i in 1:k for j in 1:k if abs(i-j) > 1]...,
+            [Comm(a0, A[i]) for i in 1:k if i != 4]...,
+            [A[i]*A[(i+1)]*A[i]*inv(A[i+1]*A[i]*A[i+1]) for i in 1:k-1]...,
+            A[4]*a0*A[4]*inv(a0*A[4]*a0)
+            ]
+
+        # 3-chain relation
+        c = prod(reverse(A[1:4]))*prod(A[1:4])
+        b0 = c*a0*inv(c)
+        push!(relations, (A[1]*A[2]*A[3])^4*inv(a0*b0))
+
+        # Lantern relation
+        b1 = inv(A[4]*A[5]*A[3]*A[4])*a0*(A[4]*A[5]*A[3]*A[4])
+        b2 = inv(A[2]*A[3]*A[1]*A[2])*b1*(A[2]*A[3]*A[1]*A[2])
+        u  = inv(A[6]*A[5])*b1*(A[6]*A[5])
+        x  = prod(reverse(A[2:6]))*u*prod(inv.(A[1:4]))
+        b3 = x*a0*inv(x)
+        push!(relations, a0*b2*b1*inv(A[1]*A[3]*A[5]*b3))
+
+        # Hyperelliptic relation
+        X = prod(reverse(A))*prod(A)
+
+        function n(i::Int, b=b0)
+            if i == 1
+                return A[1]
+            elseif i == 2
+                return b
+            else
+                return w(i-2)*n(i-2)*w(i-2)
+            end
+        end
+
+        function w(i::Int)
+           (A[2i+4]*A[2i+3]*A[2i+2]* n(i+1))*(A[2i+1]*A[2i]  *A[2i+2]*A[2i+1])*
+           (A[2i+3]*A[2i+2]*A[2i+4]*A[2i+3])*( n(i+1)*A[2i+2]*A[2i+1]*A[2i]  )
+        end
+
+        # push!(relations, X*n(N)*inv(n(N)*X))
+
+        relations = [relations; [inv(rel) for rel in relations]]
+        Groups.add_rels!(MCGroup, Dict(rel => MCGroup() for rel in relations))
+
+        return MCGroup
+    end
+end

groups/speciallinear.jl

@@ -0,0 +1,62 @@
+struct SpecialLinearGroup{N} <: SymmetrizedGroup
+    group::AbstractAlgebra.Group
+    p::Int
+    X::Bool
+end
+
+function SpecialLinearGroup(args::Dict)
+    N = args["SL"]
+    p = args["p"]
+    X = args["X"]
+
+    if p == 0
+        G = MatrixSpace(Nemo.ZZ, N, N)
+    else
+        R = Nemo.NmodRing(UInt(p))
+        G = MatrixSpace(R, N, N)
+    end
+    return SpecialLinearGroup{N}(G, p, X)
+end
+
+function name(G::SpecialLinearGroup{N}) where N
+    if G.p == 0
+       R = (G.X ? "Z[x]" : "Z")
+    else
+       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)
+    m[i,j] = val
+    return m
+end
+
+function autS(G::SpecialLinearGroup{N}) where N
+    return WreathProduct(PermutationGroup(2), PermutationGroup(N))
+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

@@ -0,0 +1,61 @@
+using PropertyT
+
+include("FPGroups_GAP.jl")
+
+include("groups/Allgroups.jl")
+using PropertyTGroups
+
+import PropertyT.Settings
+
+function summarize(sett::PropertyT.Settings)
+    info("Threads:     $(Threads.nthreads())")
+    info("Workers:     $(workers())")
+    info("GroupDir:    $(PropertyT.prepath(sett))")
+    info(string(sett.G))
+    info("with generating set of size $(length(sett.S))")
+
+    info("Radius:      $(sett.radius)")
+    info("Precision:   $(sett.tol)")
+    info("Upper bound: $(sett.upper_bound)")
+    info("Solver:      $(sett.solver)")
+end
+
+function Settings(Gr::PropertyTGroup, args, solver)
+    r = get(args, "radius", 2)
+    gr_name = PropertyTGroups.name(Gr)*"_r$r"
+    G = PropertyTGroups.group(Gr)
+    S = PropertyTGroups.generatingset(Gr)
+
+    sol = solver
+    ub = get(args,"upper-bound", Inf)
+    tol = get(args,"tol", 1e-10)
+    ws = get(args, "warmstart", false)
+
+    if get(args, "nosymmetry", false)
+        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
+
+function main(::PropertyTGroup, sett::PropertyT.Settings)
+    isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
+
+    summarize(sett)
+
+    return PropertyT.check_property_T(sett)
+end
+
+function main(::GAPGroup, sett::PropertyT.Settings)
+    isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
+
+    summarize(sett)
+
+    S = [s for s in sett.S if s.symbols[1].pow == 1]
+    relations = [k*inv(v) for (k,v) in sett.G.rels]
+
+    prepare_pm_delta(PropertyT.prepath(sett), GAP_groupcode(S, relations), sett.radius)
+
+    return PropertyT.check_property_T(sett)
+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