diff --git a/positivity/check_positivity.jl b/positivity/check_positivity.jl index c5b4a74..d24c04a 100644 --- a/positivity/check_positivity.jl +++ b/positivity/check_positivity.jl @@ -129,40 +129,48 @@ orbit_data = PropertyT.decimate(orbit_data); elt = adj+3op; -SDP_problem, varλ, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=sett.upper_bound) +const SOLUTION_FILE = PropertyT.filename(sett, :solution) -begin - using SCS - scs_solver = SCS.SCSSolver(linearsolver=SCS.Direct, - eps=sett.tol, - max_iters=args["iterations"], - alpha=1.95, - acceleration_lookback=1) +if !isfile(SOLUTION_FILE) + + SDP_problem, varλ, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=sett.upper_bound) - JuMP.setsolver(SDP_problem, scs_solver) -end + begin + using SCS + scs_solver = SCS.SCSSolver(linearsolver=SCS.Direct, + eps=sett.tol, + max_iters=args["iterations"], + alpha=1.95, + acceleration_lookback=1) -λ = Ps = nothing -ws = PropertyT.warmstart(sett) - -# using ProgressMeter - -# @showprogress "Running SDP optimization step... " for i in 1:args["repetitions"] -while true - λ, Ps, ws = PropertyT.solve(PropertyT.filename(sett, :solverlog), - SDP_problem, varλ, varP, ws); - - if all((!isnan).(ws[1])) - save(WARMSTART_FILE, "warmstart", ws, "λ", λ, "Ps", Ps) - save(WARMSTART_FILE[1:end-4]*"_$(now()).jld", "warmstart", ws, "λ", λ, "Ps", Ps) - else - warn("No valid solution was saved!") + JuMP.setsolver(SDP_problem, scs_solver) end + + λ = Ps = nothing + ws = PropertyT.warmstart(sett) + + # using ProgressMeter + + # @showprogress "Running SDP optimization step... " for i in 1:args["repetitions"] + while true + λ, Ps, ws = PropertyT.solve(PropertyT.filename(sett, :solverlog), + SDP_problem, varλ, varP, ws); + + if all((!isnan).(ws[1])) + save(WARMSTART_FILE, "warmstart", ws, "λ", λ, "Ps", Ps) + save(WARMSTART_FILE[1:end-4]*"_$(now()).jld", "warmstart", ws, "λ", λ, "Ps", Ps) + else + warn("No valid solution was saved!") + end + end + + info("Reconstructing P...") + @time P = PropertyT.reconstruct(Ps, orbit_data); + save(SOLUTION_FILE, "λ", λ, "P", P) end -info("Reconstructing P...") -@time P = PropertyT.reconstruct(Ps, orbit_data); -save(PropertyT.filename(sett, :solution), "λ", λ, "P", P) +P, λ = load(SOLUTION_FILE, "P", "λ") +@show λ; @time const Q = real(sqrtm(P)); @@ -172,14 +180,12 @@ function SOS_residual(eoi::GroupRingElem, Q::Matrix) return eoi - sos end -addprocs(4) -@everywhere using PropertyT - -@show λ; +info("Floating Point arithmetic:") EOI = elt - λ*Δ b = SOS_residual(EOI, Q) @show norm(b, 1); +info("Interval arithmetic:") using IntervalArithmetic Qint = PropertyT.augIdproj(Q); @assert all([zero(eltype(Q)) in sum(view(Qint, :, i)) for i in 1:size(Qint, 2)]) @@ -190,4 +196,4 @@ Q_int = PropertyT.augIdproj(Q); b_int = SOS_residual(EOI_int, Q_int) @show norm(b_int, 1); -@interval(λ) - 2^2*norm(b_int,1) +info("λ is certified to be > ", (@interval(λ) - 2^2*norm(b_int,1)).lo)