kalmar
/
PropertyT.jl
mirror of https://github.com/kalmarek/PropertyT.jl.git
1
0
Fork 0
Code Issues Releases Wiki Activity

Merge branch 'enh/rework-logging' of https://git.wmi.amu.edu.pl/kalmar/PropertyT.jl into enh/rework-logging

This commit is contained in:
kalmarek 2018-01-04 21:46:00 +01:00
parent 596b688c1f cb3c6a0ef5
commit 8413254917
6 changed files with 680 additions and 706 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

140
src/CheckSolution.jl
View File

@ -23,106 +23,112 @@ end
function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int) function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int)
# result = zeros(eltype(Q), l) # result = zeros(eltype(Q), l)
# r = similar(result) # r = similar(result)
# for i in 1:size(Q,2) # for i in 1:size(Q,2)
# print(" $i") # print(" $i")
# result += GroupRings.mul!(r, view(Q,:,i), view(Q,:,i), pm) # result += GroupRings.mul!(r, view(Q,:,i), view(Q,:,i), pm)
# end # end
@everywhere groupring_square = PropertyT.groupring_square @everywhere groupring_square = PropertyT.groupring_square
result = @parallel (+) for i in 1:size(Q,2) result = @parallel (+) for i in 1:size(Q,2)
groupring_square(Q[:,i], l, pm) groupring_square(Q[:,i], l, pm)
end end
println("") return result
return result
end end
function compute_SOS(Q::AbstractArray, RG::GroupRing, l::Int) function compute_SOS(Q::AbstractArray, RG::GroupRing, l::Int)
result = compute_SOS(Q, RG.pm, l) result = compute_SOS(Q, RG.pm, l)
return GroupRingElem(result, RG) return GroupRingElem(result, RG)
end end
function distance_to_cone{T<:Interval}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int) function distances_to_cone(elt::GroupRingElem, wlen::Int)
SOS = compute_SOS(Q, parent(elt), length(elt.coeffs)) ɛ_dist = GroupRings.augmentation(elt)
SOS_diff = elt - SOS
ɛ_dist = GroupRings.augmentation(SOS_diff) eoi_SOS_L1_dist = norm(elt,1)
info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
eoi_SOS_L1_dist = norm(SOS_diff,1) dist = 2^(wlen-1)*eoi_SOS_L1_dist
info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)") return dist, ɛ_dist, eoi_SOS_L1_dist
dist = 2^(wlen-1)*eoi_SOS_L1_dist
return dist
end
function distance_to_cone{T}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int)
SOS = compute_SOS(Q, parent(elt), length(elt.coeffs))
SOS_diff = elt - SOS
ɛ_dist = GroupRings.augmentation(SOS_diff)
info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
eoi_SOS_L1_dist = norm(SOS_diff,1)
info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L1_dist))")
dist = 2^(wlen-1)*eoi_SOS_L1_dist
return dist
end end
function augIdproj{T, I<:AbstractInterval}(S::Type{I}, Q::AbstractArray{T,2}) function augIdproj{T, I<:AbstractInterval}(S::Type{I}, Q::AbstractArray{T,2})
l = size(Q, 2) l = size(Q, 2)
R = zeros(S, (l,l)) R = zeros(S, (l,l))
Threads.@threads for j in 1:l Threads.@threads for j in 1:l
col = sum(view(Q, :,j))/l col = sum(view(Q, :,j))/l
for i in 1:l for i in 1:l
R[i,j] = Q[i,j] - col ± eps(0.0) R[i,j] = Q[i,j] - col ± eps(0.0)
end end
end end
return R return R
end end
function augIdproj{T}(Q::AbstractArray{T,2}, logger) function augIdproj{T}(Q::AbstractArray{T,2}, logger)
info(logger, "Projecting columns of Q to the augmentation ideal...") info(logger, "Projecting columns of Q to the augmentation ideal...")
@logtime logger Q = augIdproj(Interval{T}, Q) @logtime logger Q = augIdproj(Interval{T}, Q)
info(logger, "Checking that sum of every column contains 0.0... ") info(logger, "Checking that sum of every column contains 0.0... ")
check = all([0.0 in sum(view(Q, :, i)) for i in 1:size(Q, 2)]) check = all([0.0 in sum(view(Q, :, i)) for i in 1:size(Q, 2)])
info(logger, (check? "They do." : "FAILED!")) info(logger, (check? "They do." : "FAILED!"))
@assert check @assert check
return Q return Q
end end
function distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen::Int) function distance_to_cone(elt::GroupRingElem, λ::T, Q::AbstractArray{T,2}, wlen::Int, logger) where {T<:AbstractFloat}
info(logger, "------------------------------------------------------------") info(logger, "------------------------------------------------------------")
info(logger, "λ = $λ") info(logger, "λ = $λ")
info(logger, "Checking in floating-point arithmetic...") info(logger, "Checking in floating-point arithmetic...")
Δ²_λΔ = EOI(Δ, λ) @logtime logger SOS_diff = elt - compute_SOS(Q, parent(elt), length(elt.coeffs))
@logtime logger fp_distance = λ - distance_to_cone(Δ²_λΔ, Q, wlen) dist, ɛ_dist, eoi_SOS_L1_dist = distances_to_cone(SOS_diff, wlen)
info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))") info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L1_dist))")
fp_distance = λ - dist
info(logger, "Floating point distance (to positive cone) ≈")
info(logger, "$(@sprintf("%.10f", fp_distance))")
info(logger, "")
return fp_distance
end
function distance_to_cone(elt::GroupRingElem, λ::T, Q::AbstractArray{T,2}, wlen::Int, logger) where {T<:AbstractInterval}
info(logger, "------------------------------------------------------------") info(logger, "------------------------------------------------------------")
info(logger, "λ = $λ")
info(logger, "Checking in interval arithmetic...")
@logtime logger SOS_diff = elt - compute_SOS(Q, parent(elt), length(elt.coeffs))
dist, ɛ_dist, eoi_SOS_L1_dist = distances_to_cone(SOS_diff, wlen)
info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)")
int_distance = λ - dist
info(logger, "The Augmentation-projected actual distance (to positive cone) ∈")
info(logger, "$(int_distance)")
info(logger, "")
return int_distance
end
function check_distance_to_cone(Δ::GroupRingElem, λ, Q, wlen::Int, logger)
fp_distance = distance_to_cone(EOI(Δ, λ), λ, Q, wlen, logger)
if fp_distance ≤ 0 if fp_distance ≤ 0
return fp_distance return fp_distance
end end
info(logger, "")
Q = augIdproj(Q, logger)
info(logger, "Checking in interval arithmetic")
λ = @interval(λ) λ = @interval(λ)
Δ = GroupRingElem([@interval(c) for c in Δ.coeffs], parent(Δ)) Δ = GroupRingElem([@interval(c) for c in Δ.coeffs], parent(Δ))
Δ²_λΔ = EOI(Δ, λ) Q = augIdproj(Q, logger)
@logtime logger Interval_dist_to_ΣSq = λ - distance_to_cone(Δ²_λΔ, Q, wlen) int_distance = distance_to_cone(EOI(Δ, λ), λ, Q, wlen, logger)
info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_ΣSq)")
info(logger, "------------------------------------------------------------")
return Interval_dist_to_ΣSq return int_distance.lo
end end

326
src/Orbit-wise.jl
View File

@ -4,16 +4,17 @@ using SCS
export Settings, OrbitData export Settings, OrbitData
immutable Settings{T<:AbstractMathProgSolver} immutable Settings{T<:AbstractMathProgSolver}
name::String name::String
N::Int N::Int
G::Group G::Group
S::Vector S::Vector
autS::Group autS::Group
radius::Int radius::Int
solver::T solver::T
upper_bound::Float64 upper_bound::Float64
tol::Float64 tol::Float64
warmstart::Bool warmstart::Bool
logger
end end
prefix(s::Settings) = s.name prefix(s::Settings) = s.name
@ -22,43 +23,43 @@ prepath(s::Settings) = prefix(s)
fullpath(s::Settings) = joinpath(prefix(s), suffix(s)) fullpath(s::Settings) = joinpath(prefix(s), suffix(s))
immutable OrbitData{T<:AbstractArray{Float64, 2}, LapType <:AbstractVector{Float64}} immutable OrbitData{T<:AbstractArray{Float64, 2}, LapType <:AbstractVector{Float64}}
name::String name::String
Us::Vector{T} Us::Vector{T}
Ps::Vector{Array{JuMP.Variable,2}} Ps::Vector{Array{JuMP.Variable,2}}
cnstr::Vector{SparseMatrixCSC{Float64, Int}} cnstr::Vector{SparseMatrixCSC{Float64, Int}}
laplacian::LapType laplacian::LapType
laplacianSq::LapType laplacianSq::LapType
dims::Vector{Int} dims::Vector{Int}
end end
function OrbitData(sett::Settings) function OrbitData(sett::Settings)
splap = load(joinpath(prepath(sett), "delta.jld"), "Δ"); splap = load(filename(prepath(sett), :Δ), "Δ");
pm = load(joinpath(prepath(sett), "pm.jld"), "pm"); pm = load(filename(prepath(sett), :pm), "pm");
cnstr = PropertyT.constraints(pm); cnstr = PropertyT.constraints(pm);
splap² = similar(splap) splap² = similar(splap)
splap² = GroupRings.mul!(splap², splap, splap, pm); splap² = GroupRings.mul!(splap², splap, splap, pm);
Uπs = load(joinpath(prepath(sett), "U_pis.jld"), "Uπs") Uπs = load(filename(prepath(sett), :Uπs), "Uπs")
nzros = [i for i in 1:length(Uπs) if size(Uπs[i],2) !=0] nzros = [i for i in 1:length(Uπs) if size(Uπs[i],2) !=0]
Uπs = Uπs[nzros] Uπs = Uπs[nzros]
Uπs = sparsify!.(Uπs, sett.tol, check=true, verbose=true) Uπs = sparsify!.(Uπs, sett.tol, check=true, verbose=true)
#dimensions of the corresponding πs: #dimensions of the corresponding πs:
dims = load(joinpath(prepath(sett), "U_pis.jld"), "dims")[nzros] dims = load(filename(prepath(sett), :Uπs), "dims")[nzros]
m, P = init_model(size(Uπs,1), [size(U,2) for U in Uπs]); m, P = init_model(size(Uπs,1), [size(U,2) for U in Uπs]);
orbits = load(joinpath(prepath(sett), "orbits.jld"), "orbits"); orbits = load(filename(prepath(sett), :orb), "orbits");
n = size(Uπs[1],1) n = size(Uπs[1],1)
orb_spcnstrm = [orbit_constraint(cnstr[collect(orb)], n) for orb in orbits] orb_spcnstrm = [orbit_constraint(cnstr[collect(orb)], n) for orb in orbits]
orb_splap = orbit_spvector(splap, orbits) orb_splap = orbit_spvector(splap, orbits)
orb_splap² = orbit_spvector(splap², orbits) orb_splap² = orbit_spvector(splap², orbits)
orbData = OrbitData(fullpath(sett), Uπs, P, orb_spcnstrm, orb_splap, orb_splap², dims); orbData = OrbitData(fullpath(sett), Uπs, P, orb_spcnstrm, orb_splap, orb_splap², dims);
# orbData = OrbitData(name, Uπs, P, orb_spcnstrm, splap, splap², dims); # orbData = OrbitData(name, Uπs, P, orb_spcnstrm, splap, splap², dims);
return m, orbData return m, orbData
end end
include("OrbitDecomposition.jl") include("OrbitDecomposition.jl")
@ -67,152 +68,139 @@ dens(M::SparseMatrixCSC) = length(M.nzval)/length(M)
dens(M::AbstractArray) = length(findn(M)[1])/length(M) dens(M::AbstractArray) = length(findn(M)[1])/length(M)
function sparsify!{Tv,Ti}(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv); verbose=false) function sparsify!{Tv,Ti}(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv); verbose=false)
n = nnz(M) n = nnz(M)
densM = dens(M) densM = dens(M)
for i in eachindex(M.nzval) for i in eachindex(M.nzval)
if abs(M.nzval[i]) < eps if abs(M.nzval[i]) < eps
M.nzval[i] = zero(Tv) M.nzval[i] = zero(Tv)
end end
end end
dropzeros!(M) dropzeros!(M)
m = nnz(M) m = nnz(M)
if verbose if verbose
info(logger, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M), 20)) info("Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M), 20))
end end
return M return M
end end
function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); check=false, verbose=false) function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); check=false, verbose=false)
densM = dens(M) densM = dens(M)
rankM = rank(M) rankM = rank(M)
M[abs.(M) .< eps] .= zero(T) M[abs.(M) .< eps] .= zero(T)
if check && rankM != rank(M) if check && rankM != rank(M)
warn(logger, "Sparsification decreased the rank!") warn("Sparsification decreased the rank!")
end end
if verbose if verbose
info(logger, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M),20)) info("Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M),20))
end end
return sparse(M) return sparse(M)
end end
sparsify{T}(U::AbstractArray{T}, tol=eps(T); check=true, verbose=false) = sparsify!(deepcopy(U), tol, check=check, verbose=verbose) sparsify{T}(U::AbstractArray{T}, tol=eps(T); check=true, verbose=false) = sparsify!(deepcopy(U), tol, check=check, verbose=verbose)
function init_orbit_data(logger, sett::Settings; radius=2)
ex(fname) = isfile(joinpath(prepath(sett), fname))
files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld", "preps.jld"])
if !all(files_exists)
compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.autS, radius=radius)
end
return 0
end
function transform(U::AbstractArray, V::AbstractArray; sparse=true) function transform(U::AbstractArray, V::AbstractArray; sparse=true)
if sparse if sparse
return sparsify!(U'*V*U) return sparsify!(U'*V*U)
else else
return U'*V*U return U'*V*U
end end
end end
A(data::OrbitData, π, t) = data.dims[π].*transform(data.Us[π], data.cnstr[t]) A(data::OrbitData, π, t) = data.dims[π].*transform(data.Us[π], data.cnstr[t])
function constrLHS(m::JuMP.Model, data::OrbitData, t) function constrLHS(m::JuMP.Model, data::OrbitData, t)
l = endof(data.Us) l = endof(data.Us)
lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l)) lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l))
return lhs return lhs
end end
function constrLHS(m::JuMP.Model, cnstr, Us, Ust, dims, vars, eps=100*eps(1.0)) function constrLHS(m::JuMP.Model, cnstr, Us, Ust, dims, vars, eps=100*eps(1.0))
M = [PropertyT.sparsify!(dims[π].*Ust[π]*cnstr*Us[π], eps) for π in 1:endof(Us)] M = [PropertyT.sparsify!(dims[π].*Ust[π]*cnstr*Us[π], eps) for π in 1:endof(Us)]
return @expression(m, sum(vecdot(M[π], vars[π]) for π in 1:endof(Us))) return @expression(m, sum(vecdot(M[π], vars[π]) for π in 1:endof(Us)))
end end
function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.laplacian); var::Symbol=:λ) function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.laplacian); var::Symbol=:λ)
λ = m[var] λ = m[var]
Ust = [U' for U in data.Us] Ust = [U' for U in data.Us]
idx = [π for π in 1:endof(data.Us) if size(data.Us[π],2) != 0] idx = [π for π in 1:endof(data.Us) if size(data.Us[π],2) != 0]
for t in 1:l for t in 1:l
if t % 100 == 0 if t % 100 == 0
print(t, ", ") print(t, ", ")
end end
# lhs = constrLHS(m, data, t) # lhs = constrLHS(m, data, t)
lhs = constrLHS(m, data.cnstr[t], data.Us[idx], Ust[idx], data.dims[idx], data.Ps[idx]) lhs = constrLHS(m, data.cnstr[t], data.Us[idx], Ust[idx], data.dims[idx], data.Ps[idx])
d, d² = data.laplacian[t], data.laplacianSq[t] d, d² = data.laplacian[t], data.laplacianSq[t]
# if lhs == zero(lhs) # if lhs == zero(lhs)
# if d == 0 && d² == 0 # if d == 0 && d² == 0
# info("Detected empty constraint") # info("Detected empty constraint")
# continue # continue
# else # else
# warn("Adding unsatisfiable constraint!") # warn("Adding unsatisfiable constraint!")
# end # end
# end # end
JuMP.@constraint(m, lhs == d² - λ*d) JuMP.@constraint(m, lhs == d² - λ*d)
end end
println("") println("")
end end
function init_model(n, sizes) function init_model(n, sizes)
m = JuMP.Model(); m = JuMP.Model();
P = Vector{Array{JuMP.Variable,2}}(n) P = Vector{Array{JuMP.Variable,2}}(n)
for (k,s) in enumerate(sizes) for (k,s) in enumerate(sizes)
P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) P[k] = JuMP.@variable(m, [i=1:s, j=1:s])
JuMP.@SDconstraint(m, P[k] >= 0.0) JuMP.@SDconstraint(m, P[k] >= 0.0)
end end
JuMP.@variable(m, λ >= 0.0) JuMP.@variable(m, λ >= 0.0)
JuMP.@objective(m, Max, λ) JuMP.@objective(m, Max, λ)
return m, P return m, P
end end
function create_SDP_problem(sett::Settings) function create_SDP_problem(sett::Settings)
info(logger, "Loading orbit data....") info(sett.logger, "Loading orbit data....")
@logtime logger SDP_problem, orb_data = OrbitData(sett); @logtime sett.logger SDP_problem, orb_data = OrbitData(sett);
if sett.upper_bound < Inf if sett.upper_bound < Inf
λ = JuMP.getvariable(SDP_problem, :λ) λ = JuMP.getvariable(SDP_problem, :λ)
JuMP.@constraint(SDP_problem, λ <= sett.upper_bound) JuMP.@constraint(SDP_problem, λ <= sett.upper_bound)
end end
t = length(orb_data.laplacian) t = length(orb_data.laplacian)
info(logger, "Adding $t constraints ... ") info(sett.logger, "Adding $t constraints ... ")
@logtime logger addconstraints!(SDP_problem, orb_data) @logtime sett.logger addconstraints!(SDP_problem, orb_data)
return SDP_problem, orb_data return SDP_problem, orb_data
end end
function λandP(m::JuMP.Model, data::OrbitData, warmstart=true) function λandP(m::JuMP.Model, data::OrbitData, warmstart=true)
varλ = m[:λ] varλ = m[:λ]
varP = data.Ps varP = data.Ps
λ, Ps = PropertyT.λandP(data.name, m, varλ, varP, warmstart) λ, Ps = PropertyT.λandP(data.name, m, varλ, varP, warmstart)
return λ, Ps return λ, Ps
end end
function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) function λandP(m::JuMP.Model, data::OrbitData, sett::Settings)
info(logger, "Solving SDP problem...") info(sett.logger, "Solving SDP problem...")
λ, Ps = λandP(m, data, sett.warmstart) @logtime sett.logger λ, Ps = λandP(m, data, sett.warmstart)
info(logger, "Reconstructing P...") info(sett.logger, "Reconstructing P...")
preps = load_preps(joinpath(prepath(sett), "preps.jld"), sett.autS) preps = load_preps(filename(prepath(sett), :preps), sett.autS)
@logtime logger recP = reconstruct_sol(preps, data.Us, Ps, data.dims) @logtime sett.logger recP = reconstruct_sol(preps, data.Us, Ps, data.dims)
fname = PropertyT.λSDPfilenames(fullpath(sett))[2] fname = filename(fullpath(sett), :P)
save(fname, "origP", Ps, "P", recP) save(fname, "origP", Ps, "P", recP)
return λ, recP return λ, recP
end end
function load_preps(fname::String, G::Nemo.Group) function load_preps(fname::String, G::Nemo.Group)
@ -229,49 +217,39 @@ end
function check_property_T(sett::Settings) function check_property_T(sett::Settings)
init_orbit_data(logger, sett, radius=sett.radius) ex(s) = exists(filename(prepath(sett), s))
if !sett.warmstart && all(isfile.(λSDPfilenames(fullpath(sett)))) files_exists = ex.([:pm, :Δ, :Uπs, :orb, :preps])
λ, P = PropertyT.λandP(fullpath(sett))
else
info(logger, "Creating SDP problem...")
SDP_problem, orb_data = create_SDP_problem(sett)
JuMP.setsolver(SDP_problem, sett.solver)
λ, P = λandP(SDP_problem, orb_data, sett) if !all(files_exists)
end compute_orbit_data(sett.logger, prepath(sett), sett.S, sett.autS, radius=sett.radius)
end
info(logger, "λ = $λ") cond1 = exists(filename(fullpath(sett), :λ))
info(logger, "sum(P) = $(sum(P))") cond2 = exists(filename(fullpath(sett), :P))
info(logger, "maximum(P) = $(maximum(P))")
info(logger, "minimum(P) = $(minimum(P))")
if λ > 0 if !sett.warmstart && cond1 && cond2
pm_fname, Δ_fname = pmΔfilenames(prepath(sett)) λ, P = PropertyT.λandP(fullpath(sett))
RG = GroupRing(sett.G, load(pm_fname, "pm")) else
Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) info(sett.logger, "Creating SDP problem...")
SDP_problem, orb_data = create_SDP_problem(sett)
JuMP.setsolver(SDP_problem, sett.solver)
info(sett.logger, Base.repr(SDP_problem))
isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || λ, P = λandP(SDP_problem, orb_data, sett)
warn("The solution matrix doesn't seem to be positive definite!") end
# @assert P == Symmetric(P)
@logtime logger Q = real(sqrtm(Symmetric(P)))
sgap = distance_to_positive_cone(Δ, λ, Q, 2*sett.radius) info(sett.logger, "λ = $λ")
if isa(sgap, Interval) info(sett.logger, "sum(P) = $(sum(P))")
sgap = sgap.lo info(sett.logger, "maximum(P) = $(maximum(P))")
end info(sett.logger, "minimum(P) = $(minimum(P))")
if sgap > 0
info(logger, "λ ≥ $(Float64(trunc(sgap,12)))") isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) ||
Kazhdan_κ = PropertyT.Kazhdan_from_sgap(sgap, length(sett.S)) warn("The solution matrix doesn't seem to be positive definite!")
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
info(logger, "κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!") if λ > 0
return true return check_λ(sett.name, sett.S, λ, P, sett.radius, sett.logger)
else end
sgap = Float64(trunc(sgap, 12)) info(sett.logger, "κ($(sett.name), S) ≥ $λ < 0: Tells us nothing about property (T)")
info(logger, "λ($(sett.name), S) ≥ $sgap: Group may NOT HAVE property (T)!") return false
return false
end
end
info(logger, "κ($(sett.name), S) ≥ $λ < 0: Tells us nothing about property (T)")
return false
end end

147
src/OrbitDecomposition.jl
View File

@ -57,7 +57,7 @@ function orbit_decomposition(G::Nemo.Group, E::Vector, rdict=GroupRings.reverse_
orbit = zeros(Int, length(elts)) orbit = zeros(Int, length(elts))
a = E[i] a = E[i]
Threads.@threads for i in 1:length(elts) Threads.@threads for i in 1:length(elts)
orbit[i] = rdict[elts[i](a)] orbit[i] = rdict[elts[i](a)]
end end
tovisit[orbit] = false tovisit[orbit] = false
push!(orbits, unique(orbit)) push!(orbits, unique(orbit))
@ -110,110 +110,111 @@ function matrix_reps{T<:GroupElem}(preps::Dict{T,perm})
end end
function perm_repr(g::GroupElem, E::Vector, E_dict) function perm_repr(g::GroupElem, E::Vector, E_dict)
p = Vector{Int}(length(E)) p = Vector{Int}(length(E))
for (i,elt) in enumerate(E) for (i,elt) in enumerate(E)
p[i] = E_dict[g(elt)] p[i] = E_dict[g(elt)]
end end
return p return p
end end
function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E)) function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
elts = collect(elements(G)) elts = collect(elements(G))
l = length(elts) l = length(elts)
preps = Vector{Generic.perm}(l) preps = Vector{Generic.perm}(l)
permG = Nemo.PermutationGroup(length(E)) permG = Nemo.PermutationGroup(length(E))
Threads.@threads for i in 1:l Threads.@threads for i in 1:l
preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict)) preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict))
end end
return Dict(elts[i]=>preps[i] for i in 1:l) return Dict(elts[i]=>preps[i] for i in 1:l)
end end
function perm_reps(S::Vector, autS::Group, radius::Int) function perm_reps(S::Vector, autS::Group, radius::Int)
E, _ = Groups.generate_balls(S, radius=radius) E, _ = Groups.generate_balls(S, radius=radius)
return perm_reps(autS, E) return perm_reps(autS, E)
end end
function reconstruct_sol{T<:GroupElem, S<:perm}(preps::Dict{T, S}, function reconstruct_sol{T<:GroupElem, S<:perm}(preps::Dict{T, S},
Us::Vector, Ps::Vector, dims::Vector) Us::Vector, Ps::Vector, dims::Vector)
l = length(Us) l = length(Us)
transfP = [dims[π].*Us[π]*Ps[π]*Us[π]' for π in 1:l] transfP = [dims[π].*Us[π]*Ps[π]*Us[π]' for π in 1:l]
tmp = [zeros(Float64, size(first(transfP))) for _ in 1:l] tmp = [zeros(Float64, size(first(transfP))) for _ in 1:l]
perms = collect(keys(preps)) perms = collect(keys(preps))
@inbounds Threads.@threads for π in 1:l @inbounds Threads.@threads for π in 1:l
for p in perms for p in perms
BLAS.axpy!(1.0, view(transfP[π], preps[p].d, preps[p].d), tmp[π]) BLAS.axpy!(1.0, view(transfP[π], preps[p].d, preps[p].d), tmp[π])
end end
end end
recP = 1/length(perms) .* sum(tmp) recP = 1/length(perms) .* sum(tmp)
recP[abs.(recP) .< eps(eltype(recP))] = zero(eltype(recP)) recP[abs.(recP) .< eps(eltype(recP))] = zero(eltype(recP))
return recP return recP
end end
function Cstar_repr(x::GroupRingElem{T}, mreps::Dict) where {T} function Cstar_repr(x::GroupRingElem{T}, mreps::Dict) where {T}
return sum(x[i].*mreps[parent(x).basis[i]] for i in findn(x.coeffs)) return sum(x[i].*mreps[parent(x).basis[i]] for i in findn(x.coeffs))
end end
function orthSVD{T}(M::AbstractMatrix{T}) function orthSVD{T}(M::AbstractMatrix{T})
M = full(M) M = full(M)
fact = svdfact(M) fact = svdfact(M)
M_rank = sum(fact[:S] .> maximum(size(M))*eps(T)) M_rank = sum(fact[:S] .> maximum(size(M))*eps(T))
return fact[:U][:,1:M_rank] return fact[:U][:,1:M_rank]
end end
function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S::Vector{T}, autS::Nemo.Group; radius=2) function compute_orbit_data{T<:GroupElem}(logger, name::String, S::Vector{T}, autS::Nemo.Group; radius=2)
isdir(name) || mkdir(name) isdir(name) || mkdir(name)
info(logger, "Generating ball of radius $(2*radius)") info(logger, "Generating ball of radius $(2*radius)")
# TODO: Fix that by multiple dispatch? # TODO: Fix that by multiple dispatch?
Id = (isa(G, Nemo.Ring) ? one(G) : G()) G = parent(first(S))
Id = (isa(G, Nemo.Ring) ? one(G) : G())
@logtime logger E_2R, sizes = Groups.generate_balls(S, Id, radius=2*radius); @logtime logger E_2R, sizes = Groups.generate_balls(S, Id, radius=2*radius);
info(logger, "Balls of sizes $sizes.") info(logger, "Balls of sizes $sizes.")
info(logger, "Reverse dict") info(logger, "Reverse dict")
@logtime logger E_rdict = GroupRings.reverse_dict(E_2R) @logtime logger E_rdict = GroupRings.reverse_dict(E_2R)
info(logger, "Product matrix") info(logger, "Product matrix")
@logtime logger pm = GroupRings.create_pm(E_2R, E_rdict, sizes[radius], twisted=true) @logtime logger pm = GroupRings.create_pm(E_2R, E_rdict, sizes[radius], twisted=true)
RG = GroupRing(G, E_2R, E_rdict, pm) RG = GroupRing(G, E_2R, E_rdict, pm)
Δ = PropertyT.splaplacian(RG, S) Δ = PropertyT.spLaplacian(RG, S)
@assert GroupRings.augmentation(Δ) == 0 @assert GroupRings.augmentation(Δ) == 0
save(joinpath(name, "delta.jld"), "Δ", Δ.coeffs) save(filename(name, :Δ), "Δ", Δ.coeffs)
save(joinpath(name, "pm.jld"), "pm", pm) save(filename(name, :pm), "pm", pm)
info(logger, "Decomposing E into orbits of $(autS)") info(logger, "Decomposing E into orbits of $(autS)")
@logtime logger orbs = orbit_decomposition(autS, E_2R, E_rdict) @logtime logger orbs = orbit_decomposition(autS, E_2R, E_rdict)
@assert sum(length(o) for o in orbs) == length(E_2R) @assert sum(length(o) for o in orbs) == length(E_2R)
info(logger, "E consists of $(length(orbs)) orbits!") info(logger, "E consists of $(length(orbs)) orbits!")
save(joinpath(name, "orbits.jld"), "orbits", orbs) save(joinpath(name, "orbits.jld"), "orbits", orbs)
info(logger, "Action matrices") info(logger, "Action matrices")
@logtime logger reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict) @logtime logger reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict)
save_preps(joinpath(name, "preps.jld"), reps) save_preps(filename(name, :preps), reps)
reps = matrix_reps(reps) reps = matrix_reps(reps)
info(logger, "Projections") info(logger, "Projections")
@logtime logger autS_mps = rankOne_projections(autS); @logtime logger autS_mps = rankOne_projections(autS);
@logtime logger π_E_projections = [Cstar_repr(p, reps) for p in autS_mps] @logtime logger π_E_projections = [Cstar_repr(p, reps) for p in autS_mps]
info(logger, "Uπs...") info(logger, "Uπs...")
@logtime logger Uπs = orthSVD.(π_E_projections) @logtime logger Uπs = orthSVD.(π_E_projections)
multiplicities = size.(Uπs,2) multiplicities = size.(Uπs,2)
info(logger, "multiplicities = $multiplicities") info(logger, "multiplicities = $multiplicities")
dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps]; dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps];
info(logger, "dimensions = $dimensions") info(logger, "dimensions = $dimensions")
@assert dot(multiplicities, dimensions) == sizes[radius] @assert dot(multiplicities, dimensions) == sizes[radius]
save(joinpath(name, "U_pis.jld"), save(joinpath(name, "U_pis.jld"),
"Uπs", Uπs, "Uπs", Uπs,
"dims", dimensions) "dims", dimensions)
return 0 return 0
end end

256
src/Projections.jl
View File

@ -7,17 +7,17 @@
abstract type AbstractCharacter end abstract type AbstractCharacter end
struct PermCharacter <: AbstractCharacter struct PermCharacter <: AbstractCharacter
p::Generic.Partition p::Generic.Partition
end end
struct DirectProdCharacter <: AbstractCharacter struct DirectProdCharacter <: AbstractCharacter
i::Int i::Int
end end
function (chi::PermCharacter)(g::Generic.perm) function (chi::PermCharacter)(g::Generic.perm)
R = Nemo.partitionseq(chi.p) R = Nemo.partitionseq(chi.p)
p = Partition(Nemo.Generic.permtype(g)) p = Partition(Nemo.Generic.permtype(g))
return Int(Nemo.Generic.MN1inner(R, p, 1, Nemo.Generic._charvalsTable)) return Int(Nemo.Generic.MN1inner(R, p, 1, Nemo.Generic._charvalsTable))
end end
Nemo.isone(p::GroupElem) = p == parent(p)() Nemo.isone(p::GroupElem) = p == parent(p)()
@ -29,23 +29,23 @@ end
## NOTE: this works only for Z/2!!!! ## NOTE: this works only for Z/2!!!!
function (chi::DirectProdCharacter)(g::DirectProductGroupElem) function (chi::DirectProdCharacter)(g::DirectProductGroupElem)
return reduce(*, 1, ((-1)^isone(g.elts[j]) for j in 1:chi.i)) return reduce(*, 1, ((-1)^isone(g.elts[j]) for j in 1:chi.i))
end end
for T in [PermCharacter, DirectProdCharacter] for T in [PermCharacter, DirectProdCharacter]
@eval begin @eval begin
function (chi::$T)(X::GroupRingElem) function (chi::$T)(X::GroupRingElem)
RG = parent(X) RG = parent(X)
z = zero(eltype(X)) z = zero(eltype(X))
result = z result = z
for i in 1:length(X.coeffs) for i in 1:length(X.coeffs)
if X.coeffs[i] != z if X.coeffs[i] != z
result += chi(RG.basis[i])*X.coeffs[i] result += chi(RG.basis[i])*X.coeffs[i]
end
end end
end return result
return result end
end end
end
end end
############################################################################### ###############################################################################
@ -55,16 +55,16 @@ end
############################################################################### ###############################################################################
function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int}) function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int})
result = RG(T) result = RG(T)
result.coeffs = full(result.coeffs) result.coeffs = full(result.coeffs)
dim = chi(RG.group()) dim = chi(RG.group())
ord = Int(order(RG.group)) ord = Int(order(RG.group))
for g in RG.basis for g in RG.basis
result[g] = convert(T, (dim//ord)*chi(g)) result[g] = convert(T, (dim//ord)*chi(g))
end end
return result return result
end end
function idempotents(RG::GroupRing{Generic.PermGroup}, T::Type=Rational{Int}) function idempotents(RG::GroupRing{Generic.PermGroup}, T::Type=Rational{Int})
@ -97,132 +97,132 @@ function idempotents(RG::GroupRing{Generic.PermGroup}, T::Type=Rational{Int})
end end
function rankOne_projection(chi::PropertyT.PermCharacter, function rankOne_projection(chi::PropertyT.PermCharacter,
idems::Vector{T}) where {T<:GroupRingElem} idems::Vector{T}) where {T<:GroupRingElem}
RG = parent(first(idems)) RG = parent(first(idems))
S = eltype(first(idems)) S = eltype(first(idems))
ids = [one(RG, S); idems] ids = [one(RG, S); idems]
zzz = zero(S) zzz = zero(S)
for (i,j,k) in Base.product(ids, ids, ids) for (i,j,k) in Base.product(ids, ids, ids)
if chi(i) == zzz || chi(j) == zzz || chi(k) == zzz if chi(i) == zzz || chi(j) == zzz || chi(k) == zzz
continue continue
end end
elt = i*j*k elt = i*j*k
elt^2 == elt || continue elt^2 == elt || continue
if chi(elt) == one(S) if chi(elt) == one(S)
return elt return elt
# return (i,j,k) # return (i,j,k)
end end
end end
throw("Couldn't find rank-one projection for $chi") throw("Couldn't find rank-one projection for $chi")
end end
function rankOne_projections(G::Generic.PermGroup, T::Type=Rational{Int}) function rankOne_projections(G::Generic.PermGroup, T::Type=Rational{Int})
if G.n == 1 if G.n == 1
return [one(GroupRing(G), T)] return [one(GroupRing(G), T)]
elseif G.n < 8 elseif G.n < 8
RG = GroupRing(G, fastm=true) RG = GroupRing(G, fastm=true)
else else
RG = GroupRing(G, fastm=false) RG = GroupRing(G, fastm=false)
end end
RGidems = idempotents(RG, T) RGidems = idempotents(RG, T)
l = length(Partitions(G.n)) l = length(Partitions(G.n))
parts = collect(Partitions(G.n)) parts = collect(Partitions(G.n))
chars = [PropertyT.PermCharacter(p) for p in parts] chars = [PropertyT.PermCharacter(p) for p in parts]
min_projs = Vector{eltype(RGidems)}(l) min_projs = Vector{eltype(RGidems)}(l)
for i in 1:l for i in 1:l
chi = PropertyT.PermCharacter(parts[i]) chi = PropertyT.PermCharacter(parts[i])
min_projs[i] = rankOne_projection(chi,RGidems)*central_projection(RG,chi) min_projs[i] = rankOne_projection(chi,RGidems)*central_projection(RG,chi)
end end
return min_projs return min_projs
end end
function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int}) function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int})
N = BN.P.n N = BN.P.n
# projections as elements of the group rings RSₙ # projections as elements of the group rings RSₙ
SNprojs_nc = [rankOne_projections(PermutationGroup(i)) for i in 1:N] SNprojs_nc = [rankOne_projections(PermutationGroup(i)) for i in 1:N]
# embedding into group ring of BN # embedding into group ring of BN
RBN = GroupRing(BN) RBN = GroupRing(BN)
RFFFF_projs = [central_projection(GroupRing(BN.N), DirectProdCharacter(i),T) RFFFF_projs = [central_projection(GroupRing(BN.N), DirectProdCharacter(i),T)
for i in 1:BN.P.n] for i in 1:BN.P.n]
e0 = central_projection(GroupRing(BN.N), DirectProdCharacter(0), T) e0 = central_projection(GroupRing(BN.N), DirectProdCharacter(0), T)
Q0 = RBN(e0, g -> BN(g)) Q0 = RBN(e0, g -> BN(g))
Qs = [RBN(q, g -> BN(g)) for q in RFFFF_projs] Qs = [RBN(q, g -> BN(g)) for q in RFFFF_projs]
all_projs = [Q0*RBN(p, g->BN(g)) for p in SNprojs_nc[N]] all_projs = [Q0*RBN(p, g->BN(g)) for p in SNprojs_nc[N]]
range = collect(1:N) range = collect(1:N)
for i in 1:N-1 for i in 1:N-1
first_emb = g->BN(Nemo.Generic.emb!(BN.P(), g, range[1:i])) first_emb = g->BN(Nemo.Generic.emb!(BN.P(), g, range[1:i]))
last_emb = g->BN(Nemo.Generic.emb!(BN.P(), g, range[i+1:end])) last_emb = g->BN(Nemo.Generic.emb!(BN.P(), g, range[i+1:end]))
Sk_first = [RBN(p, first_emb) for p in SNprojs_nc[i]] Sk_first = [RBN(p, first_emb) for p in SNprojs_nc[i]]
Sk_last = [RBN(p, last_emb) for p in SNprojs_nc[N-i]] Sk_last = [RBN(p, last_emb) for p in SNprojs_nc[N-i]]
append!(all_projs, append!(all_projs,
[Qs[i]*p1*p2 for (p1,p2) in Base.product(Sk_first,Sk_last)]) [Qs[i]*p1*p2 for (p1,p2) in Base.product(Sk_first,Sk_last)])
end end
append!(all_projs, [Qs[N]*RBN(p, g->BN(g)) for p in SNprojs_nc[N]]) append!(all_projs, [Qs[N]*RBN(p, g->BN(g)) for p in SNprojs_nc[N]])
return all_projs return all_projs
end end
############################################################################## ##############################################################################
# #
# General Groups Misc # General Groups Misc
# #
############################################################################## ##############################################################################
doc""" doc"""
products(X::Vector{GroupElem}, Y::Vector{GroupElem}, op=*) products(X::Vector{GroupElem}, Y::Vector{GroupElem}, op=*)
> Returns a vector of all possible products (or `op(x,y)`), where $x\in X$ and > Returns a vector of all possible products (or `op(x,y)`), where $x\in X$ and
> $y\in Y$ are group elements. You may specify which operation is used when > $y\in Y$ are group elements. You may specify which operation is used when
> forming 'products' by adding `op` (which is `*` by default). > forming 'products' by adding `op` (which is `*` by default).
""" """
function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=*)
result = Vector{T}() result = Vector{T}()
seen = Set{T}() seen = Set{T}()
for x in X for x in X
for y in Y for y in Y
z = op(x,y) z = op(x,y)
if !in(z, seen) if !in(z, seen)
push!(seen, z) push!(seen, z)
push!(result, z) push!(result, z)
end end
end end
end end
return result return result
end end
doc""" doc"""
generateGroup(gens::Vector{GroupElem}, r=2, Id=parent(first(gens))(), op=*) generateGroup(gens::Vector{GroupElem}, r=2, Id=parent(first(gens))(), op=*)
> Produces all elements of a group generated by elements in `gens` in ball of > Produces all elements of a group generated by elements in `gens` in ball of
> radius `r` (word-length metric induced by `gens`). > radius `r` (word-length metric induced by `gens`).
> If `r(=2)` is specified the procedure will terminate after generating ball > If `r(=2)` is specified the procedure will terminate after generating ball
> of radius `r` in the word-length metric induced by `gens`. > of radius `r` in the word-length metric induced by `gens`.
> The identity element `Id` and binary operation function `op` can be supplied > The identity element `Id` and binary operation function `op` can be supplied
> to e.g. take advantage of additive group structure. > to e.g. take advantage of additive group structure.
""" """
function generateGroup{T<:GroupElem}(gens::Vector{T}, r=2, Id::T=parent(first(gens))(), op=*) function generateGroup{T<:GroupElem}(gens::Vector{T}, r=2, Id::T=parent(first(gens))(), op=*)
n = 0 n = 0
R = 1 R = 1
elts = gens elts = gens
gens = [Id; gens] gens = [Id; gens]
while n ≠ length(elts) && R < r while n ≠ length(elts) && R < r
# @show elts # @show elts
R += 1 R += 1
n = length(elts) n = length(elts)
elts = products(elts, gens, op) elts = products(elts, gens, op)
end end
return elts return elts
end end

386
src/PropertyT.jl
View File

@ -12,19 +12,29 @@ using MathProgBase
using Memento using Memento
const logger = Memento.config("info", fmt="{msg}")
const solver_logger = Memento.config("info", fmt="{msg}")
function setup_logging(name::String) function setup_logging(name::String)
isdir(name) || mkdir(name) isdir(name) || mkdir(name)
L = Memento.config("info", fmt="{date}| {msg}")
Memento.add_handler(logger, handler = Memento.DefaultHandler(
Memento.DefaultHandler(joinpath(name,"full_$(string((now()))).log"), filename(name, :logall), Memento.DefaultFormatter("{date}| {msg}"))
Memento.DefaultFormatter("{date}| {msg}")), "full_log")
e = redirect_stderr(logger.handlers["full_log"].io) handler.levels.x = L.levels
L.handlers["all"] = handler
return logger # e = redirect_stderr(L.handlers["all"].io)
return L
end
function solverlogger(name)
logger = Memento.config("info", fmt="{msg}")
handler = DefaultHandler(
filename(name, :logsolver), DefaultFormatter("{date}| {msg}"))
handler.levels.x = logger.levels
logger.handlers["solver_log"] = handler
return logger
end end
macro logtime(logger, ex) macro logtime(logger, ex)
@ -35,8 +45,8 @@ macro logtime(logger, ex)
elapsedtime = Base.time_ns() - elapsedtime elapsedtime = Base.time_ns() - elapsedtime
local diff = Base.GC_Diff(Base.gc_num(), stats) local diff = Base.GC_Diff(Base.gc_num(), stats)
local ts = time_string(elapsedtime, diff.allocd, diff.total_time, local ts = time_string(elapsedtime, diff.allocd, diff.total_time,
Base.gc_alloc_count(diff)) Base.gc_alloc_count(diff))
esc(info(logger, ts)) $(esc(info))($(esc(logger)), ts)
val val
end end
end end
@ -66,289 +76,167 @@ function time_string(elapsedtime, bytes, gctime, allocs)
return str return str
end end
function exists(fname::String) exists(fname::String) = isfile(fname) || islink(fname)
return isfile(fname) || islink(fname)
filename(prefix, s::Symbol) = filename(prefix, Val{s})
@eval begin
for (s,n) in [
[:pm, "pm.jld"],
[:Δ, "delta.jld"],
[:λ, "lambda.jld"],
[:P, "SDPmatrix.jld"],
[:warm, "warmstart.jld"],
[:Uπs, "U_pis.jld"],
[:orb, "orbits.jld"],
[:preps,"preps.jld"],
[:logall, "full_$(string(now())).log"],
[:logsolver,"solver_$(string(now())).log"]
]
filename(prefix::String, ::Type{Val{$:(s)}}) = joinpath(prefix, :($n))
end
end end
function pmΔfilenames(prefix::String) function Laplacian(name::String, G::Group)
isdir(prefix) || mkdir(prefix) if exists(filename(name, :Δ)) && exists(filename(name, :pm))
pm_filename = joinpath(prefix, "pm.jld") RG = GroupRing(G, load(filename(name, :pm), "pm"))
Δ_coeff_filename = joinpath(prefix, "delta.jld") Δ = GroupRingElem(load(filename(name, :Δ), "Δ")[:, 1], RG)
return pm_filename, Δ_coeff_filename else
throw("You need to precompute $(filename(name, :pm)) and $(filename(name, :Δ)) to load it!")
end
return Δ
end end
function λSDPfilenames(prefix::String) function Laplacian{T<:GroupElem}(S::Vector{T}, Id::T,
isdir(prefix) || mkdir(prefix) logger=getlogger(); radius::Int=2)
λ_filename = joinpath(prefix, "lambda.jld")
SDP_filename = joinpath(prefix, "SDPmatrix.jld")
return λ_filename, SDP_filename
end
function ΔandSDPconstraints(prefix::String, G::Group) info(logger, "Generating metric ball of radius $radius...")
info(logger, "Loading precomputed pm, Δ, sdp_constraints...")
pm_fname, Δ_fname = pmΔfilenames(prefix)
product_matrix = load(pm_fname, "pm")
sdp_constraints = constraints(product_matrix)
RG = GroupRing(G, product_matrix)
Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG)
return Δ, sdp_constraints
end
function ΔandSDPconstraints{T<:GroupElem}(name::String, S::Vector{T}, Id::T; radius::Int=2)
info(logger, "Computing pm, Δ, sdp_constraints...")
Δ, sdp_constraints = ΔandSDPconstraints(S, Id, radius=radius)
pm_fname, Δ_fname = pmΔfilenames(name)
save(pm_fname, "pm", parent(Δ).pm)
save(Δ_fname, "Δ", Δ.coeffs)
return Δ, sdp_constraints
end
function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2)
info(logger, "Generating balls of sizes $sizes")
@logtime logger E_R, sizes = Groups.generate_balls(S, Id, radius=2*radius) @logtime logger E_R, sizes = Groups.generate_balls(S, Id, radius=2*radius)
info(logger, "Generated balls of sizes $sizes.")
info(logger, "Creating product matrix...") info(logger, "Creating product matrix...")
@logtime logger pm = GroupRings.create_pm(E_R, GroupRings.reverse_dict(E_R), sizes[radius]; twisted=true) @logtime logger pm = GroupRings.create_pm(E_R, GroupRings.reverse_dict(E_R), sizes[radius]; twisted=true)
info(logger, "Creating sdp_constratints...")
@logtime logger sdp_constraints = PropertyT.constraints(pm)
RG = GroupRing(parent(Id), E_R, pm) RG = GroupRing(parent(Id), E_R, pm)
Δ = splaplacian(RG, S) Δ = spLaplacian(RG, S)
return Δ, sdp_constraints return Δ
end end
function λandP(name::String) function λandP(name::String)
λ_fname, SDP_fname = λSDPfilenames(name) λ_fname = filename(name, :λ)
f₁ = exists(λ_fname) P_fname = filename(name, :P)
f₂ = exists(SDP_fname)
if f₁ && f₂ if exists(λ_fname) && exists(P_fname)
info(logger, "Loading precomputed λ, P...")
λ = load(λ_fname, "λ") λ = load(λ_fname, "λ")
P = load(SDP_fname, "P") P = load(P_fname, "P")
else else
throw(ArgumentError("You need to precompute λ and SDP matrix P to load it!")) throw("You need to precompute $λ_fname and $P_fname to load it!")
end end
return λ, P return λ, P
end end
function λandP(name::String, SDP_problem::JuMP.Model, varλ, varP, warmstart=false) function λandP(name::String, SDP::JuMP.Model, varλ, varP, warmstart=true)
add_handler(solver_logger,
DefaultHandler(joinpath(name, "solver_$(string(now())).log"),
DefaultFormatter("{date}| {msg}")),
"solver_log")
if warmstart && isfile(joinpath(name, "warmstart.jld"))
ws = load(joinpath(name, "warmstart.jld"), "warmstart")
else
ws = nothing
end
λ, P, warmstart = compute_λandP(SDP_problem, varλ, varP, warmstart=ws) if warmstart && isfile(filename(name, :warm))
ws = load(filename(name, :warm), "warmstart")
remove_handler(solver_logger, "solver_log")
λ_fname, P_fname = λSDPfilenames(name)
if λ > 0
save(λ_fname, "λ", λ)
save(P_fname, "P", P)
save(joinpath(name, "warmstart.jld"), "warmstart", warmstart)
else
throw(ErrorException("Solver did not produce a valid solution!: λ = $λ"))
end
return λ, P
end
function fillfrominternal!(m::JuMP.Model, traits)
# Copied from JuMP/src/solvers.jl:178
stat::Symbol = MathProgBase.status(m.internalModel)
numRows, numCols = length(m.linconstr), m.numCols
m.objBound = NaN
m.objVal = NaN
m.colVal = fill(NaN, numCols)
m.linconstrDuals = Array{Float64}(0)
discrete = (traits.int || traits.sos)
if stat == :Optimal
# If we think dual information might be available, try to get it
# If not, return an array of the correct length
if discrete
m.redCosts = fill(NaN, numCols)
m.linconstrDuals = fill(NaN, numRows)
else
if !traits.conic
m.redCosts = try
MathProgBase.getreducedcosts(m.internalModel)[1:numCols]
catch
fill(NaN, numCols)
end
m.linconstrDuals = try
MathProgBase.getconstrduals(m.internalModel)[1:numRows]
catch
fill(NaN, numRows)
end
elseif !traits.qp && !traits.qc
JuMP.fillConicDuals(m)
end
end
else else
# Problem was not solved to optimality, attempt to extract useful ws = nothing
# information anyway
if traits.lin
if stat == :Infeasible
m.linconstrDuals = try
infray = MathProgBase.getinfeasibilityray(m.internalModel)
@assert length(infray) == numRows
infray
catch
suppress_warnings || warn("Infeasibility ray (Farkas proof) not available")
fill(NaN, numRows)
end
elseif stat == :Unbounded
m.colVal = try
unbdray = MathProgBase.getunboundedray(m.internalModel)
@assert length(unbdray) == numCols
unbdray
catch
suppress_warnings || warn("Unbounded ray not available")
fill(NaN, numCols)
end
end
end
# conic duals (currently, SOC and SDP only)
if !discrete && traits.conic && !traits.qp && !traits.qc
if stat == :Infeasible
JuMP.fillConicDuals(m)
end
end
end end
# If the problem was solved, or if it terminated prematurely, try solver_log = solverlogger(name)
# to extract a solution anyway. This commonly occurs when a time
# limit or tolerance is set (:UserLimit) Base.Libc.flush_cstdio()
if !(stat == :Infeasible || stat == :Unbounded) o = redirect_stdout(solver_log.handlers["solver_log"].io)
try Base.Libc.flush_cstdio()
# Do a separate try since getobjval could work while getobjbound does not and vice versa
objBound = MathProgBase.getobjbound(m.internalModel) + m.obj.aff.constant λ, P, warmstart = solve_SDP(SDP, varλ, varP, warmstart=ws)
m.objBound = objBound
end Base.Libc.flush_cstdio()
try redirect_stdout(o)
objVal = MathProgBase.getobjval(m.internalModel) + m.obj.aff.constant
colVal = MathProgBase.getsolution(m.internalModel)[1:numCols] delete!(solver_log.handlers, "solver_log")
# Rescale off-diagonal terms of SDP variables
if traits.sdp if λ > 0
offdiagvars = JuMP.offdiagsdpvars(m) save(filename(name, :λ), "λ", λ)
colVal[offdiagvars] /= sqrt(2) save(filename(name, :P), "P", P)
end save(filename(name, :warm), "warmstart", warmstart)
# Don't corrupt the answers if one of the above two calls fails else
m.objVal = objVal throw(ErrorException("Solver did not produce a valid solution: λ = $λ"))
m.colVal = colVal
end
end end
return λ, P
return stat
end
function compute_λandP(m, varλ, varP; warmstart=nothing)
λ = 0.0
P = nothing
traits = JuMP.ProblemTraits(m, relaxation=false)
while λ == 0.0
try
JuMP.build(m, traits=traits)
if warmstart != nothing
p_sol, d_sol, s = warmstart
MathProgBase.SolverInterface.setwarmstart!(m.internalModel, p_sol; dual_sol = d_sol, slack=s);
end
solve_SDP(m)
λ = MathProgBase.getobjval(m.internalModel)
catch y
warn(solver_logger, y)
end
end
warmstart = (m.internalModel.primal_sol, m.internalModel.dual_sol,
m.internalModel.slack)
fillfrominternal!(m, traits)
P = JuMP.getvalue(varP)
λ = JuMP.getvalue(varλ)
return λ, P, warmstart
end end
Kazhdan_from_sgap(λ,N) = sqrt(2*λ/N) Kazhdan_from_sgap(λ,N) = sqrt(2*λ/N)
function check_λ(name, S, λ, P, radius, logger)
RG = GroupRing(parent(first(S)), load(filename(name, :pm), "pm"))
Δ = GroupRingElem(load(filename(name, :Δ), "Δ")[:, 1], RG)
@logtime logger Q = real(sqrtm(Symmetric(P)))
sgap = check_distance_to_cone(Δ, λ, Q, 2*radius, logger)
if sgap > 0
info(logger, "λ($name, S) ≥ $(Float64(trunc(sgap,12)))")
Kazhdan_κ = Kazhdan_from_sgap(sgap, length(S))
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!")
return true
else
sgap = Float64(trunc(sgap, 12))
info(logger, "λ($name, S) ≥ $sgap: Group may NOT HAVE property (T)!")
return false
end
end
function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius) function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius)
isdir(name) || mkdir(name) isdir(name) || mkdir(name)
LOGGER = Memento.getlogger()
if all(exists.(pmΔfilenames(name))) if exists(filename(name, :pm)) && exists(filename(name, :Δ))
# cached # cached
Δ, sdp_constraints = ΔandSDPconstraints(name, parent(S[1])) info(LOGGER, "Loading precomputed Δ...")
Δ = Laplacian(name, parent(S[1]))
else else
# compute # compute
Δ, sdp_constraints = ΔandSDPconstraints(name, S, Id, radius=radius) Δ = Laplacian(S, Id, LOGGER, radius=radius)
save(filename(name, :pm), "pm", parent(Δ).pm)
save(filename(name, :Δ), "Δ", Δ.coeffs)
end end
if all(exists.(λSDPfilenames(name))) fullpath = joinpath(name, string(upper_bound))
λ, P = λandP(name) isdir(fullpath) || mkdir(fullpath)
else
info(logger, "Creating SDP problem...")
SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound)
JuMP.setsolver(SDP_problem, solver)
if exists(filename(fullpath, :λ)) && exists(filename(fullpath, :P))
info(LOGGER, "Loading precomputed λ, P...")
λ, P = λandP(fullpath)
else
info(LOGGER, "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))
λ, P = λandP(name, SDP_problem, λ, P) @logtime LOGGER λ, P = λandP(fullpath, SDP_problem, varλ, varP)
end end
info(logger, "λ = $λ") info(LOGGER, "λ = $λ")
info(logger, "sum(P) = $(sum(P))") info(LOGGER, "sum(P) = $(sum(P))")
info(logger, "maximum(P) = $(maximum(P))") info(LOGGER, "maximum(P) = $(maximum(P))")
info(logger, "minimum(P) = $(minimum(P))") info(LOGGER, "minimum(P) = $(minimum(P))")
if λ > 0 isapprox(eigvals(P), abs.(eigvals(P)), atol=tol) ||
pm_fname, Δ_fname = pmΔfilenames(name) warn("The solution matrix doesn't seem to be positive definite!")
RG = GroupRing(parent(first(S)), load(pm_fname, "pm"))
Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG)
isapprox(eigvals(P), abs(eigvals(P)), atol=tol) || if λ > 0
warn("The solution matrix doesn't seem to be positive definite!") return check_λ(name, S, λ, P, radius, LOGGER)
# @assert P == Symmetric(P) end
@logtime logger Q = real(sqrtm(Symmetric(P))) info(LOGGER, "κ($name, S) ≥ $λ < 0: Tells us nothing about property (T)")
return false
sgap = distance_to_positive_cone(Δ, λ, Q, 2*radius)
if isa(sgap, Interval)
sgap = sgap.lo
end
if sgap > 0
info(logger, "λ ≥ $(Float64(trunc(sgap,12)))")
Kazhdan_κ = Kazhdan_from_sgap(sgap, length(S))
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!")
return true
else
sgap = Float64(trunc(sgap, 12))
info(logger, "λ($name, S) ≥ $sgap: Group may NOT HAVE property (T)!")
return false
end
end
info(logger, "κ($name, S) ≥ $λ < 0: Tells us nothing about property (T)")
return false
end end
include("SDPs.jl") include("SDPs.jl")

131
src/SDPs.jl
View File

@ -13,7 +13,7 @@ function constraints(pm, total_length=maximum(pm))
return constraints return constraints
end end
function splaplacian(RG::GroupRing, S, T::Type=Float64) function spLaplacian(RG::GroupRing, S, T::Type=Float64)
result = RG(T) result = RG(T)
result[RG.group()] = T(length(S)) result[RG.group()] = T(length(S))
for s in S for s in S
@ -22,7 +22,7 @@ function splaplacian(RG::GroupRing, S, T::Type=Float64)
return result return result
end end
function splaplacian{TT<:Ring}(RG::GroupRing{TT}, S, T::Type=Float64) function spLaplacian{TT<:Ring}(RG::GroupRing{TT}, S, T::Type=Float64)
result = RG(T) result = RG(T)
result[one(RG.group)] = T(length(S)) result[one(RG.group)] = T(length(S))
for s in S for s in S
@ -55,21 +55,122 @@ function create_SDP_problem(Δ::GroupRingElem, matrix_constraints; upper_bound=I
return m, λ, P return m, λ, P
end end
function solve_SDP(SDP_problem) function solve_SDP(m, varλ, varP; warmstart=nothing)
info(logger, Base.repr(SDP_problem))
o = redirect_stdout(solver_logger.handlers["solver_log"].io) traits = JuMP.ProblemTraits(m, relaxation=false)
Base.Libc.flush_cstdio()
@logtime logger solution_status = MathProgBase.optimize!(SDP_problem.internalModel) JuMP.build(m, traits=traits)
Base.Libc.flush_cstdio() if warmstart != nothing
p_sol, d_sol, s = warmstart
redirect_stdout(o) MathProgBase.SolverInterface.setwarmstart!(m.internalModel, p_sol; dual_sol = d_sol, slack=s);
if solution_status != :Optimal
warn(logger, "The solver did not solve the problem successfully!")
end end
info(logger, solution_status)
return 0 MathProgBase.optimize!(m.internalModel)
λ = MathProgBase.getobjval(m.internalModel)
warmstart = (m.internalModel.primal_sol, m.internalModel.dual_sol,
m.internalModel.slack)
fillfrominternal!(m, traits)
P = JuMP.getvalue(varP)
λ = JuMP.getvalue(varλ)
return λ, P, warmstart
end
function fillfrominternal!(m::JuMP.Model, traits)
# Copied from JuMP/src/solvers.jl:178
stat::Symbol = MathProgBase.status(m.internalModel)
numRows, numCols = length(m.linconstr), m.numCols
m.objBound = NaN
m.objVal = NaN
m.colVal = fill(NaN, numCols)
m.linconstrDuals = Array{Float64}(0)
discrete = (traits.int || traits.sos)
if stat == :Optimal
# If we think dual information might be available, try to get it
# If not, return an array of the correct length
if discrete
m.redCosts = fill(NaN, numCols)
m.linconstrDuals = fill(NaN, numRows)
else
if !traits.conic
m.redCosts = try
MathProgBase.getreducedcosts(m.internalModel)[1:numCols]
catch
fill(NaN, numCols)
end
m.linconstrDuals = try
MathProgBase.getconstrduals(m.internalModel)[1:numRows]
catch
fill(NaN, numRows)
end
elseif !traits.qp && !traits.qc
JuMP.fillConicDuals(m)
end
end
else
# Problem was not solved to optimality, attempt to extract useful
# information anyway
if traits.lin
if stat == :Infeasible
m.linconstrDuals = try
infray = MathProgBase.getinfeasibilityray(m.internalModel)
@assert length(infray) == numRows
infray
catch
suppress_warnings || warn("Infeasibility ray (Farkas proof) not available")
fill(NaN, numRows)
end
elseif stat == :Unbounded
m.colVal = try
unbdray = MathProgBase.getunboundedray(m.internalModel)
@assert length(unbdray) == numCols
unbdray
catch
suppress_warnings || warn("Unbounded ray not available")
fill(NaN, numCols)
end
end
end
# conic duals (currently, SOC and SDP only)
if !discrete && traits.conic && !traits.qp && !traits.qc
if stat == :Infeasible
JuMP.fillConicDuals(m)
end
end
end
# If the problem was solved, or if it terminated prematurely, try
# to extract a solution anyway. This commonly occurs when a time
# limit or tolerance is set (:UserLimit)
if !(stat == :Infeasible || stat == :Unbounded)
try
# Do a separate try since getobjval could work while getobjbound does not and vice versa
objBound = MathProgBase.getobjbound(m.internalModel) + m.obj.aff.constant
m.objBound = objBound
end
try
objVal = MathProgBase.getobjval(m.internalModel) + m.obj.aff.constant
colVal = MathProgBase.getsolution(m.internalModel)[1:numCols]
# Rescale off-diagonal terms of SDP variables
if traits.sdp
offdiagvars = JuMP.offdiagsdpvars(m)
colVal[offdiagvars] /= sqrt(2)
end
# Don't corrupt the answers if one of the above two calls fails
m.objVal = objVal
m.colVal = colVal
end
end
return stat
end end