mirror of
https://github.com/kalmarek/PropertyT.jl.git
synced 2024-12-25 18:25:30 +01:00
move application specific code to 1712.07167.jl
This commit is contained in:
parent
d6bb71b3cb
commit
38d80e63ed
248
src/1712.07167.jl
Normal file
248
src/1712.07167.jl
Normal file
@ -0,0 +1,248 @@
|
|||||||
|
###############################################################################
|
||||||
|
#
|
||||||
|
# Settings and filenames
|
||||||
|
#
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
struct Symmetrize end
|
||||||
|
struct Naive end
|
||||||
|
|
||||||
|
abstract type PropertyTSettings end
|
||||||
|
|
||||||
|
struct SolverSettings
|
||||||
|
sdpsolver::AbstractMathProgSolver
|
||||||
|
upper_bound::Float64
|
||||||
|
warmstart::Bool
|
||||||
|
|
||||||
|
SolverSettings(sol, ub, ws=true) = new(sol, upper_bound, ws)
|
||||||
|
end
|
||||||
|
|
||||||
|
struct Naive <: PropertyTSettings
|
||||||
|
name::String
|
||||||
|
G::Group
|
||||||
|
S::Vector{GroupElem}
|
||||||
|
radius::Int
|
||||||
|
|
||||||
|
solver::SolverSettings
|
||||||
|
end
|
||||||
|
|
||||||
|
struct Symmetrized <: PropertyTSettings
|
||||||
|
name::String
|
||||||
|
G::Group
|
||||||
|
S::Vector{GroupElem}
|
||||||
|
autS::Group
|
||||||
|
radius::Int
|
||||||
|
|
||||||
|
solver::SolverSettings
|
||||||
|
end
|
||||||
|
|
||||||
|
function Settings(name::String,
|
||||||
|
G::Group, S::Vector{GEl}, r::Integer,
|
||||||
|
sol::Solver, ub, ws=true) where {GEl<:GroupElem, Solver<:AbstractMathProgSolver}
|
||||||
|
sol_sett = SolverSettings(sol, ub, ws)
|
||||||
|
return Naive(name, G, S, r, sol_sett)
|
||||||
|
end
|
||||||
|
|
||||||
|
function Settings(name::String,
|
||||||
|
G::Group, S::Vector{GEl}, autS::Group, r::Integer,
|
||||||
|
sol::Solver, ub, ws=true) where {GEl<:GroupElem, Solver<:AbstractMathProgSolver}
|
||||||
|
sol_sett = SolverSettings(sol, ub, ws)
|
||||||
|
return Symmetrized(name, G, S, autS, r, sol_sett)
|
||||||
|
end
|
||||||
|
|
||||||
|
prefix(s::Naive) = s.name
|
||||||
|
prefix(s::Symmetrized) = "o"*s.name
|
||||||
|
suffix(s::PropertyTSettings) = "$(s.upper_bound)"
|
||||||
|
prepath(s::PropertyTSettings) = prefix(s)
|
||||||
|
fullpath(s::PropertyTSettings) = joinpath(prefix(s), suffix(s))
|
||||||
|
|
||||||
|
filename(sett::PropertyTSettings, s::Symbol) = filename(sett, Val{s})
|
||||||
|
|
||||||
|
filename(sett::PropertyTSettings, ::Type{Val{:fulllog}}) =
|
||||||
|
joinpath(fullpath(sett), "full_$(string(now())).log")
|
||||||
|
filename(sett::PropertyTSettings, ::Type{Val{:solverlog}}) =
|
||||||
|
joinpath(fullpath(sett), "solver_$(string(now())).log")
|
||||||
|
|
||||||
|
filename(sett::PropertyTSettings, ::Type{Val{:Δ}}) =
|
||||||
|
joinpath(prepath(sett), "delta.jld")
|
||||||
|
filename(sett::PropertyTSettings, ::Type{Val{:OrbitData}}) =
|
||||||
|
joinpath(prepath(sett), "OrbitData.jld")
|
||||||
|
|
||||||
|
filename(sett::PropertyTSettings, ::Type{Val{:warmstart}}) =
|
||||||
|
joinpath(fullpath(sett), "warmstart.jld")
|
||||||
|
filename(sett::PropertyTSettings, ::Type{Val{:solution}}) =
|
||||||
|
joinpath(fullpath(sett), "solution.jld")
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
#
|
||||||
|
# λandP
|
||||||
|
#
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
function warmstart(sett::PropertyTSettings)
|
||||||
|
if sett.solver.warmstart && isfile(filename(sett, :warmstart))
|
||||||
|
ws = load(filename(sett, :warmstart), "warmstart")
|
||||||
|
else
|
||||||
|
ws = nothing
|
||||||
|
end
|
||||||
|
return ws
|
||||||
|
end
|
||||||
|
|
||||||
|
function computeλandP(sett::Naive, Δ::GroupRingElem;
|
||||||
|
solverlog=tempname()*".log")
|
||||||
|
|
||||||
|
info("Creating SDP problem...")
|
||||||
|
SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, upper_bound=sett.solver.upper_bound)
|
||||||
|
JuMP.setsolver(SDP_problem, sett.solver.sdpsolver)
|
||||||
|
info(Base.repr(SDP_problem))
|
||||||
|
|
||||||
|
ws = warmstart(sett)
|
||||||
|
@time status, (λ, P, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws)
|
||||||
|
@show status
|
||||||
|
save(filename(sett, :warmstart), "warmstart", ws, "P", P, "λ", λ)
|
||||||
|
|
||||||
|
return λ, P
|
||||||
|
end
|
||||||
|
|
||||||
|
function computeλandP(sett::Symmetrized, Δ::GroupRingElem;
|
||||||
|
solverlog=tempname()*".log")
|
||||||
|
|
||||||
|
if !isfile(filename(sett, :OrbitData))
|
||||||
|
isdefined(parent(Δ), :basis) || throw("You need to define basis of Group Ring to compute orbit decomposition!")
|
||||||
|
orbit_data = OrbitData(parent(Δ), sett.autS)
|
||||||
|
save(filename(sett, :OrbitData), "OrbitData", orbit_data)
|
||||||
|
end
|
||||||
|
orbit_data = load(filename(sett, :OrbitData), "OrbitData")
|
||||||
|
orbit_data = decimate(orbit_data)
|
||||||
|
|
||||||
|
info("Creating SDP problem...")
|
||||||
|
|
||||||
|
SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, orbit_data, upper_bound=sett.solver.upper_bound)
|
||||||
|
JuMP.setsolver(SDP_problem, sett.solver.sdpsolver)
|
||||||
|
info(Base.repr(SDP_problem))
|
||||||
|
|
||||||
|
ws = warmstart(sett)
|
||||||
|
@time status, (λ, Ps, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws)
|
||||||
|
@show status
|
||||||
|
save(filename(sett, :warmstart), "warmstart", ws, "Ps", Ps, "λ", λ)
|
||||||
|
|
||||||
|
info("Reconstructing P...")
|
||||||
|
@time P = reconstruct(Ps, orbit_data)
|
||||||
|
|
||||||
|
return λ, P
|
||||||
|
end
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
#
|
||||||
|
# Checking solution
|
||||||
|
#
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
function distance_to_positive_cone(Δ::GroupRingElem, λ, Q; R::Int=2)
|
||||||
|
@info("------------------------------------------------------------")
|
||||||
|
@info("Checking in floating-point arithmetic...")
|
||||||
|
@info("λ = $λ")
|
||||||
|
eoi = Δ^2-λ*Δ
|
||||||
|
|
||||||
|
@time residual = eoi - compute_SOS(parent(eoi), augIdproj(Q))
|
||||||
|
@info("ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", aug(residual)))")
|
||||||
|
L1_norm = norm(residual,1)
|
||||||
|
@info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", L1_norm))")
|
||||||
|
|
||||||
|
distance = λ - 2.0^(2ceil(log2(R)))*L1_norm
|
||||||
|
|
||||||
|
@info("Floating point distance (to positive cone) ≈")
|
||||||
|
@info("$(@sprintf("%.10f", distance))")
|
||||||
|
@info("")
|
||||||
|
|
||||||
|
if distance ≤ 0
|
||||||
|
return distance
|
||||||
|
end
|
||||||
|
|
||||||
|
@info("------------------------------------------------------------")
|
||||||
|
@info("Checking in interval arithmetic...")
|
||||||
|
@info("λ ∈ $λ")
|
||||||
|
|
||||||
|
λ = @interval(λ)
|
||||||
|
eoi = Δ^2 - λ*Δ
|
||||||
|
|
||||||
|
@info("Projecting columns of Q to the augmentation ideal...")
|
||||||
|
@time Q, check = augIdproj(Interval, Q)
|
||||||
|
@info("Checking that sum of every column contains 0.0... ")
|
||||||
|
@info((check? "They do." : "FAILED!"))
|
||||||
|
check || @warn("The following numbers are meaningless!")
|
||||||
|
|
||||||
|
@time residual = eoi - compute_SOS(parent(eoi), Q)
|
||||||
|
@info("ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ∈ $(aug(residual))")
|
||||||
|
L1_norm = norm(residual,1)
|
||||||
|
@info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)")
|
||||||
|
|
||||||
|
distance = λ - 2.0^(2ceil(log2(R)))*L1_norm
|
||||||
|
|
||||||
|
@info("Interval distance (to positive cone) ∈")
|
||||||
|
@info("$(distance)")
|
||||||
|
@info("")
|
||||||
|
|
||||||
|
return distance.lo
|
||||||
|
end
|
||||||
|
|
||||||
|
###############################################################################
|
||||||
|
#
|
||||||
|
# Interpreting the numerical results
|
||||||
|
#
|
||||||
|
###############################################################################
|
||||||
|
|
||||||
|
Kazhdan(λ::Number, N::Integer) = sqrt(2*λ/N)
|
||||||
|
|
||||||
|
function interpret_results(sett::PropertyTSettings, sgap::Number)
|
||||||
|
|
||||||
|
if sgap > 0
|
||||||
|
Kazhdan_κ = Kazhdan(sgap, length(sett.S))
|
||||||
|
if Kazhdan_κ > 0
|
||||||
|
info("κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
||||||
|
return true
|
||||||
|
end
|
||||||
|
end
|
||||||
|
info("λ($(sett.name), S) ≥ $sgap < 0: Tells us nothing about property (T)")
|
||||||
|
return false
|
||||||
|
end
|
||||||
|
|
||||||
|
function check_property_T(sett::PropertyTSettings)
|
||||||
|
fp = PropertyT.fullpath(sett)
|
||||||
|
isdir(fp) || mkpath(fp)
|
||||||
|
|
||||||
|
if isfile(filename(sett,:Δ))
|
||||||
|
# cached
|
||||||
|
Δ = loadLaplacian(filename(sett,:Δ), sett.G)
|
||||||
|
else
|
||||||
|
# compute
|
||||||
|
Δ = Laplacian(sett.S, sett.radius)
|
||||||
|
saveLaplacian(filename(sett, :Δ), Δ)
|
||||||
|
end
|
||||||
|
|
||||||
|
if !sett.warmstart && isfile(filename(sett, :solution))
|
||||||
|
λ, P = load(filename(sett, :solution), "λ", "P")
|
||||||
|
else
|
||||||
|
λ, P = computeλandP(sett, Δ,
|
||||||
|
solverlog=filename(sett, :solverlog))
|
||||||
|
|
||||||
|
save(filename(sett, :solution), "λ", λ, "P", P)
|
||||||
|
|
||||||
|
if λ < 0
|
||||||
|
warn("Solver did not produce a valid solution!")
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
info("λ = $λ")
|
||||||
|
info("sum(P) = $(sum(P))")
|
||||||
|
info("maximum(P) = $(maximum(P))")
|
||||||
|
info("minimum(P) = $(minimum(P))")
|
||||||
|
|
||||||
|
isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) ||
|
||||||
|
warn("The solution matrix doesn't seem to be positive definite!")
|
||||||
|
|
||||||
|
@time Q = real(sqrtm((P+P')/2))
|
||||||
|
sgap = distance_to_positive_cone(Δ, λ, Q, wlen=2*sett.radius)
|
||||||
|
|
||||||
|
return interpret_results(sett, sgap)
|
||||||
|
end
|
116
src/PropertyT.jl
116
src/PropertyT.jl
@ -12,126 +12,10 @@ using JuMP
|
|||||||
|
|
||||||
import MathProgBase.SolverInterface.AbstractMathProgSolver
|
import MathProgBase.SolverInterface.AbstractMathProgSolver
|
||||||
|
|
||||||
###############################################################################
|
|
||||||
#
|
|
||||||
# Settings and filenames
|
|
||||||
#
|
|
||||||
###############################################################################
|
|
||||||
|
|
||||||
struct Symmetrize end
|
|
||||||
struct Naive end
|
|
||||||
|
|
||||||
struct Settings{T, GEl<:GroupElem}
|
|
||||||
name::String
|
|
||||||
|
|
||||||
G::Group
|
|
||||||
S::Vector{GEl}
|
|
||||||
radius::Int
|
|
||||||
|
|
||||||
solver::AbstractMathProgSolver
|
|
||||||
upper_bound::Float64
|
|
||||||
tol::Float64
|
|
||||||
warmstart::Bool
|
|
||||||
|
|
||||||
autS::Group
|
|
||||||
|
|
||||||
function Settings(name::String,
|
|
||||||
G::Group, S::Vector{GEl}, r::Int,
|
|
||||||
sol::Sol, ub, tol, ws) where {GEl<:GroupElem, Sol<:AbstractMathProgSolver}
|
|
||||||
return new{Naive, GEl}(name, G, S, r, sol, ub, tol, ws)
|
|
||||||
end
|
|
||||||
|
|
||||||
function Settings(name::String,
|
|
||||||
G::Group, S::Vector{GEl}, r::Int,
|
|
||||||
sol::Sol, ub, tol, ws, autS) where {GEl<:GroupElem, Sol<:AbstractMathProgSolver}
|
|
||||||
return new{Symmetrize, GEl}(name, G, S, r, sol, ub, tol, ws, autS)
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
prefix(s::Settings{Naive}) = s.name
|
|
||||||
prefix(s::Settings{Symmetrize}) = "o"*s.name
|
|
||||||
suffix(s::Settings) = "$(s.upper_bound)"
|
|
||||||
prepath(s::Settings) = prefix(s)
|
|
||||||
fullpath(s::Settings) = joinpath(prefix(s), suffix(s))
|
|
||||||
|
|
||||||
filename(sett::Settings, s::Symbol) = filename(sett, Val{s})
|
|
||||||
|
|
||||||
filename(sett::Settings, ::Type{Val{:fulllog}}) =
|
|
||||||
joinpath(fullpath(sett), "full_$(string(now())).log")
|
|
||||||
filename(sett::Settings, ::Type{Val{:solverlog}}) =
|
|
||||||
joinpath(fullpath(sett), "solver_$(string(now())).log")
|
|
||||||
|
|
||||||
filename(sett::Settings, ::Type{Val{:Δ}}) =
|
|
||||||
joinpath(prepath(sett), "delta.jld")
|
|
||||||
filename(sett::Settings, ::Type{Val{:OrbitData}}) =
|
|
||||||
joinpath(prepath(sett), "OrbitData.jld")
|
|
||||||
|
|
||||||
filename(sett::Settings, ::Type{Val{:warmstart}}) =
|
|
||||||
joinpath(fullpath(sett), "warmstart.jld")
|
|
||||||
filename(sett::Settings, ::Type{Val{:solution}}) =
|
|
||||||
joinpath(fullpath(sett), "solution.jld")
|
|
||||||
|
|
||||||
function check_property_T(sett::Settings)
|
|
||||||
fp = PropertyT.fullpath(sett)
|
|
||||||
isdir(fp) || mkpath(fp)
|
|
||||||
|
|
||||||
if isfile(filename(sett,:Δ))
|
|
||||||
# cached
|
|
||||||
Δ = loadLaplacian(filename(sett,:Δ), sett.G)
|
|
||||||
else
|
|
||||||
# compute
|
|
||||||
Δ = Laplacian(sett.S, sett.radius)
|
|
||||||
saveLaplacian(filename(sett, :Δ), Δ)
|
|
||||||
end
|
|
||||||
|
|
||||||
if !sett.warmstart && isfile(filename(sett, :solution))
|
|
||||||
λ, P = load(filename(sett, :solution), "λ", "P")
|
|
||||||
else
|
|
||||||
λ, P = computeλandP(sett, Δ,
|
|
||||||
solverlog=filename(sett, :solverlog))
|
|
||||||
|
|
||||||
save(filename(sett, :solution), "λ", λ, "P", P)
|
|
||||||
|
|
||||||
if λ < 0
|
|
||||||
warn("Solver did not produce a valid solution!")
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
info("λ = $λ")
|
|
||||||
info("sum(P) = $(sum(P))")
|
|
||||||
info("maximum(P) = $(maximum(P))")
|
|
||||||
info("minimum(P) = $(minimum(P))")
|
|
||||||
|
|
||||||
isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) ||
|
|
||||||
warn("The solution matrix doesn't seem to be positive definite!")
|
|
||||||
|
|
||||||
@time Q = real(sqrtm((P+P')/2))
|
|
||||||
sgap = distance_to_cone(Δ, λ, Q, wlen=2*sett.radius)
|
|
||||||
|
|
||||||
return interpret_results(sett, sgap)
|
|
||||||
end
|
|
||||||
|
|
||||||
Kazhdan(λ::Number, N::Integer) = sqrt(2*λ/N)
|
|
||||||
|
|
||||||
function interpret_results(sett::Settings, sgap::Number)
|
|
||||||
|
|
||||||
if sgap > 0
|
|
||||||
Kazhdan_κ = Kazhdan(sgap, length(sett.S))
|
|
||||||
if Kazhdan_κ > 0
|
|
||||||
info("κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
|
||||||
return true
|
|
||||||
end
|
|
||||||
end
|
|
||||||
info("λ($(sett.name), S) ≥ $sgap < 0: Tells us nothing about property (T)")
|
|
||||||
return false
|
|
||||||
end
|
|
||||||
|
|
||||||
include("laplacians.jl")
|
include("laplacians.jl")
|
||||||
include("RGprojections.jl")
|
include("RGprojections.jl")
|
||||||
include("orbitdata.jl")
|
include("orbitdata.jl")
|
||||||
include("sos_sdps.jl")
|
include("sos_sdps.jl")
|
||||||
include("checksolution.jl")
|
include("checksolution.jl")
|
||||||
|
|
||||||
|
|
||||||
end # module Property(T)
|
end # module Property(T)
|
||||||
|
@ -28,52 +28,5 @@ function augIdproj(Q::AbstractArray{T,2}) where {T<:Real}
|
|||||||
return result
|
return result
|
||||||
end
|
end
|
||||||
|
|
||||||
function distance_to_cone(Δ::GroupRingElem, λ, Q; wlen::Int=4)
|
|
||||||
info("------------------------------------------------------------")
|
|
||||||
info("Checking in floating-point arithmetic...")
|
|
||||||
info("λ = $λ")
|
|
||||||
@time sos = compute_SOS(parent(Δ), Q)
|
|
||||||
residue = Δ^2-λ*Δ - sos
|
|
||||||
info("ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", aug(residue)))")
|
|
||||||
L1_norm = norm(residue,1)
|
|
||||||
info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", L1_norm))")
|
|
||||||
|
|
||||||
distance = λ - 2^(wlen-1)*L1_norm
|
|
||||||
|
|
||||||
info("Floating point distance (to positive cone) ≈")
|
|
||||||
info("$(@sprintf("%.10f", distance))")
|
|
||||||
info("")
|
|
||||||
|
|
||||||
if distance ≤ 0
|
|
||||||
return distance
|
|
||||||
end
|
end
|
||||||
|
|
||||||
info("------------------------------------------------------------")
|
|
||||||
info("Checking in interval arithmetic...")
|
|
||||||
info("λ ∈ $λ")
|
|
||||||
|
|
||||||
λ = @interval(λ)
|
|
||||||
eoi = Δ^2 - λ*Δ
|
|
||||||
info("Projecting columns of Q to the augmentation ideal...")
|
|
||||||
T = eltype(Q)
|
|
||||||
@time Q = augIdproj(Q)
|
|
||||||
|
|
||||||
info("Checking that sum of every column contains 0.0... ")
|
|
||||||
check = all([zero(T) in sum(view(Q, :, i)) for i in 1:size(Q, 2)])
|
|
||||||
info((check? "They do." : "FAILED!"))
|
|
||||||
|
|
||||||
@assert check
|
|
||||||
|
|
||||||
@time sos = compute_SOS(parent(Δ), Q)
|
|
||||||
residue = Δ^2-λ*Δ - sos
|
|
||||||
info("ɛ(∑ξᵢ*ξᵢ) ∈ $(aug(residue))")
|
|
||||||
L1_norm = norm(residue,1)
|
|
||||||
info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)")
|
|
||||||
|
|
||||||
distance = λ - 2^(wlen-1)*L1_norm
|
|
||||||
info("The Augmentation-projected distance (to positive cone) ∈")
|
|
||||||
info("$(distance)")
|
|
||||||
info("")
|
|
||||||
|
|
||||||
return distance.lo
|
|
||||||
end
|
end
|
||||||
|
@ -61,61 +61,3 @@ function loadGRElem(fname::String, G::Group)
|
|||||||
end
|
end
|
||||||
return Δ
|
return Δ
|
||||||
end
|
end
|
||||||
|
|
||||||
###############################################################################
|
|
||||||
#
|
|
||||||
# λandP
|
|
||||||
#
|
|
||||||
###############################################################################
|
|
||||||
|
|
||||||
function computeλandP(sett::Settings{Naive}, Δ::GroupRingElem;
|
|
||||||
solverlog=tempname()*".log")
|
|
||||||
|
|
||||||
info("Creating SDP problem...")
|
|
||||||
SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, upper_bound=sett.upper_bound)
|
|
||||||
JuMP.setsolver(SDP_problem, sett.solver)
|
|
||||||
info(Base.repr(SDP_problem))
|
|
||||||
|
|
||||||
ws = warmstart(sett)
|
|
||||||
@time status, (λ, P, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws)
|
|
||||||
save(filename(sett, :warmstart), "warmstart", ws)
|
|
||||||
|
|
||||||
return λ, P
|
|
||||||
end
|
|
||||||
|
|
||||||
function computeλandP(sett::Settings{Symmetrize}, Δ::GroupRingElem;
|
|
||||||
solverlog=tempname()*".log")
|
|
||||||
|
|
||||||
if isfile(filename(sett, :OrbitData))
|
|
||||||
orbit_data = load(filename(sett, :OrbitData), "OrbitData")
|
|
||||||
else
|
|
||||||
isdefined(parent(Δ), :basis) || throw("You need to define basis of Group Ring to compute orbit decomposition!")
|
|
||||||
orbit_data = OrbitData(parent(Δ), sett.autS)
|
|
||||||
save(filename(sett, :OrbitData), "OrbitData", orbit_data)
|
|
||||||
end
|
|
||||||
orbit_data = decimate(orbit_data)
|
|
||||||
|
|
||||||
info("Creating SDP problem...")
|
|
||||||
|
|
||||||
SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, orbit_data, upper_bound=sett.upper_bound)
|
|
||||||
JuMP.setsolver(SDP_problem, sett.solver)
|
|
||||||
info(Base.repr(SDP_problem))
|
|
||||||
|
|
||||||
ws = warmstart(sett)
|
|
||||||
@time status, (λ, Ps, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws)
|
|
||||||
save(filename(sett, :warmstart), "warmstart", ws, "Ps", Ps, "λ", λ)
|
|
||||||
|
|
||||||
info("Reconstructing P...")
|
|
||||||
@time P = reconstruct(Ps, orbit_data)
|
|
||||||
|
|
||||||
return λ, P
|
|
||||||
end
|
|
||||||
|
|
||||||
function warmstart(sett::Settings)
|
|
||||||
if sett.warmstart && isfile(filename(sett, :warmstart))
|
|
||||||
ws = load(filename(sett, :warmstart), "warmstart")
|
|
||||||
else
|
|
||||||
ws = nothing
|
|
||||||
end
|
|
||||||
return ws
|
|
||||||
end
|
|
||||||
|
Loading…
Reference in New Issue
Block a user