2017-06-22 14:12:35 +02:00
|
|
|
using JuMP
|
|
|
|
using SCS
|
|
|
|
|
2017-06-22 15:15:43 +02:00
|
|
|
export Settings, OrbitData
|
2017-06-22 14:12:35 +02:00
|
|
|
|
2017-10-09 19:45:28 +02:00
|
|
|
|
2018-08-20 03:54:03 +02:00
|
|
|
struct OrbitData{T<:AbstractArray{Float64, 2}}
|
|
|
|
orbits::Vector{Vector{Int}}
|
|
|
|
Uπs::Vector{T}
|
2018-01-01 14:06:33 +01:00
|
|
|
dims::Vector{Int}
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2017-10-09 19:46:08 +02:00
|
|
|
function OrbitData(sett::Settings)
|
2018-08-20 03:54:03 +02:00
|
|
|
info("Loading Uπs, dims, orbits...")
|
2018-01-02 02:55:53 +01:00
|
|
|
Uπs = load(filename(prepath(sett), :Uπs), "Uπs")
|
2018-01-01 14:06:33 +01:00
|
|
|
nzros = [i for i in 1:length(Uπs) if size(Uπs[i],2) !=0]
|
|
|
|
Uπs = Uπs[nzros]
|
2017-11-13 10:57:57 +01:00
|
|
|
|
2018-08-20 03:54:03 +02:00
|
|
|
Uπs = map(x -> sparsify!(x, sett.tol/100, verbose=true), Uπs)
|
2018-01-01 14:06:33 +01:00
|
|
|
#dimensions of the corresponding πs:
|
2018-01-02 02:55:53 +01:00
|
|
|
dims = load(filename(prepath(sett), :Uπs), "dims")[nzros]
|
2017-06-22 15:09:46 +02:00
|
|
|
|
2018-08-20 03:54:03 +02:00
|
|
|
orbits = load(filename(prepath(sett), :orbits), "orbits");
|
2017-06-22 15:09:46 +02:00
|
|
|
|
2018-08-20 03:54:03 +02:00
|
|
|
return OrbitData(orbits, Uπs, dims)
|
2017-06-22 15:09:46 +02:00
|
|
|
end
|
|
|
|
|
2017-06-22 14:12:35 +02:00
|
|
|
include("OrbitDecomposition.jl")
|
|
|
|
|
2018-07-31 10:18:03 +02:00
|
|
|
dens(M::SparseMatrixCSC) = nnz(M)/length(M)
|
|
|
|
dens(M::AbstractArray) = countnz(M)/length(M)
|
2017-07-26 10:29:11 +02:00
|
|
|
|
2017-08-27 21:23:42 +02:00
|
|
|
function sparsify!{Tv,Ti}(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv); verbose=false)
|
2018-01-01 14:06:33 +01:00
|
|
|
|
|
|
|
densM = dens(M)
|
|
|
|
for i in eachindex(M.nzval)
|
|
|
|
if abs(M.nzval[i]) < eps
|
|
|
|
M.nzval[i] = zero(Tv)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
dropzeros!(M)
|
|
|
|
|
|
|
|
if verbose
|
2018-07-31 10:18:03 +02:00
|
|
|
info("Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M), 20), " ($(nnz(M)) non-zeros)")
|
2018-01-01 14:06:33 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
return M
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2018-07-31 10:18:03 +02:00
|
|
|
function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); verbose=false)
|
2018-01-01 14:06:33 +01:00
|
|
|
densM = dens(M)
|
2018-07-31 10:18:03 +02:00
|
|
|
if verbose
|
|
|
|
info("Sparsifying $(size(M))-matrix... ")
|
|
|
|
end
|
2017-08-27 18:58:01 +02:00
|
|
|
|
2018-07-31 10:18:03 +02:00
|
|
|
for n in eachindex(M)
|
|
|
|
if abs(M[n]) < eps
|
|
|
|
M[n] = zero(T)
|
|
|
|
end
|
2018-01-01 14:06:33 +01:00
|
|
|
end
|
2017-08-27 18:58:01 +02:00
|
|
|
|
2018-01-01 14:06:33 +01:00
|
|
|
if verbose
|
2018-07-31 10:18:03 +02:00
|
|
|
info("$(rpad(densM, 20)) → $(rpad(dens(M),20))), ($(countnz(M)) non-zeros)")
|
2018-01-01 14:06:33 +01:00
|
|
|
end
|
2017-08-27 21:23:14 +02:00
|
|
|
|
2018-01-01 14:06:33 +01:00
|
|
|
return sparse(M)
|
2017-08-27 18:58:01 +02:00
|
|
|
end
|
|
|
|
|
2018-07-31 10:18:03 +02:00
|
|
|
sparsify{T}(U::AbstractArray{T}, tol=eps(T); verbose=false) = sparsify!(deepcopy(U), tol, verbose=verbose)
|
2017-08-27 18:58:01 +02:00
|
|
|
|
2017-08-27 18:59:16 +02:00
|
|
|
function constrLHS(m::JuMP.Model, cnstr, Us, Ust, dims, vars, eps=100*eps(1.0))
|
2018-01-01 14:06:33 +01:00
|
|
|
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)))
|
2017-07-26 18:02:21 +02:00
|
|
|
end
|
|
|
|
|
2018-08-20 04:00:14 +02:00
|
|
|
function addconstraints!(m::JuMP.Model, X::GroupRingElem, orderunit::GroupRingElem, λ::JuMP.Variable, P, data::OrbitData)
|
|
|
|
|
|
|
|
orderunit_orb = orbit_spvector(orderunit.coeffs, data.orbits)
|
|
|
|
X_orb = orbit_spvector(X.coeffs, data.orbits)
|
|
|
|
Ust = [U' for U in data.Uπs]
|
|
|
|
n = size(parent(X).pm, 1)
|
|
|
|
|
|
|
|
for t in 1:length(X_orb)
|
|
|
|
x, u = X_orb[t], orderunit_orb[t]
|
|
|
|
cnstrs = [constraint(parent(X).pm, o) for o in data.orbits[t]]
|
|
|
|
lhs = constrLHS(m, orbit_constraint(cnstrs,n), data.Uπs, Ust, data.dims, P)
|
|
|
|
|
|
|
|
JuMP.@constraint(m, lhs == x - λ*u)
|
2018-01-01 14:06:33 +01:00
|
|
|
end
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2018-08-20 03:59:39 +02:00
|
|
|
function init_model(m, sizes)
|
|
|
|
P = Vector{Array{JuMP.Variable,2}}(length(sizes))
|
2017-11-05 20:55:53 +01:00
|
|
|
|
2018-01-01 14:06:33 +01:00
|
|
|
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
|
2017-11-05 20:55:53 +01:00
|
|
|
|
2018-08-20 03:59:39 +02:00
|
|
|
return P
|
2017-06-22 15:11:14 +02:00
|
|
|
end
|
|
|
|
|
2018-08-20 03:59:39 +02:00
|
|
|
function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData; upper_bound=Inf)
|
|
|
|
m = JuMP.Model();
|
|
|
|
P = init_model(m, size.(data.Uπs,2))
|
2017-06-22 14:12:35 +02:00
|
|
|
|
2018-08-20 03:59:39 +02:00
|
|
|
λ = JuMP.@variable(m, λ)
|
|
|
|
if upper_bound < Inf
|
2018-09-03 10:23:08 +02:00
|
|
|
JuMP.@constraint(m, λ <= upper_bound)
|
2018-01-01 14:06:33 +01:00
|
|
|
end
|
2017-06-22 14:12:35 +02:00
|
|
|
|
2018-08-20 03:59:39 +02:00
|
|
|
info("Adding $(length(data.orbits)) constraints... ")
|
2017-06-22 14:12:35 +02:00
|
|
|
|
2018-08-20 03:59:39 +02:00
|
|
|
@time addconstraints!(m, X, orderunit, λ, P, data)
|
2017-06-22 14:12:35 +02:00
|
|
|
|
2018-08-20 03:59:39 +02:00
|
|
|
JuMP.@objective(m, Max, λ)
|
|
|
|
return m, λ, P
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
|
|
|
|
2017-08-27 19:06:17 +02:00
|
|
|
|
2018-08-20 04:00:58 +02:00
|
|
|
preps = load_preps(filename(prepath(sett), :preps), sett.autS)
|
|
|
|
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2018-07-31 10:21:54 +02:00
|
|
|
function load_preps(fname::String, G::Group)
|
2017-11-08 09:33:00 +01:00
|
|
|
lded_preps = load(fname, "perms_d")
|
|
|
|
permG = PermutationGroup(length(first(lded_preps)))
|
|
|
|
@assert length(lded_preps) == order(G)
|
|
|
|
return Dict(k=>permG(v) for (k,v) in zip(elements(G), lded_preps))
|
|
|
|
end
|
|
|
|
|
|
|
|
function save_preps(fname::String, preps)
|
|
|
|
autS = parent(first(keys(preps)))
|
2018-08-20 04:02:44 +02:00
|
|
|
save(fname, "perms_d", [preps[elt].d for elt in elements(autS)])
|
2017-11-08 09:33:00 +01:00
|
|
|
end
|
|
|
|
|
2018-08-20 03:57:33 +02:00
|
|
|
|
2018-01-02 03:19:57 +01:00
|
|
|
|
2018-09-05 09:06:00 +02:00
|
|
|
|
|
|
|
|
2018-01-02 03:19:57 +01:00
|
|
|
end
|
2018-01-01 14:06:33 +01:00
|
|
|
|
|
|
|
|
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|