From b9de894fb1471b707a2d58e4fbfb05019d6e8319 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 5 Sep 2018 09:18:38 +0200 Subject: [PATCH] move SDP-related functions together --- src/Orbit-wise.jl | 43 ++------------------- src/OrbitDecomposition.jl | 22 ----------- src/SDPs.jl | 79 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 83 insertions(+), 61 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 9fe7395..fb48ef9 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -33,53 +33,18 @@ 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))) + + + + end -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) - end end -function init_model(m, sizes) - P = Vector{Array{JuMP.Variable,2}}(length(sizes)) - 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 - return P -end - -function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData; upper_bound=Inf) - m = JuMP.Model(); - P = init_model(m, size.(data.Uπs,2)) - - λ = JuMP.@variable(m, λ) - if upper_bound < Inf - JuMP.@constraint(m, λ <= upper_bound) - end - - info("Adding $(length(data.orbits)) constraints... ") - - @time addconstraints!(m, X, orderunit, λ, P, data) - - JuMP.@objective(m, Max, λ) - return m, λ, P end diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index b150f25..d3505da 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -96,28 +96,6 @@ function perm_reps(S::Vector, autS::Group, radius::Int) return perm_reps(autS, E) end -function reconstruct_sol(preps::Dict{T, S}, Us::Vector, Ps::Vector, dims::Vector) where {T<:GroupElem, S<:perm} - - 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 - - recP = 1/length(perms) .* sum(tmp) - for i in eachindex(recP) - if abs(recP[i]) .< eps(eltype(recP))*100 - recP[i] = zero(eltype(recP)) - end - end - return recP -end - function Cstar_repr(x::GroupRingElem{T}, mreps::Dict) where {T} nzeros = findn(x.coeffs) return sum(x[i].*mreps[parent(x).basis[i]] for i in nzeros) diff --git a/src/SDPs.jl b/src/SDPs.jl index 192f38b..c706602 100644 --- a/src/SDPs.jl +++ b/src/SDPs.jl @@ -41,10 +41,89 @@ function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem; upper_bound=Inf return m, λ, P end +############################################################################### +# +# Symmetrized SDP +# +############################################################################### +function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData; upper_bound=Inf) + m = JuMP.Model(); + + sizes = size.(data.Uπs, 2) + P = Vector{Matrix{JuMP.Variable}}(length(sizes)) + + 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, λ) + if upper_bound < Inf + JuMP.@constraint(m, λ <= upper_bound) + end + info("Adding $(length(data.orbits)) constraints... ") + + addconstraints!(m,P,λ,X,orderunit, data) + + JuMP.@objective(m, Max, λ) + return m, λ, P +end + +function constraintLHS!(M, cnstr, Us, Ust, dims, eps=1000*eps(eltype(first(M)))) + for π in 1:endof(Us) + M[π] = PropertyT.sparsify!(dims[π].*Ust[π]*cnstr*Us[π], eps) + end +end + +function addconstraints!(m::JuMP.Model, + P::Vector{Matrix{JuMP.Variable}}, λ::JuMP.Variable, + X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData) + + orderunit_orb = orbit_spvector(orderunit.coeffs, data.orbits) + X_orb = orbit_spvector(X.coeffs, data.orbits) + UπsT = [U' for U in data.Uπs] + + cnstrs = constraints(parent(X).pm) + orb_cnstr = spzeros(Float64, size(parent(X).pm)...) + + M = [Array{Float64}(n,n) for n in size.(UπsT,1)] + + for t in 1:length(X_orb) + orbit_constraint!(orb_cnstr, cnstrs[data.orbits[t]]) + constraintLHS!(M, orb_cnstr, data.Uπs, UπsT, data.dims) + + lhs = @expression(m, sum(vecdot(M[π], P[π]) for π in 1:endof(data.Uπs))) + x, u = X_orb[t], orderunit_orb[t] + JuMP.@constraint(m, lhs == x - λ*u) + end +end + +function reconstruct(Ps::Vector{Matrix{F}}, data::OrbitData) where F + return reconstruct(Ps, data.preps, data.Uπs, data.dims) +end + +function reconstruct(Ps::Vector{M}, preps::Dict{GEl, P}, Uπs::Vector{U}, dims::Vector{Int}) where {M<:AbstractMatrix, GEl<:GroupElem, P<:perm, U<:AbstractMatrix} + + l = length(Uπs) + transfP = [dims[π].*Uπs[π]*Ps[π]*Uπs[π]' for π in 1:l] + tmp = [zeros(Float64, size(first(transfP))) for _ in 1:l] + perms = collect(keys(preps)) + + 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) + # for i in eachindex(recP) + # if abs(recP[i]) .< eps(eltype(recP))*100 + # recP[i] = zero(eltype(recP)) + # end + # end + return recP end ###############################################################################