move SDP-related functions together

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
 
 

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)

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
 
 ###############################################################################