fix indentation to 4 spaces
This commit is contained in:
parent
fc54803b58
commit
47f6d3637e
|
|
@ -23,80 +23,80 @@ 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("")
|
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 distance_to_cone{T<:Interval}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int)
|
||||||
SOS = compute_SOS(Q, parent(elt), length(elt.coeffs))
|
SOS = compute_SOS(Q, parent(elt), length(elt.coeffs))
|
||||||
SOS_diff = elt - SOS
|
SOS_diff = elt - SOS
|
||||||
|
|
||||||
ɛ_dist = GroupRings.augmentation(SOS_diff)
|
ɛ_dist = GroupRings.augmentation(SOS_diff)
|
||||||
info(LOGGER, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
|
info(LOGGER, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
|
||||||
|
|
||||||
eoi_SOS_L1_dist = norm(SOS_diff,1)
|
eoi_SOS_L1_dist = norm(SOS_diff,1)
|
||||||
info(LOGGER, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)")
|
info(LOGGER, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)")
|
||||||
|
|
||||||
dist = 2^(wlen-1)*eoi_SOS_L1_dist
|
dist = 2^(wlen-1)*eoi_SOS_L1_dist
|
||||||
return dist
|
return dist
|
||||||
end
|
end
|
||||||
|
|
||||||
function distance_to_cone{T}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int)
|
function distance_to_cone{T}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int)
|
||||||
SOS = compute_SOS(Q, parent(elt), length(elt.coeffs))
|
SOS = compute_SOS(Q, parent(elt), length(elt.coeffs))
|
||||||
SOS_diff = elt - SOS
|
SOS_diff = elt - SOS
|
||||||
|
|
||||||
ɛ_dist = GroupRings.augmentation(SOS_diff)
|
ɛ_dist = GroupRings.augmentation(SOS_diff)
|
||||||
info(LOGGER, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
|
info(LOGGER, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
|
||||||
|
|
||||||
eoi_SOS_L1_dist = norm(SOS_diff,1)
|
eoi_SOS_L1_dist = norm(SOS_diff,1)
|
||||||
info(LOGGER, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L1_dist))")
|
info(LOGGER, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L1_dist))")
|
||||||
|
|
||||||
dist = 2^(wlen-1)*eoi_SOS_L1_dist
|
dist = 2^(wlen-1)*eoi_SOS_L1_dist
|
||||||
return 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, logger)
|
function distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen::Int, logger)
|
||||||
|
|
|
||||||
|
|
@ -4,16 +4,16 @@ 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
|
||||||
end
|
end
|
||||||
|
|
||||||
prefix(s::Settings) = s.name
|
prefix(s::Settings) = s.name
|
||||||
|
|
@ -22,43 +22,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(joinpath(prepath(sett), "delta.jld"), "Δ");
|
||||||
pm = load(joinpath(prepath(sett), "pm.jld"), "pm");
|
pm = load(joinpath(prepath(sett), "pm.jld"), "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(joinpath(prepath(sett), "U_pis.jld"), "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(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");
|
orbits = load(joinpath(prepath(sett), "orbits.jld"), "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 +67,152 @@ 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(LOGGER, "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(LOGGER, "Sparsification decreased the rank!")
|
||||||
end
|
end
|
||||||
|
|
||||||
if verbose
|
if verbose
|
||||||
info(LOGGER, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M),20))
|
info(LOGGER, "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)
|
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)
|
if !all(files_exists)
|
||||||
compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.autS, radius=radius)
|
compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.autS, radius=radius)
|
||||||
end
|
end
|
||||||
|
|
||||||
return 0
|
return 0
|
||||||
end
|
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(LOGGER, "Loading orbit data....")
|
||||||
@logtime LOGGER SDP_problem, orb_data = OrbitData(sett);
|
@logtime 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(LOGGER, "Adding $t constraints ... ")
|
||||||
@logtime LOGGER addconstraints!(SDP_problem, orb_data)
|
@logtime 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(LOGGER, "Solving SDP problem...")
|
||||||
λ, Ps = λandP(m, data, sett.warmstart)
|
λ, 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]
|
fname = PropertyT.λSDPfilenames(fullpath(sett))[2]
|
||||||
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 +229,49 @@ end
|
||||||
|
|
||||||
function check_property_T(sett::Settings)
|
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))))
|
if !sett.warmstart && all(isfile.(λSDPfilenames(fullpath(sett))))
|
||||||
λ, P = PropertyT.λandP(fullpath(sett))
|
λ, P = PropertyT.λandP(fullpath(sett))
|
||||||
else
|
else
|
||||||
info(LOGGER, "Creating SDP problem...")
|
info(LOGGER, "Creating SDP problem...")
|
||||||
SDP_problem, orb_data = create_SDP_problem(sett)
|
SDP_problem, orb_data = create_SDP_problem(sett)
|
||||||
JuMP.setsolver(SDP_problem, sett.solver)
|
JuMP.setsolver(SDP_problem, sett.solver)
|
||||||
|
|
||||||
λ, P = λandP(SDP_problem, orb_data, sett)
|
λ, P = λandP(SDP_problem, orb_data, sett)
|
||||||
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
|
if λ > 0
|
||||||
pm_fname, Δ_fname = pmΔfilenames(prepath(sett))
|
pm_fname, Δ_fname = pmΔfilenames(prepath(sett))
|
||||||
RG = GroupRing(sett.G, load(pm_fname, "pm"))
|
RG = GroupRing(sett.G, load(pm_fname, "pm"))
|
||||||
Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG)
|
Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG)
|
||||||
|
|
||||||
isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) ||
|
isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) ||
|
||||||
warn("The solution matrix doesn't seem to be positive definite!")
|
warn("The solution matrix doesn't seem to be positive definite!")
|
||||||
# @assert P == Symmetric(P)
|
# @assert P == Symmetric(P)
|
||||||
@logtime LOGGER Q = real(sqrtm(Symmetric(P)))
|
@logtime LOGGER Q = real(sqrtm(Symmetric(P)))
|
||||||
|
|
||||||
sgap = distance_to_positive_cone(Δ, λ, Q, 2*sett.radius, LOGGER)
|
sgap = distance_to_positive_cone(Δ, λ, Q, 2*sett.radius, LOGGER)
|
||||||
if isa(sgap, Interval)
|
if isa(sgap, Interval)
|
||||||
sgap = sgap.lo
|
sgap = sgap.lo
|
||||||
end
|
end
|
||||||
if sgap > 0
|
if sgap > 0
|
||||||
info(LOGGER, "λ ≥ $(Float64(trunc(sgap,12)))")
|
info(LOGGER, "λ ≥ $(Float64(trunc(sgap,12)))")
|
||||||
Kazhdan_κ = PropertyT.Kazhdan_from_sgap(sgap, length(sett.S))
|
Kazhdan_κ = PropertyT.Kazhdan_from_sgap(sgap, length(sett.S))
|
||||||
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
|
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
|
||||||
info(LOGGER, "κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
info(LOGGER, "κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
||||||
return true
|
return true
|
||||||
else
|
else
|
||||||
sgap = Float64(trunc(sgap, 12))
|
sgap = Float64(trunc(sgap, 12))
|
||||||
info(LOGGER, "λ($(sett.name), S) ≥ $sgap: Group may NOT HAVE property (T)!")
|
info(LOGGER, "λ($(sett.name), S) ≥ $sgap: Group may NOT HAVE property (T)!")
|
||||||
return false
|
return false
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
info(LOGGER, "κ($(sett.name), S) ≥ $λ < 0: Tells us nothing about property (T)")
|
info(LOGGER, "κ($(sett.name), S) ≥ $λ < 0: Tells us nothing about property (T)")
|
||||||
return false
|
return false
|
||||||
end
|
end
|
||||||
|
|
|
||||||
|
|
@ -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,110 @@ 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, 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?
|
# TODO: Fix that by multiple dispatch?
|
||||||
Id = (isa(G, Nemo.Ring) ? one(G) : G())
|
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(joinpath(name, "delta.jld"), "Δ", Δ.coeffs)
|
||||||
save(joinpath(name, "pm.jld"), "pm", pm)
|
save(joinpath(name, "pm.jld"), "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(joinpath(name, "preps.jld"), 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
|
||||||
|
|
|
||||||
|
|
@ -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
|
||||||
|
|
|
||||||
141
src/PropertyT.jl
141
src/PropertyT.jl
|
|
@ -68,22 +68,20 @@ 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)
|
|
||||||
end
|
|
||||||
|
|
||||||
function pmΔfilenames(prefix::String)
|
function pmΔfilenames(prefix::String)
|
||||||
isdir(prefix) || mkdir(prefix)
|
isdir(prefix) || mkdir(prefix)
|
||||||
pm_filename = joinpath(prefix, "pm.jld")
|
pm_filename = joinpath(prefix, "pm.jld")
|
||||||
Δ_coeff_filename = joinpath(prefix, "delta.jld")
|
Δ_coeff_filename = joinpath(prefix, "delta.jld")
|
||||||
return pm_filename, Δ_coeff_filename
|
return pm_filename, Δ_coeff_filename
|
||||||
end
|
end
|
||||||
|
|
||||||
function λSDPfilenames(prefix::String)
|
function λSDPfilenames(prefix::String)
|
||||||
isdir(prefix) || mkdir(prefix)
|
isdir(prefix) || mkdir(prefix)
|
||||||
λ_filename = joinpath(prefix, "lambda.jld")
|
λ_filename = joinpath(prefix, "lambda.jld")
|
||||||
SDP_filename = joinpath(prefix, "SDPmatrix.jld")
|
SDP_filename = joinpath(prefix, "SDPmatrix.jld")
|
||||||
return λ_filename, SDP_filename
|
return λ_filename, SDP_filename
|
||||||
end
|
end
|
||||||
|
|
||||||
function ΔandSDPconstraints(prefix::String, G::Group)
|
function ΔandSDPconstraints(prefix::String, G::Group)
|
||||||
|
|
@ -100,12 +98,12 @@ function ΔandSDPconstraints(prefix::String, G::Group)
|
||||||
end
|
end
|
||||||
|
|
||||||
function ΔandSDPconstraints{T<:GroupElem}(name::String, S::Vector{T}, Id::T; radius::Int=2)
|
function ΔandSDPconstraints{T<:GroupElem}(name::String, S::Vector{T}, Id::T; radius::Int=2)
|
||||||
info(LOGGER, "Computing pm, Δ, sdp_constraints...")
|
info(LOGGER, "Computing pm, Δ, sdp_constraints...")
|
||||||
Δ, sdp_constraints = ΔandSDPconstraints(S, Id, radius=radius)
|
Δ, sdp_constraints = ΔandSDPconstraints(S, Id, radius=radius)
|
||||||
pm_fname, Δ_fname = pmΔfilenames(name)
|
pm_fname, Δ_fname = pmΔfilenames(name)
|
||||||
save(pm_fname, "pm", parent(Δ).pm)
|
save(pm_fname, "pm", parent(Δ).pm)
|
||||||
save(Δ_fname, "Δ", Δ.coeffs)
|
save(Δ_fname, "Δ", Δ.coeffs)
|
||||||
return Δ, sdp_constraints
|
return Δ, sdp_constraints
|
||||||
end
|
end
|
||||||
|
|
||||||
function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2)
|
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
|
handler.levels.x = LOGGER_SOLVER.levels
|
||||||
LOGGER_SOLVER.handlers["solver_log"] = handler
|
LOGGER_SOLVER.handlers["solver_log"] = handler
|
||||||
|
|
||||||
if warmstart && isfile(joinpath(name, "warmstart.jld"))
|
if warmstart && isfile(joinpath(name, "warmstart.jld"))
|
||||||
ws = load(joinpath(name, "warmstart.jld"), "warmstart")
|
ws = load(joinpath(name, "warmstart.jld"), "warmstart")
|
||||||
else
|
else
|
||||||
ws = nothing
|
ws = nothing
|
||||||
end
|
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)
|
λ_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
|
|
||||||
|
|
||||||
|
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
|
end
|
||||||
|
|
||||||
function fillfrominternal!(m::JuMP.Model, traits)
|
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)
|
Δ, sdp_constraints = ΔandSDPconstraints(name, S, Id, radius=radius)
|
||||||
end
|
end
|
||||||
|
|
||||||
if all(exists.(λSDPfilenames(name)))
|
if all(exists.(λSDPfilenames(name)))
|
||||||
λ, P = λandP(name)
|
λ, P = λandP(name)
|
||||||
else
|
else
|
||||||
info(LOGGER, "Creating SDP problem...")
|
info(LOGGER, "Creating SDP problem...")
|
||||||
SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound)
|
SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound)
|
||||||
JuMP.setsolver(SDP_problem, solver)
|
JuMP.setsolver(SDP_problem, solver)
|
||||||
|
|
||||||
|
λ, P = λandP(name, SDP_problem, λ, P)
|
||||||
|
end
|
||||||
|
|
||||||
λ, P = λandP(name, SDP_problem, λ, P)
|
info(LOGGER, "λ = $λ")
|
||||||
end
|
info(LOGGER, "sum(P) = $(sum(P))")
|
||||||
|
info(LOGGER, "maximum(P) = $(maximum(P))")
|
||||||
|
info(LOGGER, "minimum(P) = $(minimum(P))")
|
||||||
|
|
||||||
info(LOGGER, "λ = $λ")
|
if λ > 0
|
||||||
info(LOGGER, "sum(P) = $(sum(P))")
|
pm_fname, Δ_fname = pmΔfilenames(name)
|
||||||
info(LOGGER, "maximum(P) = $(maximum(P))")
|
RG = GroupRing(parent(first(S)), load(pm_fname, "pm"))
|
||||||
info(LOGGER, "minimum(P) = $(minimum(P))")
|
Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG)
|
||||||
|
|
||||||
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"))
|
@logtime LOGGER Q = real(sqrtm(Symmetric(P)))
|
||||||
Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG)
|
|
||||||
|
|
||||||
isapprox(eigvals(P), abs(eigvals(P)), atol=tol) ||
|
sgap = distance_to_positive_cone(Δ, λ, Q, 2*radius, LOGGER)
|
||||||
warn("The solution matrix doesn't seem to be positive definite!")
|
if isa(sgap, Interval)
|
||||||
# @assert P == Symmetric(P)
|
sgap = sgap.lo
|
||||||
@logtime LOGGER Q = real(sqrtm(Symmetric(P)))
|
end
|
||||||
|
if sgap > 0
|
||||||
sgap = distance_to_positive_cone(Δ, λ, Q, 2*radius, LOGGER)
|
info(LOGGER, "λ ≥ $(Float64(trunc(sgap,12)))")
|
||||||
if isa(sgap, Interval)
|
Kazhdan_κ = Kazhdan_from_sgap(sgap, length(S))
|
||||||
sgap = sgap.lo
|
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
|
||||||
end
|
info(LOGGER, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
||||||
if sgap > 0
|
return true
|
||||||
info(LOGGER, "λ ≥ $(Float64(trunc(sgap,12)))")
|
else
|
||||||
Kazhdan_κ = Kazhdan_from_sgap(sgap, length(S))
|
sgap = Float64(trunc(sgap, 12))
|
||||||
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
|
info(LOGGER, "λ($name, S) ≥ $sgap: Group may NOT HAVE property (T)!")
|
||||||
info(LOGGER, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
return false
|
||||||
return true
|
end
|
||||||
else
|
end
|
||||||
sgap = Float64(trunc(sgap, 12))
|
info(LOGGER, "κ($name, S) ≥ $λ < 0: Tells us nothing about property (T)")
|
||||||
info(LOGGER, "λ($name, S) ≥ $sgap: Group may NOT HAVE property (T)!")
|
return false
|
||||||
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")
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue