.gitignore

@@ -11,5 +11,3 @@ SL*_*
 *.gws
 .*
 tests*
-*.py
-*.pyc

AutFn.jl

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

@@ -0,0 +1,62 @@
+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")
+include("FPGroups_GAP.jl")
+
+if PARSEDARGS["Caprace"]
+    G = PropertyTGroups.CapraceGroup(PARSEDARGS)
+elseif PARSEDARGS["Higman"]
+    G = PropertyTGroups.HigmanGroup(PARSEDARGS)
+elseif PARSEDARGS["MCG"] != nothing
+    G = PropertyTGroups.MappingClassGroup(PARSEDARGS)
+else
+    throw("You need to specify one of the --Higman, --Caprace, --MCG N")
+end
+
+main(G)

README.md

@@ -1,10 +1,4 @@
-# DEPRECATED!
+This repository contains code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
-
-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

@@ -0,0 +1,66 @@
+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,19 +1,10 @@
 module PropertyTGroups
 
-using PropertyT
 using AbstractAlgebra
 using Nemo
 using Groups
-using GroupRings
 
-export PropertyTGroup, SymmetrizedGroup, GAPGroup,
+export PropertyTGroup, SymmetrizedGroup, GAPGroup
-    SpecialLinearGroup,
-    SpecialAutomorphismGroup,
-    HigmanGroup,
-    CapraceGroup,
-    MappingClassGroup
-
-export PropertyTGroup
 
 abstract type PropertyTGroup end
 
@@ -21,36 +12,15 @@ 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)
+generatingset(G::GAPGroup) = gens(group(G))
-    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

@@ -1,92 +0,0 @@
-function (p::perm)(A::GroupRingElem)
-    RG = parent(A)
-    result = zero(RG, eltype(A.coeffs))
-    
-    for (idx, c) in enumerate(A.coeffs)
-        if c!= zero(eltype(A.coeffs))
-            result[p(RG.basis[idx])] = c
-        end
-    end
-    return result
-end
-
-###############################################################################
-#
-#  Action of WreathProductElems on Nemo.MatElem
-#
-###############################################################################
-
-function matrix_emb(n::DirectProductGroupElem, p::perm)
-    Id = parent(n.elts[1])()
-    elt = diagm([(-1)^(el == Id ? 0 : 1) for el in n.elts])
-    return elt[:, p.d]
-end
-
-function (g::WreathProductElem)(A::MatElem)
-    g_inv = inv(g)
-    G = matrix_emb(g.n, g_inv.p)
-    G_inv = matrix_emb(g_inv.n, g.p)
-    M = parent(A)
-    return M(G)*A*M(G_inv)
-end
-
-import Base.*
-
-doc"""
-    *(x::AbstractAlgebra.MatElem, P::Generic.perm)
-> Apply the pemutation $P$ to the rows of the matrix $x$ and return the result.
-"""
-function *(x::AbstractAlgebra.MatElem, P::Generic.perm)
-   z = similar(x)
-   m = rows(x)
-   n = cols(x)
-   for i = 1:m
-      for j = 1:n
-         z[i, j] = x[i,P[j]]
-      end
-   end
-   return z
-end
-
-function (p::perm)(A::MatElem)
-    length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
-    return p*A*inv(p)
-end
-
-###############################################################################
-#
-#  Action of WreathProductElems on AutGroupElem
-#
-###############################################################################
-
-function AutFG_emb(A::AutGroup, g::WreathProductElem)
-    isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
-    parent(g).P.n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
-    elt = A()
-    Id = parent(g.n.elts[1])()
-    flips = Groups.AutSymbol[Groups.flip_autsymbol(i) for i in 1:length(g.p.d) if g.n.elts[i] != Id]
-    Groups.r_multiply!(elt, flips, reduced=false)
-    Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])
-    return elt
-end
-
-function AutFG_emb(A::AutGroup, p::perm)
-    isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
-    parent(p).n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(p)) into $A")
-    return A(Groups.perm_autsymbol(p))
-end
-
-function (g::WreathProductElem)(a::Groups.Automorphism)
-    A = parent(a)
-    g = AutFG_emb(A,g)
-    res = A()
-    Groups.r_multiply!(res, g.symbols, reduced=false)
-    Groups.r_multiply!(res, a.symbols, reduced=false)
-    Groups.r_multiply!(res, [inv(s) for s in reverse!(g.symbols)])
-    return res
-end
-
-function (p::perm)(a::Groups.Automorphism)
-    g = AutFG_emb(parent(a),p)
-    return g*a*inv(g)
-end

groups/autfreegroup.jl

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

groups/higman.jl

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

groups/mappingclassgroup.jl

@@ -1,23 +1,26 @@
-struct MappingClassGroup{N} <: GAPGroup end
+struct MappingClassGroup <: GAPGroup
+    args::Dict{String,Any}
+end
 
-MappingClassGroup(args::Dict) = MappingClassGroup{args["MCG"]}()
+function name(G::MappingClassGroup)
-
+    N = G.args["MCG"]
-name(G::MappingClassGroup{N}) where N = "MCG(N)"
+    return "MCG($(N))"
-
+end
-function group(G::MappingClassGroup{N}) where N
 
+function group(G::MappingClassGroup)
+    N = G.args["MCG"]
     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:G.n-1]...,
+           [S[i]*S[i+1]*S[i]*inv(S[i+1]*S[i]*S[i+1]) for i in 1:N-1]...,
            (S[1]*S[2]*S[3])^4*inv(S[5])^2,
            Comm(A, S[1]),
            A^2

groups/speciallinear.jl

@@ -1,44 +1,61 @@
-struct SpecialLinearGroup{N} <: SymmetrizedGroup
+struct SpecialLinearGroup <: SymmetrizedGroup
+    args::Dict{String,Any}
     group::AbstractAlgebra.Group
-    p::Int
-    X::Bool
-end
 
     function SpecialLinearGroup(args::Dict)
-    N = args["SL"]
+        n = args["N"]
         p = args["p"]
         X = args["X"]
 
         if p == 0
-        G = MatrixSpace(Nemo.ZZ, N, N)
+            G = MatrixSpace(Nemo.ZZ, n, n)
         else
             R = Nemo.NmodRing(UInt(p))
-        G = MatrixSpace(R, N, N)
+            G = MatrixSpace(R, n, n)
+        end
+        return new(args, G)
     end
-    return SpecialLinearGroup{N}(G, p, X)
 end
 
-function name(G::SpecialLinearGroup{N}) where N
+function name(G::SpecialLinearGroup)
-    if G.p == 0
+    N = G.args["N"]
-       R = (G.X ? "Z[x]" : "Z")
+    p = G.args["p"]
+    X = G.args["X"]
+
+    if p == 0
+       R = (X ? "Z[x]" : "Z")
     else
-       R = "F$(G.p)"
+       R = "F$p"
+    end
+    if G.args["nosymmetry"]
+        return "SL($N,$R)"
+    else
+        return "oSL($N,$R)"
     end
-    return SL($(G.N),$R)
 end
 
 group(G::SpecialLinearGroup) = G.group
 
-function generatingset(G::SpecialLinearGroup{N}) where N
+function E(i::Int, j::Int, M::MatSpace, val=one(M.base_ring))
-    G.p > 0 && G.X && throw("SL(n, F_p[x]) not implemented")
+    @assert i≠j
-    SL = group(G)
+    m = one(M)
-    return generatingset(SL, G.X)
+    m[i,j] = val
+    return m
 end
 
-# r is the injectivity radius of
+function generatingset(G::SpecialLinearGroup)
-# SL(n, Z[X]) -> SL(n, Z) induced by X -> 100
+    n = G.args["N"]
+    p = G.args["p"]
+    X = G.args["X"]
+    p > 0 && X && throw("SL(n, F_p[x]) not implemented")
+    SL = group(G)
+    r = G.args["radius"]
+    return generatingset(SL, r, X)
+end
+
+
+function generatingset(SL::MatSpace, radius::Integer, X::Bool=false)
 
-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]
 
@@ -50,13 +67,50 @@ function generatingset(SL::MatSpace, X::Bool=false, r=5)
     return unique([S; inv.(S)])
 end
 
-function E(i::Int, j::Int, M::MatSpace, val=one(M.base_ring))
+function autS(G::SpecialLinearGroup)
-    @assert i≠j
+    N = G.args["N"]
-    m = one(M)
-    m[i,j] = val
-    return m
-end
-
-function autS(G::SpecialLinearGroup{N}) where N
     return WreathProduct(PermutationGroup(2), PermutationGroup(N))
 end
+
+###############################################################################
+#
+#  Action of WreathProductElems on Nemo.MatElem
+#
+###############################################################################
+
+function matrix_emb(n::DirectProductGroupElem, p::perm)
+    Id = parent(n.elts[1])()
+    elt = diagm([(-1)^(el == Id ? 0 : 1) for el in n.elts])
+    return elt[:, p.d]
+end
+
+function (g::WreathProductElem)(A::MatElem)
+    g_inv = inv(g)
+    G = matrix_emb(g.n, g_inv.p)
+    G_inv = matrix_emb(g_inv.n, g.p)
+    M = parent(A)
+    return M(G)*A*M(G_inv)
+end
+
+import Base.*
+
+doc"""
+    *(x::AbstractAlgebra.MatElem, P::Generic.perm)
+> Apply the pemutation $P$ to the rows of the matrix $x$ and return the result.
+"""
+function *(x::AbstractAlgebra.MatElem, P::Generic.perm)
+   z = similar(x)
+   m = rows(x)
+   n = cols(x)
+   for i = 1:m
+      for j = 1:n
+         z[i, j] = x[i,P[j]]
+      end
+   end
+   return z
+end
+
+function (p::perm)(A::MatElem)
+    length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
+    return p*A*inv(p)
+end

logging.jl

@@ -1,6 +1,6 @@
 using Memento
 
-function setup_logging(filename::String, handlername::Symbol=:log)
+function setup_logging(filename::String, handlername::Symbol)
     isdir(dirname(filename)) || mkdir(dirname(filename))
     logger = Memento.config!("info", fmt="{date}| {msg}")
     handler = DefaultHandler(filename, DefaultFormatter("{date}| {msg}"))

main.jl

@@ -1,61 +1,140 @@
-using PropertyT
+include("logging.jl")
 
-include("FPGroups_GAP.jl")
+using AbstractAlgebra
+using Nemo
+using PropertyT
+using Groups
+
+using SCS.SCSSolver
+# using Mosek
+# using CSDP
+# using SDPA
 
 include("groups/Allgroups.jl")
 using PropertyTGroups
 
-import PropertyT.Settings
+struct Symmetrize end
+struct Standard end
 
-function summarize(sett::PropertyT.Settings)
+function summarize(groupdir, iterations, tol, upper_bound, radius, G, S)
+    info("Group:       $groupdir")
+    info("Iterations:  $iterations")
+    info("Precision:   $tol")
+    info("Upper bound: $upper_bound")
+    info("Radius:      $radius")
     info("Threads:     $(Threads.nthreads())")
     info("Workers:     $(workers())")
-    info("GroupDir:    $(PropertyT.prepath(sett))")
+    info(string(G))
-    info(string(sett.G))
+    info("with generating set of size $(length(S))")
-    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)
+function params(Gr::SymmetrizedGroup)
-    r = get(args, "radius", 2)
+    radius = Gr.args["radius"]
-    gr_name = PropertyTGroups.name(Gr)*"_r$r"
+    tol = Gr.args["tol"]
+    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
+
+scs_solver(tol, iterations) = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct, alpha=1.95, 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)
+
+function main(Gr::PropertyTGroup)
+    r = Gr.args["radius"]
+    ub = Gr.args["upper-bound"]
+    groupdir = "$(PropertyTGroups.name(Gr))_r$r"
+    isdir(groupdir) || mkdir(groupdir)
+    logfile = PropertyT.filename(joinpath(groupdir, string(ub)), :fulllog)
+
+    logger=setup_logging(logfile, :fulllog)
+
+    if Gr.args["nosymmetry"]
+        return main(Naive, Gr, dir=groupdir)
+    else
+        return main(Symmetrize, Gr, dir=groupdir)
+    end
+end
     G = PropertyTGroups.group(Gr)
     S = PropertyTGroups.generatingset(Gr)
 
-    sol = solver
+    summarize(dir, iterations, tol, upper_bound, radius, G, S)
-    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)
+    info("Symmetrising with $(autS)")
-    end
-end
 
-function main(::PropertyTGroup, sett::PropertyT.Settings)
+    solver = scs_solver(tol, iterations)
-    isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
-
-    summarize(sett)
 
+    sett = Settings(dir, N, G, S, autS,
+        radius, solver, upper_bound, tol, warm)
     return PropertyT.check_property_T(sett)
 end
 
-function main(::GAPGroup, sett::PropertyT.Settings)
+function main(::Type{Standard}, Gr::SymmetrizedGroup)
-    isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
 
-    summarize(sett)
+    radius, tol, iterations, upper_bound, warm, _ = params(Gr)
 
-    S = [s for s in sett.S if s.symbols[1].pow == 1]
+    groupdir = "$(PropertyTGroups.name(Gr))_r$radius"
-    relations = [k*inv(v) for (k,v) in sett.G.rels]
+    isdir(groupdir) || mkdir(groupdir)
+    logger = PropertyT.setup_logging(joinpath(groupdir, "$(upper_bound)"), :fulllog)
 
-    prepare_pm_delta(PropertyT.prepath(sett), GAP_groupcode(S, relations), sett.radius)
+    G = PropertyTGroups.group(Gr)
+    S = PropertyTGroups.generatingset(Gr)
 
-    return PropertyT.check_property_T(sett)
+    summarize(dir, iterations, tol, upper_bound, radius, G, S)
+
+    solver = scs_solver(tol, iterations)
+
+    if G isa AbstractAlgebra.Ring
+        Id = one(G)
+    else
+        Id = G()
+    end
+
+    return PropertyT.check_property_T(groupdir, S, Id,
+        solver, upper_bound, tol, radius, warm)
+
+end
+
+
+function main(Gr::GAPGroup)
+
+    radius, tol, iterations, upper_bound, warm = 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)
+
+    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 = scs_solver(tol, iterations)
+
+    return PropertyT.check_property_T(groupdir, S, G(),
+        solver, upper_bound, tol, radius, warm)
 end

positivity/check_positivity.jl

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

run.jl

@@ -1,108 +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
-
-        "--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

@@ -3,156 +3,78 @@ using Base.Test
 include("main.jl")
 
 testdir = "tests_"*string(now())
-mkdir(testdir)
-include("logging.jl")
-logger=setup_logging(joinpath(testdir, "tests.log"))
 info(testdir)
-
+mkdir(testdir)
 cd(testdir)
 
-# groupname = name(G)
+function SL_tests(::Type{T}, args) where {T<:Union{Standard, Symmetrize}}
-# ub = PARSEDARGS["upper-bound"]
-#
-# fullpath = joinpath(groupname, string(ub))
-# isdir(fullpath) || mkpath(fullpath)
 
-separator(n=60) = info("\n"*("\n"*"="^n*"\n"^3)*"\n")
+    G = PropertyTGroups.SpecialLinearGroup(args)
-
-
-function SL_tests(args)
-
-
-    args["SL"] = 2
     args["p"] = 3
-    G = PropertyTGroup(args)
+    @test main(T, G) == true
-    @test main(G) == true
-    separator()
 
-    let args = args
+    println("\n"*"="^30*"\n")
-        args["SL"] = 2
+
+    begin
         args["p"] = 5
-        G = PropertyTGroup(args)
+        G = PropertyTGroups.SpecialLinearGroup(args)
-        @test main(G) == false
+        @test main(T, G) == false
-        separator()
 
         args["warmstart"] = true
-        G = PropertyTGroup(args)
+        G = PropertyTGroups.SpecialLinearGroup(args)
-        @test main(G) == false
+        @test main(T, G) == false
-        separator()
-
-        args["upper-bound"] = 0.1
-        G = PropertyTGroup(args)
-        @test main(G) == true
-        separator()
     end
 
-    args["SL"] = 2
+    println("\n"*"="^30*"\n")
+
+    begin
         args["p"] = 7
-    G = PropertyTGroup(args)
+        G = PropertyTGroups.SpecialLinearGroup(args)
-    @test main(G) == false
+        @test main(T, G) == false
-    separator()
 
-    args["SL"] = 3
+        println("\n"*"="^30*"\n")
-    args["p"] = 7
-    G = PropertyTGroup(args)
-    @test main(G) == true
-    separator()
 
-    # begin
+        args["upper-bound"] = 0.25
-    #     args["iterations"] = 25000
+        G = PropertyTGroups.SpecialLinearGroup(args)
-    #     args["N"] = 3
+        @test main(T, G) == false
-    #     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)
+    println("\n"*"="^30*"\n")
 
-    G = PropertyTGroup(args)
+    args["N"] = 3
-    @test main(G) == false
+    G = PropertyTGroups.SpecialLinearGroup(args)
-    separator()
+    @test main(T, G) == true
 
+    println("\n"*"="^30*"\n")
+
+    begin
+        args["p"] = 0
+        args["iterations"] = 50000
+        args["upper-bound"] = Inf
+
+        G = PropertyTGroups.SpecialLinearGroup(args)
+        @test main(T, G) == false
+
+        args["upper-bound"] = 0.27
         args["warmstart"] = true
-    G = PropertyTGroup(args)
+        G = PropertyTGroups.SpecialLinearGroup(args)
-    @test main(G) == false
+        @test main(T, G) == false
-    separator()
 
-    args["upper-bound"] = 0.1
+        G = PropertyTGroups.SpecialLinearGroup(args)
-    G = PropertyTGroup(args)
+        @test main(T, G) == true
-    @test main(G) == false
+        G = PropertyTGroups.SpecialLinearGroup(args)
-    separator()
+        @test main(T, G) == true
-
-    return 0
     end
 
-@testset "Groups with(out) (T)" begin
+    return main(T, G)
-
-    @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,
+            "N" => 2,
             "p" => 3,
             "X" => false,
             "iterations"=>50000,
@@ -164,16 +86,16 @@ end
             "nosymmetry"=>true,
         )
 
-            SL_tests(args)
+        @time SL_tests(Standard, args)
     end
 
     @testset "Symmetrized" begin
 
         args = Dict(
-                "SL" => 2,
+            "N" => 2,
             "p" => 3,
             "X" => false,
-                "iterations"=>20000,
+            "iterations"=>50000,
             "tol"=>1e-7,
             "upper-bound"=>Inf,
             "cpus"=>2,
@@ -182,41 +104,6 @@ end
             "nosymmetry"=>false,
         )
 
-            SL_tests(args)
+        @time SL_tests(Symmetrize, 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