diff --git a/FPGroups_GAP.jl b/FPGroups_GAP.jl new file mode 100644 index 0000000..8be68d5 --- /dev/null +++ b/FPGroups_GAP.jl @@ -0,0 +1,155 @@ +using JLD + +function GAP_code(group_code, dir, R; maxeqns=10_000, infolevel=2) + code = """ +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.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, "\\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(`gap -q $(GAP_file)`) +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) + o = GAP_execute(gap_code, name) + return o +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 diff --git a/MCG.jl b/MCG.jl new file mode 100644 index 0000000..c15300e --- /dev/null +++ b/MCG.jl @@ -0,0 +1,115 @@ +using ArgParse +using JLD + +using Nemo +import SCS.SCSSolver +using PropertyT + +using Groups + +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() + args = ArgParseSettings() + + @add_arg_table args 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 mapping class group of surface of genus N" + arg_type = Int + default = 2 + "--radius" + help = "Radius of ball B_r(e,S) to find solution over" + arg_type = Int + default = 4 + "--warmstart" + help = "Use warmstart.jl as the initial guess for SCS" + action = :store_true + end + + return parse_args(args) +end + +include("FPGroups_GAP.jl") + +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 + + tol = parsed_args["tol"] + iterations = parsed_args["iterations"] + upper_bound = parsed_args["upper-bound"] + radius = parsed_args["radius"] + N = parsed_args["N"] + + solver = SCSSolver(eps=tol, max_iters=iterations, linearsolver=SCS.Direct) + + prefix = "MCG($N)" + name = "$(prefix)" + + isdir(name) || mkdir(name) + + prepare_pm_delta(prefix, name, radius) + + logger = PropertyT.setup_logging(name) + + info(logger, "Group: $name") + info(logger, "Iterations: $iterations") + info(logger, "Precision: $tol") + info(logger, "Upper bound: $upper_bound") + + MCGroup = Groups.FPGroup(["a1","a2","a3","a4", "a5"]); + S = Nemo.gens(MCGroup) + Comm(x,y) = x*y*x^-1*y^-1 + k = length(S) + + relations = [[Comm(S[i], S[j]) for i in 1:k for j in 1:k 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:k-1]..., + (S[1]*S[2]*S[3])^4*inv(S[5])^5, + Comm(prod(reverse(S))*prod(S), S[1]), + (prod(reverse(S))*prod(S))^2 + ]; + + relations = [relations; [inv(rel) for rel in relations]] + + Groups.add_rels!(MCGroup, Dict(rel => MCGroup() for rel in relations)) + + S = gens(MCGroup) + S = unique([S; [inv(s) for s in S]]) + Id = MCGroup() + + @time PropertyT.check_property_T(name, S, Id, solver, upper_bound, tol, radius) + return 0 +end + +main()