using JLD using JuMP using SCS using GroupRings using PropertyT using ValidatedNumerics import Nemo: Group, GroupElem using ArgParse include("OrbitDecomposition.jl") immutable OrbitData name::String Us::Vector Ps::Vector{Array{JuMP.Variable,2}} cnstr::Vector laplacian::Vector laplacianSq::Vector dims::Vector{Int} end immutable Settings name::String N::Int G::Group S::Vector AutS::Group radius::Int solver::SCSSolver upper_bound::Float64 tol::Float64 end 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 U = sparse(U) dropzeros!(U) return U end function sparsify!{T}(U::SparseMatrixCSC{T}, eps=eps(T)) U[abs.(U) .< eps] = zero(T) dropzeros!(U) return 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, [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); return m, orbData end function transform{T}(U::AbstractArray{T,2}, V::AbstractArray{T,2}, eps=eps(T); sparse=true) w = U'*V*U sparse && sparsify!(w, eps) 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.cnstr); var::Symbol = :λ) λ = m[var] for t in 1:l d, d² = data.laplacian[t], data.laplacianSq[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 reconstructP(m::JuMP.Model, data::OrbitData) 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....") 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 λ, 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 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) 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