250 lines
6.8 KiB
Julia
250 lines
6.8 KiB
Julia
using JLD
|
|
using JuMP
|
|
using SCS
|
|
|
|
using GroupRings
|
|
using PropertyT
|
|
|
|
using ValidatedNumerics
|
|
using ArgParse
|
|
|
|
import Nemo: Group, GroupElem
|
|
|
|
immutable Settings
|
|
name::String
|
|
N::Int
|
|
G::Group
|
|
S::Vector
|
|
AutS::Group
|
|
radius::Int
|
|
solver::SCSSolver
|
|
upper_bound::Float64
|
|
tol::Float64
|
|
end
|
|
|
|
immutable OrbitData
|
|
name::String
|
|
Us::Vector
|
|
Ps::Vector{Array{JuMP.Variable,2}}
|
|
cnstr::Vector
|
|
laplacian::Vector
|
|
laplacianSq::Vector
|
|
dims::Vector{Int}
|
|
end
|
|
|
|
include("OrbitDecomposition.jl")
|
|
function sparsify{T}(U::Array{T}, eps=eps(T))
|
|
n = rank(U)
|
|
W = deepcopy(U)
|
|
W[abs.(W) .< eps] = zero(T)
|
|
if rank(W) != n
|
|
warn("Sparsification would decrease the rank!")
|
|
W = U
|
|
end
|
|
W = sparse(W)
|
|
dropzeros!(W)
|
|
return W
|
|
end
|
|
|
|
function sparsify!{T}(U::SparseMatrixCSC{T}, eps=eps(T))
|
|
U[abs.(U) .< eps] = zero(T)
|
|
dropzeros!(U)
|
|
return U
|
|
end
|
|
|
|
sparsify{T}(U::SparseMatrixCSC{T}, eps=eps(T)) = sparsify!(deepcopy(U), eps)
|
|
|
|
function init_model(Uπs)
|
|
m = JuMP.Model();
|
|
l = size(Uπs,1)
|
|
P = Vector{Array{JuMP.Variable,2}}(l)
|
|
|
|
for k in 1:l
|
|
s = size(Uπs[k],2)
|
|
P[k] = JuMP.@variable(m, [i=1:s, j=1:s])
|
|
JuMP.@SDconstraint(m, P[k] >= 0.0)
|
|
end
|
|
|
|
JuMP.@variable(m, λ >= 0.0)
|
|
JuMP.@objective(m, Max, λ)
|
|
return m, P
|
|
end
|
|
|
|
function init_OrbitData(name::String)
|
|
splap = load(joinpath(name, "delta.jld"), "Δ");
|
|
pm = load(joinpath(name, "pm.jld"), "pm");
|
|
cnstr = PropertyT.constraints_from_pm(pm);
|
|
splap² = GroupRings.mul(splap, splap, pm);
|
|
|
|
Uπs = load(joinpath(name, "U_pis.jld"), "Uπs");
|
|
# Uπs = sparsify.(Uπs);
|
|
#dimensions of the corresponding πs:
|
|
dims = load(joinpath(name, "U_pis.jld"), "dims")
|
|
|
|
m, P = init_model(Uπs);
|
|
|
|
orbits = load(joinpath(name, "orbits.jld"), "orbits");
|
|
n = size(Uπs[1],1)
|
|
orb_spcnstrm = [orbit_constraint(cnstr[collect(orb)], n) for orb in orbits]
|
|
orb_splap = orbit_spvector(splap, orbits)
|
|
orb_splap² = orbit_spvector(splap², orbits)
|
|
|
|
orbData = OrbitData(name, Uπs, P, orb_spcnstrm, orb_splap, orb_splap², dims);
|
|
|
|
# orbData = OrbitData(name, Uπs, P, orb_spcnstrm, splap, splap², dims);
|
|
|
|
return m, orbData
|
|
end
|
|
|
|
function transform(U::AbstractArray, V::AbstractArray; sparse=false)
|
|
w = U'*V*U
|
|
sparse && sparsify!(w)
|
|
return w
|
|
end
|
|
|
|
A(data::OrbitData, π, t) = data.dims[π]*transform(data.Us[π], data.cnstr[t])
|
|
|
|
function constrLHS(m::JuMP.Model, data::OrbitData, t)
|
|
l = endof(data.Us)
|
|
lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l))
|
|
return lhs
|
|
end
|
|
|
|
function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.laplacian); var::Symbol = :λ)
|
|
λ = m[var]
|
|
# orbits = load(joinpath(data.name, "orbits.jld"), "orbits");
|
|
# locate(t, orb=orbits) = findfirst(x->t in x, orb)
|
|
for t in 1:l
|
|
d, d² = data.laplacian[t], data.laplacianSq[t]
|
|
# lhs = constrLHS(m, data, locate(t))
|
|
lhs = constrLHS(m, data, t)
|
|
if lhs == zero(lhs)
|
|
if d == 0 && d² == 0
|
|
info("Detected empty constraint")
|
|
continue
|
|
else
|
|
warn("Adding unsatisfiable constraint!")
|
|
end
|
|
end
|
|
JuMP.@constraint(m, lhs == d² - λ*d)
|
|
end
|
|
end
|
|
|
|
function create_SDP_problem(name::String; upper_bound=Inf)
|
|
info(PropertyT.logger, "Loading data....")
|
|
t = @timed SDP_problem, orb_data = init_OrbitData(name);
|
|
info(PropertyT.logger, PropertyT.timed_msg(t))
|
|
|
|
if upper_bound < Inf
|
|
λ = JuMP.getvariable(SDP_problem, :λ)
|
|
JuMP.@constraint(SDP_problem, λ <= upper_bound)
|
|
end
|
|
|
|
info(PropertyT.logger, "Adding constraints... ")
|
|
t = @timed addconstraints!(SDP_problem, orb_data)
|
|
info(PropertyT.logger, PropertyT.timed_msg(t))
|
|
|
|
return SDP_problem, orb_data
|
|
end
|
|
|
|
function λandP(m::JuMP.Model, data::OrbitData)
|
|
info(PropertyT.logger, "Solving SDP problem...")
|
|
varλ = m[:λ]
|
|
varP = data.Ps
|
|
λ, Ps = PropertyT.λandP(data.name, m, varλ, varP)
|
|
return λ, Ps
|
|
end
|
|
|
|
function matrix_reps{T<:GroupElem}(G::Group, S::Vector{T}, AutS::Group, radius::Int)
|
|
Id = (isa(G, Nemo.Ring) ? one(G) : G())
|
|
E2, _ = Groups.generate_balls(S, Id, radius=radius)
|
|
Edict = GroupRings.reverse_dict(E2)
|
|
|
|
mreps = Dict(g => matrix_repr(g, E2, Edict) for g in elements(AutS))
|
|
return mreps
|
|
end
|
|
|
|
function reconstruct(mreps::Dict, Us::Vector, Ps::Vector, dims::Vector)
|
|
recP = zeros(size(Us[1],1), size(Us[1],1))
|
|
|
|
for g in keys(mreps)
|
|
for π in 1:endof(Us)
|
|
recP .+= dims[π] .* mreps[g]*transform(Us[π]', Ps[π])*mreps[inv(g)]
|
|
end
|
|
end
|
|
recP .*= 1/length(collect(keys(mreps)))
|
|
|
|
return recP
|
|
end
|
|
|
|
function λandP(m::JuMP.Model, data::OrbitData, sett::Settings)
|
|
info(PropertyT.logger, "Solving SDP problem...")
|
|
λ, Ps = λandP(m, data)
|
|
|
|
info(PropertyT.logger, "Reconstructing P...")
|
|
mreps = matrix_reps(sett.G, sett.S, sett.AutS, sett.radius)
|
|
|
|
recP = reconstruct(mreps, data.Us, Ps, data.dims)
|
|
|
|
fname = PropertyT.λSDPfilenames(data.name)[2]
|
|
save(fname, "origP", Ps, "P", recP)
|
|
return λ, recP
|
|
end
|
|
|
|
function init_orbit_data(logger, sett::Settings; radius=2)
|
|
|
|
ex(fname) = isfile(joinpath(sett.name, fname))
|
|
|
|
files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld"])
|
|
|
|
if !all(files_exists)
|
|
compute_orbit_data(logger, sett.name, sett.G, sett.S, sett.AutS, radius=radius)
|
|
end
|
|
|
|
return 0
|
|
end
|
|
|
|
function orbit_check_propertyT(logger, sett::Settings)
|
|
|
|
init_orbit_data(logger, sett, radius=sett.radius)
|
|
|
|
Δ = PropertyT.ΔandSDPconstraints(sett.name, sett.G)[1]
|
|
|
|
fnames = PropertyT.λSDPfilenames(sett.name)
|
|
|
|
if all(isfile.(fnames))
|
|
λ, P = PropertyT.λandP(sett.name)
|
|
else
|
|
info(logger, "Creating SDP problem...")
|
|
SDP_problem, orb_data = create_SDP_problem(sett.name, upper_bound=sett.upper_bound)
|
|
JuMP.setsolver(SDP_problem, sett.solver)
|
|
|
|
λ, P = λandP(SDP_problem, orb_data, sett)
|
|
end
|
|
|
|
info(logger, "λ = $λ")
|
|
info(logger, "sum(P) = $(sum(P))")
|
|
info(logger, "maximum(P) = $(maximum(P))")
|
|
info(logger, "minimum(P) = $(minimum(P))")
|
|
|
|
if λ > 0
|
|
sgap = PropertyT.check_distance_to_positive_cone(Δ, λ, P, 2*sett.radius, tol=sett.tol, rational=false)
|
|
if isa(sgap, Interval)
|
|
sgap = sgap.lo
|
|
end
|
|
if sgap > 0
|
|
info(logger, "λ ≥ $(Float64(trunc(sgap,12)))")
|
|
Kazhdan_κ = PropertyT.Kazhdan_from_sgap(sgap, length(sett.S))
|
|
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
|
|
info(logger, "κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
|
return true
|
|
else
|
|
sgap = Float64(trunc(sgap, 12))
|
|
info(logger, "λ($(sett.name), S) ≥ $sgap: Group may NOT HAVE property (T)!")
|
|
return false
|
|
end
|
|
end
|
|
info(logger, "κ($(sett.name), S) ≥ $λ < 0: Tells us nothing about property (T)")
|
|
return false
|
|
end
|