GroupsWithPropertyT/Orb_AutFN.jl

207 lines
5.7 KiB
Julia
Raw Normal View History

2017-06-06 12:15:30 +02:00
using JLD
using JuMP
using SCS
using GroupRings
using PropertyT
2017-06-06 18:47:10 +02:00
import Nemo: Group, GroupElem
2017-06-06 12:15:30 +02:00
using ArgParse
2017-06-06 17:54:05 +02:00
include("OrbitDecomposition.jl")
2017-06-08 21:38:14 +02:00
immutable OrbitData
2017-06-06 12:15:30 +02:00
name::String
2017-06-08 21:37:15 +02:00
Us::Vector
2017-06-06 12:15:30 +02:00
Ps::Vector{Array{JuMP.Variable,2}}
2017-06-08 21:37:15 +02:00
cnstr::Vector
laplacian::Vector
laplacianSq::Vector
2017-06-06 12:15:30 +02:00
dims::Vector{Int}
end
2017-06-08 21:35:27 +02:00
immutable Settings
name::String
N::Int
G::Group
S::Vector
AutS::Group
radius::Int
solver::SCSSolver
upper_bound::Float64
tol::Float64
end
2017-06-06 12:15:30 +02:00
function sparsify!{T}(U::AbstractArray{T}, eps=eps(T))
# n = rank(U)
U[abs.(U) .< eps] = zero(T)
# @assert rank(U) == n
return sparse(U)
end
sparsify{T}(U::AbstractArray{T}, eps=eps(T)) = sparsify!(deepcopy(U), eps)
small_to_zero!{T}(A::AbstractArray{T}, eps=eps(T)) = A[abs(A) .< eps] = zero(T)
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, P[k][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
2017-06-08 21:38:14 +02:00
function init_OrbitData(name::String)
splap = load(joinpath(name, "delta.jld"), "Δ");
2017-06-06 12:15:30 +02:00
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);
2017-06-06 17:53:33 +02:00
#dimensions of the corresponding πs:
2017-06-06 18:48:25 +02:00
dims = load(joinpath(name, "U_pis.jld"), "dims")
2017-06-06 12:15:30 +02:00
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)
2017-06-08 21:38:14 +02:00
orbData = OrbitData(name, Uπs, P, orb_spcnstrm, orb_splap, orb_splap², dims);
2017-06-06 12:15:30 +02:00
2017-06-08 21:38:14 +02:00
return m, orbData
2017-06-06 12:15:30 +02:00
end
function transform{T}(U::AbstractArray{T,2}, V::AbstractArray{T,2}, eps=eps(T))
w = U'*V*U
sparsify!(w, eps)
dropzeros!(w)
return w
end
2017-06-08 21:38:14 +02:00
A(data::OrbitData, π, t) = data.dims[π]*transform(data.Us[π], data.cnstr[t])
2017-06-06 12:15:30 +02:00
2017-06-08 21:38:14 +02:00
function constrLHS(m::JuMP.Model, data::OrbitData, t)
2017-06-06 12:15:30 +02:00
l = endof(data.Us)
lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l))
return lhs
end
2017-06-08 21:38:14 +02:00
function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.cnstr); var::Symbol = )
2017-06-06 12:15:30 +02:00
λ = getvariable(m, var)
for t in 1:l
d, = data.laplacian[t], data.laplacianSq[t]
lhs = constrLHS(m, data, t)
if lhs == zero(lhs)
if d == 0 && == 0
info("Detected empty constraint")
continue
else
warn("Adding unsatisfiable constraint!")
end
end
JuMP.@constraint(m, lhs == - λ*d)
end
end
2017-06-08 21:38:14 +02:00
function reconstructP(m::JuMP.Model, data::OrbitData)
2017-06-06 12:15:30 +02:00
computedPs = [getvalue(P) for P in data.Ps]
return sum(data.dims[π]*data.Us[π]*computedPs[π]*data.Us[π]' for π in 1:endof(data.Ps))
end
function create_SDP_problem(name::String; upper_bound=Inf)
info(PropertyT.logger, "Loading data....")
2017-06-08 21:38:14 +02:00
t = @timed SDP_problem, orb_data = init_OrbitData(name);
2017-06-06 12:15:30 +02:00
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
2017-06-08 21:38:14 +02:00
function λandP(m::JuMP.Model, data::OrbitData)
2017-06-06 12:15:30 +02:00
info(PropertyT.logger, "Solving SDP problem...")
varλ = JuMP.getvariable(m, )
varP = data.Ps
λ, P = PropertyT.λandP(data.name, m, varλ, varP)
recP = reconstructP(m, data)
fname = PropertyT.λSDPfilenames(data.name)[2]
save(fname, "origP", P, "P", recP)
return λ, recP
end
2017-06-08 21:39:45 +02:00
function init_orbit_data(logger, sett::Settings; radius=2)
2017-06-08 21:39:45 +02:00
ex(fname) = isfile(joinpath(sett.name, fname))
files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld"])
if !all(files_exists)
2017-06-08 21:39:45 +02:00
compute_orbit_data(logger, sett.name, sett.G, sett.S, sett.AutS, radius=radius)
end
return 0
end
2017-06-08 21:35:27 +02:00
function orbit_check_propertyT(logger, sett::Settings)
2017-06-06 12:15:30 +02:00
2017-06-08 21:35:27 +02:00
init_orbit_data(logger, sett, radius=sett.radius)
2017-06-08 21:35:27 +02:00
Δ = float(PropertyT.ΔandSDPconstraints(sett.name, sett.G)[1])
2017-06-06 12:15:30 +02:00
2017-06-08 21:35:27 +02:00
fnames = PropertyT.λSDPfilenames(sett.name)
2017-06-08 20:07:34 +02:00
2017-06-06 18:02:54 +02:00
if all(isfile.(fnames))
2017-06-08 21:35:27 +02:00
λ, P = PropertyT.λandP(sett.name)
2017-06-06 12:15:30 +02:00
else
info(logger, "Creating SDP problem...")
2017-06-08 21:35:27 +02:00
SDP_problem, orb_data = create_SDP_problem(sett.name, upper_bound=sett.upper_bound)
JuMP.setsolver(SDP_problem, sett.solver)
2017-06-06 12:15:30 +02:00
λ, P = λandP(SDP_problem, orb_data)
end
info(PropertyT.logger, "λ = ")
info(PropertyT.logger, "sum(P) = $(sum(P))")
info(PropertyT.logger, "maximum(P) = $(maximum(P))")
info(PropertyT.logger, "minimum(P) = $(minimum(P))")
if λ > 0
2017-06-08 21:35:27 +02:00
sgap = PropertyT.check_distance_to_positive_cone(Δ, λ, P, tol=sett.tol, rational=false, len=2*sett.radius)
if isa(sgap, Interval)
sgap = sgap.lo
end
if sgap > 0
info(logger, "λ ≥ $(Float64(trunc(sgap,12)))")
Kazhdan_κ = PropertyT.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
2017-06-06 12:15:30 +02:00
end