diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index dc2aee7..768397c 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -23,80 +23,80 @@ end function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int) - # result = zeros(eltype(Q), l) - # r = similar(result) - # for i in 1:size(Q,2) - # print(" $i") - # result += GroupRings.mul!(r, view(Q,:,i), view(Q,:,i), pm) - # end + # result = zeros(eltype(Q), l) + # r = similar(result) + # for i in 1:size(Q,2) + # print(" $i") + # result += GroupRings.mul!(r, view(Q,:,i), view(Q,:,i), pm) + # end - @everywhere groupring_square = PropertyT.groupring_square + @everywhere groupring_square = PropertyT.groupring_square - result = @parallel (+) for i in 1:size(Q,2) - groupring_square(Q[:,i], l, pm) - end + result = @parallel (+) for i in 1:size(Q,2) + groupring_square(Q[:,i], l, pm) +end - println("") +println("") - return result +return result end function compute_SOS(Q::AbstractArray, RG::GroupRing, l::Int) - result = compute_SOS(Q, RG.pm, l) - return GroupRingElem(result, RG) + result = compute_SOS(Q, RG.pm, l) + return GroupRingElem(result, RG) end function distance_to_cone{T<:Interval}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int) - SOS = compute_SOS(Q, parent(elt), length(elt.coeffs)) - SOS_diff = elt - SOS + SOS = compute_SOS(Q, parent(elt), length(elt.coeffs)) + SOS_diff = elt - SOS - ɛ_dist = GroupRings.augmentation(SOS_diff) - info(LOGGER, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)") + ɛ_dist = GroupRings.augmentation(SOS_diff) + info(LOGGER, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)") - eoi_SOS_L1_dist = norm(SOS_diff,1) - info(LOGGER, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)") + eoi_SOS_L1_dist = norm(SOS_diff,1) + info(LOGGER, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)") - dist = 2^(wlen-1)*eoi_SOS_L1_dist - return 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 + SOS = compute_SOS(Q, parent(elt), length(elt.coeffs)) + SOS_diff = elt - SOS - ɛ_dist = GroupRings.augmentation(SOS_diff) - info(LOGGER, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))") + ɛ_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))") + 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 + dist = 2^(wlen-1)*eoi_SOS_L1_dist + return dist end function augIdproj{T, I<:AbstractInterval}(S::Type{I}, Q::AbstractArray{T,2}) - l = size(Q, 2) - R = zeros(S, (l,l)) - Threads.@threads for j in 1:l - col = sum(view(Q, :,j))/l - for i in 1:l - R[i,j] = Q[i,j] - col ± eps(0.0) - end - end - return R + l = size(Q, 2) + R = zeros(S, (l,l)) + Threads.@threads for j in 1:l + col = sum(view(Q, :,j))/l + for i in 1:l + R[i,j] = Q[i,j] - col ± eps(0.0) + end + end + return R end function augIdproj{T}(Q::AbstractArray{T,2}, logger) - info(logger, "Projecting columns of Q to the augmentation ideal...") - @logtime logger Q = augIdproj(Interval{T}, Q) + info(logger, "Projecting columns of Q to the augmentation ideal...") + @logtime logger Q = augIdproj(Interval{T}, Q) - 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)]) - info(logger, (check? "They do." : "FAILED!")) + 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)]) + info(logger, (check? "They do." : "FAILED!")) - @assert check + @assert check - return Q + return Q end function distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen::Int, logger) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 71c25f0..352c4cf 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -4,16 +4,16 @@ using SCS export Settings, OrbitData immutable Settings{T<:AbstractMathProgSolver} - name::String - N::Int - G::Group - S::Vector - autS::Group - radius::Int - solver::T - upper_bound::Float64 - tol::Float64 - warmstart::Bool + name::String + N::Int + G::Group + S::Vector + autS::Group + radius::Int + solver::T + upper_bound::Float64 + tol::Float64 + warmstart::Bool end prefix(s::Settings) = s.name @@ -22,43 +22,43 @@ prepath(s::Settings) = prefix(s) fullpath(s::Settings) = joinpath(prefix(s), suffix(s)) immutable OrbitData{T<:AbstractArray{Float64, 2}, LapType <:AbstractVector{Float64}} - name::String - Us::Vector{T} - Ps::Vector{Array{JuMP.Variable,2}} - cnstr::Vector{SparseMatrixCSC{Float64, Int}} - laplacian::LapType - laplacianSq::LapType - dims::Vector{Int} + name::String + Us::Vector{T} + Ps::Vector{Array{JuMP.Variable,2}} + cnstr::Vector{SparseMatrixCSC{Float64, Int}} + laplacian::LapType + laplacianSq::LapType + dims::Vector{Int} end function OrbitData(sett::Settings) - splap = load(joinpath(prepath(sett), "delta.jld"), "Δ"); - pm = load(joinpath(prepath(sett), "pm.jld"), "pm"); - cnstr = PropertyT.constraints(pm); - splap² = similar(splap) - splap² = GroupRings.mul!(splap², splap, splap, pm); + splap = load(joinpath(prepath(sett), "delta.jld"), "Δ"); + pm = load(joinpath(prepath(sett), "pm.jld"), "pm"); + cnstr = PropertyT.constraints(pm); + splap² = similar(splap) + splap² = GroupRings.mul!(splap², splap, splap, pm); - Uπs = load(joinpath(prepath(sett), "U_pis.jld"), "Uπs") - nzros = [i for i in 1:length(Uπs) if size(Uπs[i],2) !=0] - Uπs = Uπs[nzros] - Uπs = sparsify!.(Uπs, sett.tol, check=true, verbose=true) + Uπs = load(joinpath(prepath(sett), "U_pis.jld"), "Uπs") + nzros = [i for i in 1:length(Uπs) if size(Uπs[i],2) !=0] + Uπs = Uπs[nzros] + Uπs = sparsify!.(Uπs, sett.tol, check=true, verbose=true) - #dimensions of the corresponding πs: - dims = load(joinpath(prepath(sett), "U_pis.jld"), "dims")[nzros] + #dimensions of the corresponding πs: + dims = load(joinpath(prepath(sett), "U_pis.jld"), "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"); - 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) + orbits = load(joinpath(prepath(sett), "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(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 include("OrbitDecomposition.jl") @@ -67,152 +67,152 @@ dens(M::SparseMatrixCSC) = length(M.nzval)/length(M) dens(M::AbstractArray) = length(findn(M)[1])/length(M) function sparsify!{Tv,Ti}(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv); verbose=false) - n = nnz(M) + n = nnz(M) - densM = dens(M) - for i in eachindex(M.nzval) - if abs(M.nzval[i]) < eps - M.nzval[i] = zero(Tv) - end - end - dropzeros!(M) - m = nnz(M) + densM = dens(M) + for i in eachindex(M.nzval) + if abs(M.nzval[i]) < eps + M.nzval[i] = zero(Tv) + end + end + dropzeros!(M) + m = nnz(M) - if verbose - info(LOGGER, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M), 20)) - end + if verbose + info(LOGGER, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M), 20)) + end - return M + return M end function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); check=false, verbose=false) - densM = dens(M) - rankM = rank(M) - M[abs.(M) .< eps] .= zero(T) + densM = dens(M) + rankM = rank(M) + M[abs.(M) .< eps] .= zero(T) - if check && rankM != rank(M) - warn(LOGGER, "Sparsification decreased the rank!") - end + if check && rankM != rank(M) + warn(LOGGER, "Sparsification decreased the rank!") + end - if verbose - info(LOGGER, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M),20)) - end + if verbose + info(LOGGER, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M),20)) + end - return sparse(M) + return sparse(M) end 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)) + ex(fname) = isfile(joinpath(prepath(sett), fname)) - files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld", "preps.jld"]) + 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 + if !all(files_exists) + compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.autS, radius=radius) + end - return 0 + return 0 end function transform(U::AbstractArray, V::AbstractArray; sparse=true) - if sparse - return sparsify!(U'*V*U) - else - return U'*V*U - end + if sparse + return sparsify!(U'*V*U) + else + return U'*V*U + end 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 + l = endof(data.Us) + lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l)) + return lhs end 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)] - return @expression(m, sum(vecdot(M[π], vars[π]) 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))) end function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.laplacian); var::Symbol=:λ) - λ = m[var] - Ust = [U' for U in data.Us] - idx = [π for π in 1:endof(data.Us) if size(data.Us[π],2) != 0] + λ = m[var] + Ust = [U' for U in data.Us] + idx = [π for π in 1:endof(data.Us) if size(data.Us[π],2) != 0] - for t in 1:l - if t % 100 == 0 - print(t, ", ") - end - # lhs = constrLHS(m, data, t) - lhs = constrLHS(m, data.cnstr[t], data.Us[idx], Ust[idx], data.dims[idx], data.Ps[idx]) + for t in 1:l + if t % 100 == 0 + print(t, ", ") + end + # lhs = constrLHS(m, data, t) + 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] - # 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 - println("") + d, d² = data.laplacian[t], data.laplacianSq[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 + println("") end function init_model(n, sizes) - m = JuMP.Model(); - P = Vector{Array{JuMP.Variable,2}}(n) + m = JuMP.Model(); + P = Vector{Array{JuMP.Variable,2}}(n) - for (k,s) in enumerate(sizes) - P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) - JuMP.@SDconstraint(m, P[k] >= 0.0) - end + for (k,s) in enumerate(sizes) + 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 + JuMP.@variable(m, λ >= 0.0) + JuMP.@objective(m, Max, λ) + return m, P end function create_SDP_problem(sett::Settings) - info(LOGGER, "Loading orbit data....") - @logtime LOGGER SDP_problem, orb_data = OrbitData(sett); + info(LOGGER, "Loading orbit data....") + @logtime LOGGER SDP_problem, orb_data = OrbitData(sett); - if sett.upper_bound < Inf - λ = JuMP.getvariable(SDP_problem, :λ) - JuMP.@constraint(SDP_problem, λ <= sett.upper_bound) - end + if sett.upper_bound < Inf + λ = JuMP.getvariable(SDP_problem, :λ) + JuMP.@constraint(SDP_problem, λ <= sett.upper_bound) + end - t = length(orb_data.laplacian) - info(LOGGER, "Adding $t constraints ... ") - @logtime LOGGER addconstraints!(SDP_problem, orb_data) + t = length(orb_data.laplacian) + info(LOGGER, "Adding $t constraints ... ") + @logtime LOGGER addconstraints!(SDP_problem, orb_data) - return SDP_problem, orb_data + return SDP_problem, orb_data end function λandP(m::JuMP.Model, data::OrbitData, warmstart=true) - varλ = m[:λ] - varP = data.Ps - λ, Ps = PropertyT.λandP(data.name, m, varλ, varP, warmstart) - return λ, Ps + varλ = m[:λ] + varP = data.Ps + λ, Ps = PropertyT.λandP(data.name, m, varλ, varP, warmstart) + return λ, Ps end function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) - info(LOGGER, "Solving SDP problem...") - λ, Ps = λandP(m, data, sett.warmstart) + info(LOGGER, "Solving SDP problem...") + λ, Ps = λandP(m, data, sett.warmstart) - info(LOGGER, "Reconstructing P...") + info(LOGGER, "Reconstructing P...") - preps = load_preps(joinpath(prepath(sett), "preps.jld"), sett.autS) + preps = load_preps(joinpath(prepath(sett), "preps.jld"), sett.autS) - @logtime LOGGER recP = reconstruct_sol(preps, data.Us, Ps, data.dims) + @logtime LOGGER recP = reconstruct_sol(preps, data.Us, Ps, data.dims) - fname = PropertyT.λSDPfilenames(fullpath(sett))[2] - save(fname, "origP", Ps, "P", recP) - return λ, recP + fname = PropertyT.λSDPfilenames(fullpath(sett))[2] + save(fname, "origP", Ps, "P", recP) + return λ, recP end function load_preps(fname::String, G::Nemo.Group) @@ -229,49 +229,49 @@ end function check_property_T(sett::Settings) - init_orbit_data(LOGGER, sett, radius=sett.radius) + init_orbit_data(LOGGER, sett, radius=sett.radius) - if !sett.warmstart && all(isfile.(λSDPfilenames(fullpath(sett)))) - λ, 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) + if !sett.warmstart && all(isfile.(λSDPfilenames(fullpath(sett)))) + λ, 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) - end + λ, 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))") + info(LOGGER, "λ = $λ") + info(LOGGER, "sum(P) = $(sum(P))") + info(LOGGER, "maximum(P) = $(maximum(P))") + info(LOGGER, "minimum(P) = $(minimum(P))") - if λ > 0 - pm_fname, Δ_fname = pmΔfilenames(prepath(sett)) - RG = GroupRing(sett.G, load(pm_fname, "pm")) - Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) + if λ > 0 + pm_fname, Δ_fname = pmΔfilenames(prepath(sett)) + RG = GroupRing(sett.G, load(pm_fname, "pm")) + Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) - isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || - warn("The solution matrix doesn't seem to be positive definite!") - # @assert P == Symmetric(P) - @logtime LOGGER Q = real(sqrtm(Symmetric(P))) + isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || + warn("The solution matrix doesn't seem to be positive definite!") + # @assert P == Symmetric(P) + @logtime LOGGER Q = real(sqrtm(Symmetric(P))) - sgap = distance_to_positive_cone(Δ, λ, Q, 2*sett.radius, LOGGER) - if isa(sgap, Interval) - sgap = sgap.lo - end - if sgap > 0 - info(LOGGER, "λ ≥ $(Float64(trunc(sgap,12)))") + sgap = distance_to_positive_cone(Δ, λ, Q, 2*sett.radius, LOGGER) + 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 + 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 diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index 1c5d357..93a122b 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -57,7 +57,7 @@ function orbit_decomposition(G::Nemo.Group, E::Vector, rdict=GroupRings.reverse_ orbit = zeros(Int, length(elts)) a = E[i] Threads.@threads for i in 1:length(elts) - orbit[i] = rdict[elts[i](a)] + orbit[i] = rdict[elts[i](a)] end tovisit[orbit] = false push!(orbits, unique(orbit)) @@ -110,110 +110,110 @@ function matrix_reps{T<:GroupElem}(preps::Dict{T,perm}) end function perm_repr(g::GroupElem, E::Vector, E_dict) - p = Vector{Int}(length(E)) - for (i,elt) in enumerate(E) - p[i] = E_dict[g(elt)] - end - return p + p = Vector{Int}(length(E)) + for (i,elt) in enumerate(E) + p[i] = E_dict[g(elt)] + end + return p end function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E)) - elts = collect(elements(G)) - l = length(elts) - preps = Vector{Generic.perm}(l) + elts = collect(elements(G)) + l = length(elts) + preps = Vector{Generic.perm}(l) - permG = Nemo.PermutationGroup(length(E)) + permG = Nemo.PermutationGroup(length(E)) - Threads.@threads for i in 1:l - preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict)) - end + Threads.@threads for i in 1:l + preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict)) + end - return Dict(elts[i]=>preps[i] for i in 1:l) + return Dict(elts[i]=>preps[i] for i in 1:l) end function perm_reps(S::Vector, autS::Group, radius::Int) - E, _ = Groups.generate_balls(S, radius=radius) - return perm_reps(autS, E) + E, _ = Groups.generate_balls(S, radius=radius) + return perm_reps(autS, E) end 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) - transfP = [dims[π].*Us[π]*Ps[π]*Us[π]' for π in 1:l] - tmp = [zeros(Float64, size(first(transfP))) for _ in 1:l] - perms = collect(keys(preps)) + l = length(Us) + transfP = [dims[π].*Us[π]*Ps[π]*Us[π]' for π in 1:l] + tmp = [zeros(Float64, size(first(transfP))) for _ in 1:l] + perms = collect(keys(preps)) - @inbounds Threads.@threads for π in 1:l - for p in perms - BLAS.axpy!(1.0, view(transfP[π], preps[p].d, preps[p].d), tmp[π]) - end - end + @inbounds Threads.@threads for π in 1:l + for p in perms + BLAS.axpy!(1.0, view(transfP[π], preps[p].d, preps[p].d), tmp[π]) + end + end - recP = 1/length(perms) .* sum(tmp) - recP[abs.(recP) .< eps(eltype(recP))] = zero(eltype(recP)) - return recP + recP = 1/length(perms) .* sum(tmp) + recP[abs.(recP) .< eps(eltype(recP))] = zero(eltype(recP)) + return recP end 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 function orthSVD{T}(M::AbstractMatrix{T}) - M = full(M) - fact = svdfact(M) - M_rank = sum(fact[:S] .> maximum(size(M))*eps(T)) - return fact[:U][:,1:M_rank] + M = full(M) + fact = svdfact(M) + M_rank = sum(fact[:S] .> maximum(size(M))*eps(T)) + return fact[:U][:,1:M_rank] end function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, 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? - Id = (isa(G, Nemo.Ring) ? one(G) : G()) + # TODO: Fix that by multiple dispatch? + Id = (isa(G, Nemo.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) + @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) - info(logger, "Product matrix") - @logtime logger 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.augmentation(Δ) == 0 - save(joinpath(name, "delta.jld"), "Δ", Δ.coeffs) - save(joinpath(name, "pm.jld"), "pm", pm) + info(logger, "Product matrix") + @logtime logger 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.augmentation(Δ) == 0 + save(joinpath(name, "delta.jld"), "Δ", Δ.coeffs) + save(joinpath(name, "pm.jld"), "pm", pm) - info(logger, "Decomposing E into orbits of $(autS)") - @logtime logger 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!") - save(joinpath(name, "orbits.jld"), "orbits", orbs) + info(logger, "Decomposing E into orbits of $(autS)") + @logtime logger 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!") + save(joinpath(name, "orbits.jld"), "orbits", orbs) - info(logger, "Action matrices") - @logtime logger reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict) - save_preps(joinpath(name, "preps.jld"), reps) - reps = matrix_reps(reps) + info(logger, "Action matrices") + @logtime logger reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict) + save_preps(joinpath(name, "preps.jld"), reps) + reps = matrix_reps(reps) - info(logger, "Projections") - @logtime logger autS_mps = rankOne_projections(autS); + info(logger, "Projections") + @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...") - @logtime logger Uπs = orthSVD.(π_E_projections) + info(logger, "Uπs...") + @logtime logger Uπs = orthSVD.(π_E_projections) - multiplicities = size.(Uπs,2) - info(logger, "multiplicities = $multiplicities") - dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps]; - info(logger, "dimensions = $dimensions") - @assert dot(multiplicities, dimensions) == sizes[radius] + multiplicities = size.(Uπs,2) + info(logger, "multiplicities = $multiplicities") + dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps]; + info(logger, "dimensions = $dimensions") + @assert dot(multiplicities, dimensions) == sizes[radius] - save(joinpath(name, "U_pis.jld"), - "Uπs", Uπs, - "dims", dimensions) - return 0 + save(joinpath(name, "U_pis.jld"), + "Uπs", Uπs, + "dims", dimensions) + return 0 end diff --git a/src/Projections.jl b/src/Projections.jl index ae468d1..4506846 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -7,17 +7,17 @@ abstract type AbstractCharacter end struct PermCharacter <: AbstractCharacter - p::Generic.Partition + p::Generic.Partition end struct DirectProdCharacter <: AbstractCharacter - i::Int + i::Int end function (chi::PermCharacter)(g::Generic.perm) - R = Nemo.partitionseq(chi.p) - p = Partition(Nemo.Generic.permtype(g)) - return Int(Nemo.Generic.MN1inner(R, p, 1, Nemo.Generic._charvalsTable)) + R = Nemo.partitionseq(chi.p) + p = Partition(Nemo.Generic.permtype(g)) + return Int(Nemo.Generic.MN1inner(R, p, 1, Nemo.Generic._charvalsTable)) end Nemo.isone(p::GroupElem) = p == parent(p)() @@ -29,23 +29,23 @@ end ## NOTE: this works only for Z/2!!!! 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 for T in [PermCharacter, DirectProdCharacter] - @eval begin - function (chi::$T)(X::GroupRingElem) - RG = parent(X) - z = zero(eltype(X)) - result = z - for i in 1:length(X.coeffs) - if X.coeffs[i] != z - result += chi(RG.basis[i])*X.coeffs[i] + @eval begin + function (chi::$T)(X::GroupRingElem) + RG = parent(X) + z = zero(eltype(X)) + result = z + for i in 1:length(X.coeffs) + if X.coeffs[i] != z + result += chi(RG.basis[i])*X.coeffs[i] + end end - end - return result - end - end + return result + end + end end ############################################################################### @@ -55,16 +55,16 @@ end ############################################################################### function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int}) - result = RG(T) - result.coeffs = full(result.coeffs) - dim = chi(RG.group()) - ord = Int(order(RG.group)) + result = RG(T) + result.coeffs = full(result.coeffs) + dim = chi(RG.group()) + ord = Int(order(RG.group)) - for g in RG.basis - result[g] = convert(T, (dim//ord)*chi(g)) - end + for g in RG.basis + result[g] = convert(T, (dim//ord)*chi(g)) + end - return result + return result end 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 function rankOne_projection(chi::PropertyT.PermCharacter, - idems::Vector{T}) where {T<:GroupRingElem} + idems::Vector{T}) where {T<:GroupRingElem} - RG = parent(first(idems)) - S = eltype(first(idems)) + RG = parent(first(idems)) + S = eltype(first(idems)) - ids = [one(RG, S); idems] - zzz = zero(S) + ids = [one(RG, S); idems] + zzz = zero(S) - for (i,j,k) in Base.product(ids, ids, ids) - if chi(i) == zzz || chi(j) == zzz || chi(k) == zzz - continue - end - elt = i*j*k - elt^2 == elt || continue - if chi(elt) == one(S) - return elt - # return (i,j,k) - end - end - throw("Couldn't find rank-one projection for $chi") + for (i,j,k) in Base.product(ids, ids, ids) + if chi(i) == zzz || chi(j) == zzz || chi(k) == zzz + continue + end + elt = i*j*k + elt^2 == elt || continue + if chi(elt) == one(S) + return elt + # return (i,j,k) + end + end + throw("Couldn't find rank-one projection for $chi") end function rankOne_projections(G::Generic.PermGroup, T::Type=Rational{Int}) - if G.n == 1 - return [one(GroupRing(G), T)] - elseif G.n < 8 - RG = GroupRing(G, fastm=true) - else - RG = GroupRing(G, fastm=false) - end + if G.n == 1 + return [one(GroupRing(G), T)] + elseif G.n < 8 + RG = GroupRing(G, fastm=true) + else + RG = GroupRing(G, fastm=false) + end - RGidems = idempotents(RG, T) - l = length(Partitions(G.n)) + RGidems = idempotents(RG, T) + l = length(Partitions(G.n)) - parts = collect(Partitions(G.n)) - chars = [PropertyT.PermCharacter(p) for p in parts] - min_projs = Vector{eltype(RGidems)}(l) + parts = collect(Partitions(G.n)) + chars = [PropertyT.PermCharacter(p) for p in parts] + min_projs = Vector{eltype(RGidems)}(l) - for i in 1:l - chi = PropertyT.PermCharacter(parts[i]) - min_projs[i] = rankOne_projection(chi,RGidems)*central_projection(RG,chi) - end + for i in 1:l + chi = PropertyT.PermCharacter(parts[i]) + min_projs[i] = rankOne_projection(chi,RGidems)*central_projection(RG,chi) + end - return min_projs + return min_projs end function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int}) - N = BN.P.n - # projections as elements of the group rings RSₙ - SNprojs_nc = [rankOne_projections(PermutationGroup(i)) for i in 1:N] + N = BN.P.n + # projections as elements of the group rings RSₙ + SNprojs_nc = [rankOne_projections(PermutationGroup(i)) for i in 1:N] - # embedding into group ring of BN - RBN = GroupRing(BN) - RFFFF_projs = [central_projection(GroupRing(BN.N), DirectProdCharacter(i),T) - for i in 1:BN.P.n] + # embedding into group ring of BN + RBN = GroupRing(BN) + RFFFF_projs = [central_projection(GroupRing(BN.N), DirectProdCharacter(i),T) + for i in 1:BN.P.n] - e0 = central_projection(GroupRing(BN.N), DirectProdCharacter(0), T) - Q0 = RBN(e0, g -> BN(g)) - Qs = [RBN(q, g -> BN(g)) for q in RFFFF_projs] + e0 = central_projection(GroupRing(BN.N), DirectProdCharacter(0), T) + Q0 = RBN(e0, g -> BN(g)) + 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) - for i in 1:N-1 - 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])) + range = collect(1:N) + for i in 1:N-1 + 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])) - 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_first = [RBN(p, first_emb) for p in SNprojs_nc[i]] + Sk_last = [RBN(p, last_emb) for p in SNprojs_nc[N-i]] - append!(all_projs, - [Qs[i]*p1*p2 for (p1,p2) in Base.product(Sk_first,Sk_last)]) - end + append!(all_projs, + [Qs[i]*p1*p2 for (p1,p2) in Base.product(Sk_first,Sk_last)]) + 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 -end + return all_projs + end -############################################################################## -# -# General Groups Misc -# -############################################################################## + ############################################################################## + # + # General Groups Misc + # + ############################################################################## -doc""" + doc""" products(X::Vector{GroupElem}, Y::Vector{GroupElem}, op=*) -> 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 -> forming 'products' by adding `op` (which is `*` by default). -""" -function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) - result = Vector{T}() - seen = Set{T}() - for x in X - for y in Y - z = op(x,y) - if !in(z, seen) - push!(seen, z) - push!(result, z) - end - end - end - return result -end + > 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 + > forming 'products' by adding `op` (which is `*` by default). + """ + function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) + result = Vector{T}() + seen = Set{T}() + for x in X + for y in Y + z = op(x,y) + if !in(z, seen) + push!(seen, z) + push!(result, z) + end + end + end + return result + end -doc""" + doc""" 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 -> radius `r` (word-length metric induced by `gens`). -> If `r(=2)` is specified the procedure will terminate after generating ball -> of radius `r` in the word-length metric induced by `gens`. -> The identity element `Id` and binary operation function `op` can be supplied -> to e.g. take advantage of additive group structure. -""" -function generateGroup{T<:GroupElem}(gens::Vector{T}, r=2, Id::T=parent(first(gens))(), op=*) - n = 0 - R = 1 - elts = gens - gens = [Id; gens] - while n ≠ length(elts) && R < r - # @show elts - R += 1 - n = length(elts) - elts = products(elts, gens, op) - end - return elts -end + > Produces all elements of a group generated by elements in `gens` in ball of + > radius `r` (word-length metric induced by `gens`). + > If `r(=2)` is specified the procedure will terminate after generating ball + > of radius `r` in the word-length metric induced by `gens`. + > The identity element `Id` and binary operation function `op` can be supplied + > to e.g. take advantage of additive group structure. + """ + function generateGroup{T<:GroupElem}(gens::Vector{T}, r=2, Id::T=parent(first(gens))(), op=*) + n = 0 + R = 1 + elts = gens + gens = [Id; gens] + while n ≠ length(elts) && R < r + # @show elts + R += 1 + n = length(elts) + elts = products(elts, gens, op) + end + return elts + end diff --git a/src/PropertyT.jl b/src/PropertyT.jl index a73cc45..6239b53 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -68,22 +68,20 @@ function time_string(elapsedtime, bytes, gctime, allocs) return str end -function exists(fname::String) - return isfile(fname) || islink(fname) -end +exists(fname::String) = isfile(fname) || islink(fname) function pmΔfilenames(prefix::String) - isdir(prefix) || mkdir(prefix) - pm_filename = joinpath(prefix, "pm.jld") - Δ_coeff_filename = joinpath(prefix, "delta.jld") - return pm_filename, Δ_coeff_filename + isdir(prefix) || mkdir(prefix) + pm_filename = joinpath(prefix, "pm.jld") + Δ_coeff_filename = joinpath(prefix, "delta.jld") + return pm_filename, Δ_coeff_filename end function λSDPfilenames(prefix::String) - isdir(prefix) || mkdir(prefix) - λ_filename = joinpath(prefix, "lambda.jld") - SDP_filename = joinpath(prefix, "SDPmatrix.jld") - return λ_filename, SDP_filename + isdir(prefix) || mkdir(prefix) + λ_filename = joinpath(prefix, "lambda.jld") + SDP_filename = joinpath(prefix, "SDPmatrix.jld") + return λ_filename, SDP_filename end function ΔandSDPconstraints(prefix::String, G::Group) @@ -100,12 +98,12 @@ function ΔandSDPconstraints(prefix::String, G::Group) 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 + 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) @@ -148,27 +146,26 @@ function λandP(name::String, SDP_problem::JuMP.Model, varλ, varP, warmstart=fa handler.levels.x = LOGGER_SOLVER.levels LOGGER_SOLVER.handlers["solver_log"] = handler - if warmstart && isfile(joinpath(name, "warmstart.jld")) - ws = load(joinpath(name, "warmstart.jld"), "warmstart") - else - ws = nothing - end + 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) + λ, P, warmstart = compute_λandP(SDP_problem, varλ, varP, warmstart=ws) - delete!(LOGGER_SOLVER.handlers, "solver_log") + delete!(LOGGER_SOLVER.handlers, "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 + λ_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) @@ -311,50 +308,48 @@ function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius) Δ, sdp_constraints = ΔandSDPconstraints(name, S, Id, radius=radius) end - if all(exists.(λSDPfilenames(name))) - λ, P = λandP(name) - else - info(LOGGER, "Creating SDP problem...") - SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound) - JuMP.setsolver(SDP_problem, solver) + if all(exists.(λSDPfilenames(name))) + λ, P = λandP(name) + else + info(LOGGER, "Creating SDP problem...") + SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound) + JuMP.setsolver(SDP_problem, solver) + λ, P = λandP(name, SDP_problem, λ, P) + end - λ, P = λandP(name, SDP_problem, λ, P) - end + info(LOGGER, "λ = $λ") + info(LOGGER, "sum(P) = $(sum(P))") + info(LOGGER, "maximum(P) = $(maximum(P))") + info(LOGGER, "minimum(P) = $(minimum(P))") - info(LOGGER, "λ = $λ") - info(LOGGER, "sum(P) = $(sum(P))") - info(LOGGER, "maximum(P) = $(maximum(P))") - info(LOGGER, "minimum(P) = $(minimum(P))") + if λ > 0 + pm_fname, Δ_fname = pmΔfilenames(name) + RG = GroupRing(parent(first(S)), load(pm_fname, "pm")) + Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) - if λ > 0 - pm_fname, Δ_fname = pmΔfilenames(name) - RG = GroupRing(parent(first(S)), load(pm_fname, "pm")) - Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) + isapprox(eigvals(P), abs(eigvals(P)), atol=tol) || + warn("The solution matrix doesn't seem to be positive definite!") + @logtime LOGGER Q = real(sqrtm(Symmetric(P))) - isapprox(eigvals(P), abs(eigvals(P)), atol=tol) || - warn("The solution matrix doesn't seem to be positive definite!") - # @assert P == Symmetric(P) - @logtime LOGGER Q = real(sqrtm(Symmetric(P))) - - sgap = distance_to_positive_cone(Δ, λ, Q, 2*radius, LOGGER) - 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 + sgap = distance_to_positive_cone(Δ, λ, Q, 2*radius, LOGGER) + 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 include("SDPs.jl")