diff --git a/REQUIRE b/REQUIRE index f730591..9513fd8 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,9 +1,8 @@ julia -JuMP 0.18.0 -SCS -IntervalArithmetic -JLD -Memento AbstractAlgebra Groups GroupRings +IntervalArithmetic +JuMP 0.18.0 +JLD +SCS diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index dfe0e9f..96fba73 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -43,56 +43,52 @@ function augIdproj(Q::AbstractArray{T,2}) where {T<:Real} return R end -function augIdproj(Q::AbstractArray{T,2}, logger) where {T<:Real} - info(logger, "Projecting columns of Q to the augmentation ideal...") - @logtime logger Q = augIdproj(Q) - - info(logger, "Checking that sum of every column contains 0.0... ") - check = all([zero(T) in sum(view(Q, :, i)) for i in 1:size(Q, 2)]) - info(logger, (check? "They do." : "FAILED!")) - - @assert check - - return Q -end - -function distance_to_cone(Δ::GroupRingElem, λ, Q; wlen::Int=4, logger=getlogger()) - info(logger, "------------------------------------------------------------") - info(logger, "Checking in floating-point arithmetic...") - info(logger, "λ = $λ") - @logtime logger sos = compute_SOS(parent(Δ), Q) +function distance_to_cone(Δ::GroupRingElem, λ, Q; wlen::Int=4) + info("------------------------------------------------------------") + info("Checking in floating-point arithmetic...") + info("λ = $λ") + @time sos = compute_SOS(parent(Δ), Q) residue = Δ^2-λ*Δ - sos - info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", aug(residue)))") + info("ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", aug(residue)))") L1_norm = norm(residue,1) - info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", L1_norm))") + info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", L1_norm))") distance = λ - 2^(wlen-1)*L1_norm - info(logger, "Floating point distance (to positive cone) ≈") - info(logger, "$(@sprintf("%.10f", distance))") - info(logger, "") + info("Floating point distance (to positive cone) ≈") + info("$(@sprintf("%.10f", distance))") + info("") if distance ≤ 0 return distance end + info("------------------------------------------------------------") + info("Checking in interval arithmetic...") + info("λ ∈ $λ") + λ = @interval(λ) eoi = Δ^2 - λ*Δ - Q = augIdproj(Q, logger) + info("Projecting columns of Q to the augmentation ideal...") + T = eltype(Q) + @time Q = augIdproj(Q) - info(logger, "------------------------------------------------------------") - info(logger, "Checking in interval arithmetic...") - info(logger, "λ ∈ $λ") - @logtime logger sos = compute_SOS(parent(Δ), Q) + info("Checking that sum of every column contains 0.0... ") + check = all([zero(T) in sum(view(Q, :, i)) for i in 1:size(Q, 2)]) + info((check? "They do." : "FAILED!")) + + @assert check + + @time sos = compute_SOS(parent(Δ), Q) residue = Δ^2-λ*Δ - sos - info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(aug(residue))") + info("ɛ(∑ξᵢ*ξᵢ) ∈ $(aug(residue))") L1_norm = norm(residue,1) - info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)") + info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)") distance = λ - 2^(wlen-1)*L1_norm - info(logger, "The Augmentation-projected distance (to positive cone) ∈") - info(logger, "$(distance)") - info(logger, "") + info("The Augmentation-projected distance (to positive cone) ∈") + info("$(distance)") + info("") return distance.lo end diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 97779fd..bdfc1df 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -14,7 +14,6 @@ immutable Settings{T<:AbstractMathProgSolver} upper_bound::Float64 tol::Float64 warmstart::Bool - logger end prefix(s::Settings) = s.name @@ -167,8 +166,8 @@ function init_model(n, sizes) end function create_SDP_problem(sett::Settings) - info(sett.logger, "Loading orbit data....") - @logtime sett.logger SDP_problem, orb_data = OrbitData(sett); + info("Loading orbit data....") + @time SDP_problem, orb_data = OrbitData(sett); if sett.upper_bound < Inf λ = JuMP.getvariable(SDP_problem, :λ) @@ -176,8 +175,8 @@ function create_SDP_problem(sett::Settings) end t = length(orb_data.laplacian) - info(sett.logger, "Adding $t constraints ... ") - @logtime sett.logger addconstraints!(SDP_problem, orb_data) + info("Adding $t constraints ... ") + @time addconstraints!(SDP_problem, orb_data) return SDP_problem, orb_data end @@ -190,14 +189,14 @@ function λandP(m::JuMP.Model, data::OrbitData, warmstart=true) end function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) - info(sett.logger, "Solving SDP problem...") - @logtime sett.logger λ, Ps = λandP(m, data, sett.warmstart) + info("Solving SDP problem...") + @time λ, Ps = λandP(m, data, sett.warmstart) - info(sett.logger, "Reconstructing P...") + info("Reconstructing P...") preps = load_preps(filename(prepath(sett), :preps), sett.autS) - @logtime sett.logger recP = reconstruct_sol(preps, data.Us, Ps, data.dims) + @time recP = reconstruct_sol(preps, data.Us, Ps, data.dims) fname = filename(fullpath(sett), :P) save(fname, "origP", Ps, "P", recP) @@ -223,7 +222,7 @@ function check_property_T(sett::Settings) files_exists = ex.([:pm, :Δ, :Uπs, :orb, :preps]) if !all(files_exists) - compute_orbit_data(sett.logger, prepath(sett), sett.S, sett.autS, radius=sett.radius) + compute_orbit_data(prepath(sett), sett.S, sett.autS, radius=sett.radius) end cond1 = exists(filename(fullpath(sett), :λ)) @@ -232,21 +231,21 @@ function check_property_T(sett::Settings) if !sett.warmstart && cond1 && cond2 λ, P = λandP(fullpath(sett)) else - info(sett.logger, "Creating SDP problem...") + info("Creating SDP problem...") SDP_problem, orb_data = create_SDP_problem(sett) JuMP.setsolver(SDP_problem, sett.solver) - info(sett.logger, Base.repr(SDP_problem)) + info(Base.repr(SDP_problem)) λ, P = λandP(SDP_problem, orb_data, sett) end - info(sett.logger, "λ = $λ") - info(sett.logger, "sum(P) = $(sum(P))") - info(sett.logger, "maximum(P) = $(maximum(P))") - info(sett.logger, "minimum(P) = $(minimum(P))") + info("λ = $λ") + info("sum(P) = $(sum(P))") + info("maximum(P) = $(maximum(P))") + info("minimum(P) = $(minimum(P))") isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || warn("The solution matrix doesn't seem to be positive definite!") - return interpret_results(sett.name, sett.S, sett.radius, λ, P, sett.logger) + return interpret_results(sett.name, sett.S, sett.radius, λ, P) end diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index b0c5ad4..67ec596 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -128,51 +128,51 @@ function orthSVD{T}(M::AbstractMatrix{T}) return fact[:U][:,1:M_rank] end -function compute_orbit_data{T<:GroupElem}(logger, name::String, S::Vector{T}, autS::Group; radius=2) +function compute_orbit_data{T<:GroupElem}(name::String, S::Vector{T}, autS::Group; radius=2) isdir(name) || mkdir(name) - info(logger, "Generating ball of radius $(2*radius)") + info("Generating ball of radius $(2*radius)") # TODO: Fix that by multiple dispatch? G = parent(first(S)) Id = (isa(G, Ring) ? one(G) : G()) - @logtime logger E_2R, sizes = Groups.generate_balls(S, Id, radius=2*radius); - info(logger, "Balls of sizes $sizes.") - info(logger, "Reverse dict") - @logtime logger E_rdict = GroupRings.reverse_dict(E_2R) + @time E_2R, sizes = Groups.generate_balls(S, Id, radius=2*radius); + info("Balls of sizes $sizes.") + info("Reverse dict") + @time E_rdict = GroupRings.reverse_dict(E_2R) - info(logger, "Product matrix") - @logtime logger pm = GroupRings.create_pm(E_2R, E_rdict, sizes[radius], twisted=true) + info("Product matrix") + @time pm = GroupRings.create_pm(E_2R, E_rdict, sizes[radius], twisted=true) RG = GroupRing(G, E_2R, E_rdict, pm) Δ = PropertyT.spLaplacian(RG, S) @assert GroupRings.aug(Δ) == 0 save(filename(name, :Δ), "Δ", Δ.coeffs) save(filename(name, :pm), "pm", pm) - info(logger, "Decomposing E into orbits of $(autS)") - @logtime logger orbs = orbit_decomposition(autS, E_2R, E_rdict) + info("Decomposing E into orbits of $(autS)") + @time orbs = orbit_decomposition(autS, E_2R, E_rdict) @assert sum(length(o) for o in orbs) == length(E_2R) - info(logger, "E consists of $(length(orbs)) orbits!") + info("E consists of $(length(orbs)) orbits!") save(joinpath(name, "orbits.jld"), "orbits", orbs) - info(logger, "Action matrices") - @logtime logger reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict) + info("Action matrices") + @time reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict) save_preps(filename(name, :preps), reps) reps = matrix_reps(reps) - info(logger, "Projections") - @logtime logger autS_mps = Projections.rankOne_projections(GroupRing(autS)); + info("Projections") + @time autS_mps = Projections.rankOne_projections(GroupRing(autS)); - @logtime logger π_E_projections = [Cstar_repr(p, reps) for p in autS_mps] + @time π_E_projections = [Cstar_repr(p, reps) for p in autS_mps] - info(logger, "Uπs...") - @logtime logger Uπs = orthSVD.(π_E_projections) + info("Uπs...") + @time Uπs = orthSVD.(π_E_projections) multiplicities = size.(Uπs,2) - info(logger, "multiplicities = $multiplicities") + info("multiplicities = $multiplicities") dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps]; - info(logger, "dimensions = $dimensions") + info("dimensions = $dimensions") @assert dot(multiplicities, dimensions) == sizes[radius] save(joinpath(name, "U_pis.jld"), diff --git a/src/PropertyT.jl b/src/PropertyT.jl index d27a9a6..794a299 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -1,3 +1,4 @@ +__precompile__() module PropertyT using AbstractAlgebra @@ -10,55 +11,6 @@ using JLD using JuMP using MathProgBase -using Memento - -function setup_logging(name::String, handlername::Symbol) - isdir(name) || mkdir(name) - logger = Memento.config!("info", fmt="{date}| {msg}") - handler = DefaultHandler(filename(name, Symbol(handlername)), 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 - exists(fname::String) = isfile(fname) || islink(fname) filename(prefix, s::Symbol) = filename(prefix, Val{s}) @@ -92,15 +44,14 @@ function Laplacian(name::String, G::Group) return Δ end -function Laplacian{T<:GroupElem}(S::Vector{T}, Id::T, - logger=getlogger(); radius::Int=2) +function Laplacian{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2) - info(logger, "Generating metric ball of radius $radius...") - @logtime logger E_R, sizes = Groups.generate_balls(S, Id, radius=2*radius) - info(logger, "Generated balls of sizes $sizes.") + info("Generating metric ball of radius $radius...") + @time E_R, sizes = Groups.generate_balls(S, Id, radius=2*radius) + info("Generated balls of sizes $sizes.") - info(logger, "Creating product matrix...") - @logtime logger pm = GroupRings.create_pm(E_R, GroupRings.reverse_dict(E_R), sizes[radius]; twisted=true) + info("Creating product matrix...") + @time pm = GroupRings.create_pm(E_R, GroupRings.reverse_dict(E_R), sizes[radius]; twisted=true) RG = GroupRing(parent(Id), E_R, pm) @@ -155,15 +106,14 @@ Kazhdan(λ::Number,N::Integer) = sqrt(2*λ/N) function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius, warm::Bool=false) isdir(name) || mkdir(name) - logger = Memento.getlogger() if exists(filename(name, :pm)) && exists(filename(name, :Δ)) # cached - info(logger, "Loading precomputed Δ...") + info("Loading precomputed Δ...") Δ = Laplacian(name, parent(S[1])) else # compute - Δ = Laplacian(S, Id, logger, radius=radius) + Δ = Laplacian(S, Id, radius=radius) save(filename(name, :pm), "pm", parent(Δ).pm) save(filename(name, :Δ), "Δ", Δ.coeffs) end @@ -171,48 +121,61 @@ function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius, fullpath = joinpath(name, string(upper_bound)) isdir(fullpath) || mkdir(fullpath) - cond1 = exists(filename(fullpath, :λ)) - cond2 = exists(filename(fullpath, :P)) + files_exist = exists(filename(fullpath, :λ)) && exists(filename(fullpath, :P)) - if !(warm) && cond1 && cond2 - info(logger, "Loading precomputed λ, P...") + if !(warm) && files_exist + info("Loading precomputed λ, P...") λ, P = λandP(fullpath) else - info(logger, "Creating SDP problem...") + info("Creating SDP problem...") SDP_problem, varλ, varP = create_SDP_problem(Δ, constraints(parent(Δ).pm), upper_bound=upper_bound) JuMP.setsolver(SDP_problem, solver) - info(logger, Base.repr(SDP_problem)) + info(Base.repr(SDP_problem)) - @logtime logger λ, P = λandP(fullpath, SDP_problem, varλ, varP) + if warm && isfile(filename(name, :warm)) + ws = load(filename(name, :warm), "warmstart") + else + ws = nothing + end + + @time λ, P, ws = λandP(SDP_problem, varλ, varP, warmstart=ws, solverlog=filename(name, :solverlog)) + + if λ > 0 + save(filename(name, :λ), "λ", λ) + save(filename(name, :P), "P", P) + save(filename(name, :warm), "warmstart", ws) + else + throw(ErrorException("Solver did not produce a valid solution: λ = $λ")) + end end - info(logger, "λ = $λ") - info(logger, "sum(P) = $(sum(P))") - info(logger, "maximum(P) = $(maximum(P))") - info(logger, "minimum(P) = $(minimum(P))") + info("λ = $λ") + info("sum(P) = $(sum(P))") + info("maximum(P) = $(maximum(P))") + info("minimum(P) = $(minimum(P))") isapprox(eigvals(P), abs.(eigvals(P)), atol=tol) || warn("The solution matrix doesn't seem to be positive definite!") - return interpret_results(name, S, radius, λ, P, logger) + return interpret_results(name, S, radius, λ, P) end -function interpret_results(name, S, radius, λ, P, logger) +function interpret_results(name, S, radius, λ, P) RG = GroupRing(parent(first(S)), load(filename(name, :pm), "pm")) Δ = GroupRingElem(load(filename(name, :Δ), "Δ")[:, 1], RG) - @logtime logger Q = real(sqrtm(Symmetric(P))) + @time Q = real(sqrtm(Symmetric(P))) - sgap = distance_to_cone(Δ, λ, Q, wlen=2*radius, logger=logger) + sgap = distance_to_cone(Δ, λ, Q, wlen=2*radius) if sgap > 0 Kazhdan_κ = Kazhdan(sgap, length(S)) if Kazhdan_κ > 0 - info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!") + info("κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!") return true end end - info(logger, "λ($name, S) ≥ $sgap < 0: Tells us nothing about property (T)") + info("λ($name, S) ≥ $sgap < 0: Tells us nothing about property (T)") return false end