update check_positivity

This commit is contained in:
kalmarek 2018-11-17 22:42:42 +01:00
parent 523d783614
commit bf8fc67bb4

View File

@ -13,7 +13,7 @@ include("../main.jl")
using PropertyTGroups using PropertyTGroups
args = Dict("SAut" => 5, "upper-bound" => 10.0, "radius" => 2, "nosymmetry"=>false, "tol"=>1e-10, "iterations" =>200000, "repetitions"=>5) args = Dict("SAut" => 5, "upper-bound" => 50.0, "radius" => 2, "nosymmetry"=>false, "tol"=>1e-12, "iterations" =>200000, "warmstart" => true)
Gr = PropertyTGroups.PropertyTGroup(args) Gr = PropertyTGroups.PropertyTGroup(args)
sett = PropertyT.Settings(Gr, args, sett = PropertyT.Settings(Gr, args,
@ -25,19 +25,6 @@ fullpath = PropertyT.fullpath(sett)
isdir(fullpath) || mkpath(fullpath) isdir(fullpath) || mkpath(fullpath)
# setup_logging(PropertyT.filename(fullpath, :fulllog), :fulllog) # setup_logging(PropertyT.filename(fullpath, :fulllog), :fulllog)
function (p::perm)(A::GroupRingElem)
RG = parent(A)
T = eltype(A.coeffs)
result = zero(RG, T)
for (idx, c) in enumerate(A.coeffs)
if c!= zero(T)
result[p(RG.basis[idx])] = c
end
end
return result
end
function small_generating_set(RG::GroupRing{AutGroup{N}}, n) where N function small_generating_set(RG::GroupRing{AutGroup{N}}, n) where N
indexing = [(i,j) for i in 1:n for j in 1:n if i≠j] indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
@ -158,10 +145,10 @@ end
λ = Ps = nothing λ = Ps = nothing
ws = PropertyT.warmstart(sett) ws = PropertyT.warmstart(sett)
using ProgressMeter # using ProgressMeter
@showprogress "Running SDP optimization step... " for i in 1:args["repetitions"] # @showprogress "Running SDP optimization step... " for i in 1:args["repetitions"]
# while true while true
λ, Ps, ws = PropertyT.solve(PropertyT.filename(sett, :solverlog), λ, Ps, ws = PropertyT.solve(PropertyT.filename(sett, :solverlog),
SDP_problem, varλ, varP, ws); SDP_problem, varλ, varP, ws);
@ -173,28 +160,34 @@ using ProgressMeter
end end
end end
function check_SOS_precision(Q::Matrix, eoi::GroupRingElem)
RG = parent(eoi)
@time sos = PropertyT.compute_SOS(RG, Q);
residue = eoi - sos
return norm(residue, 1)
end
info("Reconstructing P...") info("Reconstructing P...")
@time P = PropertyT.reconstruct(Ps, orbit_data); @time P = PropertyT.reconstruct(Ps, orbit_data);
save(PropertyT.filename(sett, :solution), "λ", λ, "P", P) save(PropertyT.filename(sett, :solution), "λ", λ, "P", P)
@time const Q = real(sqrtm(P));
function SOS_residual(eoi::GroupRingElem, Q::Matrix)
RG = parent(eoi)
@time sos = PropertyT.compute_SOS(RG, Q);
return eoi - sos
end
addprocs(4) addprocs(4)
@everywhere using PropertyT @everywhere using PropertyT
@time const Q = real(sqrtm(P)); @show λ;
const EOI = elt - λ*Δ; EOI = elt - λ*Δ
@show λ - 2^3*check_SOS_precision(Q, EOI); b = SOS_residual(EOI, Q)
@show norm(b, 1);
using IntervalArithmetic using IntervalArithmetic
Qint = PropertyT.augIdproj(Q); Qint = PropertyT.augIdproj(Q);
@assert all([zero(eltype(Q)) in sum(view(Qint, :, i)) for i in 1:size(Qint, 2)]) @assert all([zero(eltype(Q)) in sum(view(Qint, :, i)) for i in 1:size(Qint, 2)])
EOIint = GroupRingElem([@interval(c) for c in EOI.coeffs], parent(Δ)); EOI_int = elt - @interval(λ)*Δ;
b_int = check_SOS_precision(Qint, EOIint) Q_int = PropertyT.augIdproj(Q);
@show @interval(λ) - 2^3*b_int; @assert all([zero(eltype(Q)) in sum(view(Q_int, :, i)) for i in 1:size(Q_int, 2)])
b_int = SOS_residual(EOI_int, Q_int)
@show norm(b_int, 1);
@interval(λ) - 2^2*norm(b_int,1)