Compare commits

..

No commits in common. "master" and "fix/new_type_system" have entirely different histories.

17 changed files with 449 additions and 815 deletions

4
.gitignore vendored
View File

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

69
AutFn.jl Normal file
View File

@ -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

62
FPgroup.jl Normal file
View File

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

View File

@ -1,10 +1,4 @@
# DEPRECATED!
This repository has not been updated for a while!
If You are interested in replicating results for [1712.07167](https://arxiv.org/abs/1712.07167) please check [these instruction](https://kalmar.faculty.wmi.amu.edu.pl/post/1712.07176/)
Also [this notebook](https://nbviewer.jupyter.org/gist/kalmarek/03510181bc1e7c98615e86e1ec580b2a) could be of some help. If everything else fails the [zenodo dataset](https://zenodo.org/record/1133440) should contain the last-resort instructions.
This repository contains some legacy code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
This repository contains code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
# Installing
To run the code You need `julia-v0.5` (should work on `v0.6`, but with warnings).

66
SLn.jl Normal file
View File

@ -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

View File

@ -1,56 +1,26 @@
module PropertyTGroups
using PropertyT
using AbstractAlgebra
using Nemo
using Groups
using GroupRings
export PropertyTGroup, SymmetrizedGroup, GAPGroup,
SpecialLinearGroup,
SpecialAutomorphismGroup,
HigmanGroup,
CapraceGroup,
MappingClassGroup
export PropertyTGroup
export PropertyTGroup, SymmetricGroup, GAPGroup
abstract type PropertyTGroup end
abstract type SymmetrizedGroup <: PropertyTGroup end
abstract type SymmetricGroup <: 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
generatingset(G::GAPGroup) = gens(group(G))
include("mappingclassgroup.jl")
include("higman.jl")
include("caprace.jl")
include("actions.jl")
end # of module PropertyTGroups

View File

@ -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

View File

@ -1,13 +1,22 @@
struct SpecialAutomorphismGroup{N} <: SymmetrizedGroup
struct SpecialAutomorphismGroup <: SymmetricGroup
args::Dict{String,Any}
group::AutGroup
function SpecialAutomorphismGroup(args::Dict)
N = args["N"]
return new(args, AutGroup(FreeGroup(N), special=true))
end
end
function SpecialAutomorphismGroup(args::Dict)
N = args["SAut"]
return SpecialAutomorphismGroup{N}(AutGroup(FreeGroup(N), special=true))
end
function name(G::SpecialAutomorphismGroup)
N = G.args["N"]
name(G::SpecialAutomorphismGroup{N}) where N = "SAutF$(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

View File

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

View File

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

View File

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

View File

@ -1,55 +1,40 @@
struct SpecialLinearGroup{N} <: SymmetrizedGroup
struct SpecialLinearGroup <: SymmetricGroup
args::Dict{String,Any}
group::AbstractAlgebra.Group
p::Int
X::Bool
function SpecialLinearGroup(args::Dict)
n = args["N"]
p = args["p"]
X = args["X"]
if p == 0
G = MatrixSpace(Nemo.ZZ, n, n)
else
G = MatrixSpace(FiniteField(p), n, n)
end
return new(args, G)
end
end
function SpecialLinearGroup(args::Dict)
N = args["SL"]
p = args["p"]
X = args["X"]
function name(G::SpecialLinearGroup)
N = G.args["N"]
p = G.args["p"]
X = G.args["X"]
if p == 0
G = MatrixSpace(Nemo.ZZ, N, N)
R = (X ? "Z[x]" : "Z")
else
R = Nemo.NmodRing(UInt(p))
G = MatrixSpace(R, N, N)
R = "F$p"
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")
if G.args["nosymmetry"]
return "SL($N,$R)"
else
R = "F$(G.p)"
return "oSL($N,$R)"
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)
@ -57,6 +42,52 @@ function E(i::Int, j::Int, M::MatSpace, val=one(M.base_ring))
return m
end
function autS(G::SpecialLinearGroup{N}) where N
function generatingset(G::SpecialLinearGroup)
n = G.args["N"]
p = G.args["p"]
X = G.args["X"]
SL = group(G)
indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
p > 0 && X && throw("SL(n, F_p[x]) not implemented")
if !X
S = [E(idx[1],idx[2],SL) for idx in indexing]
else
r = G.args["radius"]
S = [E(i,j,SL,v) for (i,j) in indexing for v in [1, 100*r]]
end
return unique([S; inv.(S)])
end
function autS(G::SpecialLinearGroup)
N = G.args["N"]
return WreathProduct(PermutationGroup(2), PermutationGroup(N))
end
###############################################################################
#
# Action of WreathProductElems on Nemo.MatElem
#
###############################################################################
function matrix_emb(n::DirectProductGroupElem, p::perm)
Id = parent(n.elts[1])()
elt = diagm([(-1)^(el == Id ? 0 : 1) for el in n.elts])
return elt[:, p.d]
end
function (g::WreathProductElem)(A::MatElem)
g_inv = inv(g)
G = matrix_emb(g.n, g_inv.p)
G_inv = matrix_emb(g_inv.n, g.p)
M = parent(A)
res = M(G_inv)
Nemo.mul!(res, A, res)
return Nemo.mul!(res, M(G), res)
end
function (p::perm)(A::MatElem)
length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
return p*A*inv(p)
end

View File

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

140
main.jl
View File

@ -1,61 +1,127 @@
using AbstractAlgebra
using Nemo
using PropertyT
using Groups
include("FPGroups_GAP.jl")
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)
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)")
function summarize(logger, groupdir, iterations, tol, upper_bound, radius, G, S)
info(logger, "Group: $groupdir")
info(logger, "Iterations: $iterations")
info(logger, "Precision: $tol")
info(logger, "Upper bound: $upper_bound")
info(logger, "Radius: $radius")
info(logger, "Threads: $(Threads.nthreads())")
info(logger, "Workers: $(workers())")
info(logger, string(G))
info(logger, "with generating set of size $(length(S))")
end
function Settings(Gr::PropertyTGroup, args, solver)
r = get(args, "radius", 2)
gr_name = PropertyTGroups.name(Gr)*"_r$r"
function params(Gr::SymmetricGroup)
radius = Gr.args["radius"]
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=10)
main(G::SymmetricGroup) = main(Symmetrize, G)
function main(::Type{Symmetrize}, Gr::SymmetricGroup)
radius, tol, iterations, upper_bound, warm, N = params(Gr)
groupdir = "$(PropertyTGroups.name(Gr))_r$radius"
isdir(groupdir) || mkdir(groupdir)
logger = PropertyT.setup_logging(joinpath(groupdir, "$(upper_bound)"), :fulllog)
G = PropertyTGroups.group(Gr)
S = PropertyTGroups.generatingset(Gr)
sol = solver
ub = get(args,"upper-bound", Inf)
tol = get(args,"tol", 1e-10)
ws = get(args, "warmstart", false)
summarize(logger, groupdir, iterations, tol, upper_bound, radius, G, S)
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
autS = PropertyTGroups.autS(Gr)
info(logger, "Symmetrising with $(autS)")
function main(::PropertyTGroup, sett::PropertyT.Settings)
isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
solver = scs_solver(tol, iterations)
summarize(sett)
# 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(groupdir, N, G, S, autS,
radius, solver, upper_bound, tol, warm, logger)
return PropertyT.check_property_T(sett)
end
function main(::GAPGroup, sett::PropertyT.Settings)
isdir(PropertyT.fullpath(sett)) || mkpath(PropertyT.fullpath(sett))
function main(::Type{Standard}, Gr::SymmetricGroup)
summarize(sett)
radius, tol, iterations, upper_bound, warm, _ = params(Gr)
S = [s for s in sett.S if s.symbols[1].pow == 1]
relations = [k*inv(v) for (k,v) in sett.G.rels]
groupdir = "$(PropertyTGroups.name(Gr))_r$radius"
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)
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)
return PropertyT.check_property_T(sett)
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

View File

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

108
run.jl
View File

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

View File

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