kalmar/PropertyT.jl
1
0
Fork 0
mirror of https://github.com/kalmarek/PropertyT.jl.git synced 2024-11-19 07:20:28 +01:00
Code Issues Releases Wiki Activity

Change names: kappa -> lambda, A -> P

Browse Source
This commit is contained in:
kalmar 2017-04-01 15:21:57 +02:00
parent 028979bfd1
commit 5876998cba
3 changed files with 86 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

83
src/PropertyT.jl
View File

@ -20,14 +20,14 @@ function pmΔfilenames(name::String)
return pm_filename, Δ_coeff_filename return pm_filename, Δ_coeff_filename
end end
function κSDPfilenames(name::String) function λSDPfilenames(name::String)
if !isdir(name) if !isdir(name)
mkdir(name) mkdir(name)
end end
prefix = name prefix = name
κ_filename = joinpath(prefix, "kappa.jld") λ_filename = joinpath(prefix, "lambda.jld")
SDP_filename = joinpath(prefix, "SDPmatrixA.jld") SDP_filename = joinpath(prefix, "SDPmatrix.jld")
return κ_filename, SDP_filename return λ_filename, SDP_filename
end end
function ΔandSDPconstraints(name::String) function ΔandSDPconstraints(name::String)
@ -63,18 +63,19 @@ function ΔandSDPconstraints(name::String, generating_set::Function, radius::Int
end end
end end
function κandA(name::String) function λandP(name::String)
κ_fname, SDP_fname = κSDPfilenames(name) λ_fname, SDP_fname = λSDPfilenames(name)
f₁ = isfile(κ_fname) f₁ = isfile(λ_fname)
f₂ = isfile(SDP_fname) f₂ = isfile(SDP_fname)
if f₁ && f₂ if f₁ && f₂
info(logger, "Loading precomputed κ, A...") info(logger, "Loading precomputed λ, P...")
κ = load(κ_fname, "κ") λ = load(λ_fname, "λ")
A = load(SDP_fname, "A") P = load(SDP_fname, "P")
else else
throw(ArgumentError("You need to precompute κ and SDP matrix A to load it!")) throw(ArgumentError("You need to precompute λ and SDP matrix P to load it!"))
end end
return κ, A return λ, P
end end
function timed_msg(t) function timed_msg(t)
@ -87,32 +88,32 @@ function timed_msg(t)
end end
function κandA(name::String, opts...) function λandP(name::String, opts...)
try try
return κandA(name) return λandP(name)
catch err catch err
if isa(err, ArgumentError) if isa(err, ArgumentError)
if isfile("$name/solver.log") if isfile(joinpath(name, "solver.log"))
rm("$name/solver.log") rm(joinpath(name, "solver.log"))
end end
add_handler(solver_logger, DefaultHandler("./$name/solver.log", DefaultFormatter("{date}| {msg}")), "solver_log") add_handler(solver_logger, DefaultHandler(joinpath(name, "solver.log"), DefaultFormatter("{date}| {msg}")), "solver_log")
info(logger, "Creating SDP problem...") info(logger, "Creating SDP problem...")
κ, A = compute_κandA(opts...) λ, P = compute_λandP(opts...)
remove_handler(solver_logger, "solver_log") remove_handler(solver_logger, "solver_log")
κ_fname, A_fname = κSDPfilenames(name) λ_fname, P_fname = λSDPfilenames(name)
if κ > 0 if λ > 0
save(κ_fname, "κ", κ) save(λ_fname, "λ", λ)
save(A_fname, "A", A) save(P_fname, "P", P)
else else
throw(ErrorException("Solver $solver did not produce a valid solution!: κ = $κ")) throw(ErrorException("Solver $solver did not produce a valid solution!: λ = $λ"))
end end
return κ, A return λ, P
else else
# throw(err) # throw(err)
@ -121,21 +122,21 @@ function κandA(name::String, opts...)
end end
end end
function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver, upper_bound=Inf) function compute_λandP(sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver, upper_bound=Inf)
t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound) t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound)
info(logger, timed_msg(t)) info(logger, timed_msg(t))
κ = 0.0 λ = 0.0
A = nothing P = nothing
while κ == 0.0 while λ == 0.0
try try
κ, A = solve_SDP(SDP_problem, solver) λ, P = solve_SDP(SDP_problem, solver)
catch y catch y
warn(solver_logger, y) warn(solver_logger, y)
end end
end end
return κ, A return λ, P
end end
Kazhdan_from_sgap(λ,N) = sqrt(2*λ/N) Kazhdan_from_sgap(λ,N) = sqrt(2*λ/N)
@ -160,29 +161,29 @@ function check_property_T(name::String, generating_set::Function,
info(logger, "length(Δ) = $(length(Δ))") info(logger, "length(Δ) = $(length(Δ))")
info(logger, "size(Δ.product_matrix) = $(size(Δ.product_matrix))") info(logger, "size(Δ.product_matrix) = $(size(Δ.product_matrix))")
κ, A = κandA(name, sdp_constraints, Δ, solver, upper_bound) λ, P = λandP(name, sdp_constraints, Δ, solver, upper_bound)
info(logger, "κ = $κ") info(logger, "λ = $λ")
info(logger, "sum(A) = $(sum(A))") info(logger, "sum(P) = $(sum(P))")
info(logger, "maximum(A) = $(maximum(A))") info(logger, "maximum(P) = $(maximum(P))")
info(logger, "minimum(A) = $(minimum(A))") info(logger, "minimum(P) = $(minimum(P))")
if κ > 0 if λ > 0
spectral_gap = check_distance_to_positive_cone(Δ, κ, A, tol=tol, rational=false) spectral_gap = check_distance_to_positive_cone(Δ, λ, P, tol=tol, rational=false)
if isa(spectral_gap, Interval) if isa(spectral_gap, Interval)
spectral_gap = spectral_gap.lo spectral_gap = spectral_gap.lo
end end
if spectral_gap > 0 if spectral_gap > 0
@show spectral_gap @show spectral_gap
Kazhdan_const = Kazhdan_from_sgap(spectral_gap, S) Kazhdan_κ = Kazhdan_from_sgap(spectral_gap, S)
Kazhdan_const = Float64(trunc(Kazhdan_const, 12)) Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
info(logger, "κ($name, S) ≥ $Kazhdan_const: Group HAS property (T)!") info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!")
else else
spectral_gap = Float64(trunc(spectral_gap, 12)) spectral_gap = Float64(trunc(spectral_gap, 12))
info(logger, "λ($name, S) ≥ $spectral_gap: Group may NOT HAVE property (T)!") info(logger, "λ($name, S) ≥ $spectral_gap: Group may NOT HAVE property (T)!")
end end
else else
info(logger, "κ($name, S) ≥ $κ < 0: Tells us nothing about property (T)") info(logger, "κ($name, S) ≥ $λ < 0: Tells us nothing about property (T)")
end end
end end

68
src/checksolution.jl
View File

@ -6,8 +6,8 @@ ValidatedNumerics.setrounding(Interval, :correct)
ValidatedNumerics.setformat(:standard) ValidatedNumerics.setformat(:standard)
# setprecision(Interval, 53) # slightly faster than 256 # setprecision(Interval, 53) # slightly faster than 256
function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, κ::T) function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, λ::T)
return Δ*Δ - κ*Δ return Δ*Δ - λ*Δ
end end
function algebra_square(vector, elt) function algebra_square(vector, elt)
@ -62,69 +62,69 @@ end
ℚ(x, tol::Real) = rationalize(BigInt, x, tol=tol) ℚ(x, tol::Real) = rationalize(BigInt, x, tol=tol)
function distance_to_cone{T<:Rational}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T}) function distance_to_cone{T<:Rational}(λ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T})
SOS = compute_SOS(sqrt_matrix, Δ) SOS = compute_SOS(sqrt_matrix, Δ)
SOS_diff = EOI(Δ, κ) - SOS SOS_diff = EOI(Δ, λ) - SOS
eoi_SOS_L₁_dist = norm(SOS_diff,1) eoi_SOS_L₁_dist = norm(SOS_diff,1)
info(logger, "κ = $κ (≈$(@sprintf("%.10f", float(κ)))") info(logger, "λ = $λ (≈$(@sprintf("%.10f", float(λ)))")
ɛ_dist = GroupAlgebras.ɛ(SOS_diff) ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
if ɛ_dist ≠ 0//1 if ɛ_dist ≠ 0//1
warn(logger, "The SOS is not in the augmentation ideal, number below are meaningless!") warn(logger, "The SOS is not in the augmentation ideal, number below are meaningless!")
end end
info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) = $ɛ_dist") info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) = $ɛ_dist")
info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ = $(@sprintf("%.10f", float(eoi_SOS_L₁_dist)))") info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ = $(@sprintf("%.10f", float(eoi_SOS_L₁_dist)))")
distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist distance_to_cone = λ - 2^3*eoi_SOS_L₁_dist
return distance_to_cone return distance_to_cone
end end
function distance_to_cone{T<:Rational, S<:Interval}(κ::T, sqrt_matrix::Array{S,2}, Δ::GroupAlgebraElement{T}) function distance_to_cone{T<:Rational, S<:Interval}(λ::T, sqrt_matrix::Array{S,2}, Δ::GroupAlgebraElement{T})
SOS = compute_SOS(sqrt_matrix, Δ) SOS = compute_SOS(sqrt_matrix, Δ)
info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(GroupAlgebras.ɛ(SOS))") info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(GroupAlgebras.ɛ(SOS))")
κⁱⁿᵗ = @interval(κ) λⁱⁿᵗ = @interval(λ)
Δⁱⁿᵗ = GroupAlgebraElement([@interval(c) for c in Δ.coefficients], Δ.product_matrix) Δⁱⁿᵗ = GroupAlgebraElement([@interval(c) for c in Δ.coefficients], Δ.product_matrix)
SOS_diff = EOI(Δⁱⁿᵗ, κⁱⁿᵗ) - SOS SOS_diff = EOI(Δⁱⁿᵗ, λⁱⁿᵗ) - SOS
eoi_SOS_L₁_dist = norm(SOS_diff,1) eoi_SOS_L₁_dist = norm(SOS_diff,1)
info(logger, "κ = $κ (≈≥$(@sprintf("%.10f",float(κ))))") info(logger, "λ = $λ (≈≥$(@sprintf("%.10f",float(λ))))")
ɛ_dist = GroupAlgebras.ɛ(SOS_diff) ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)") info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L₁_dist)") info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L₁_dist)")
distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist distance_to_cone = λ - 2^3*eoi_SOS_L₁_dist
return distance_to_cone return distance_to_cone
end end
function distance_to_cone{T<:AbstractFloat}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T}) function distance_to_cone{T<:AbstractFloat}(λ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T})
SOS = compute_SOS(sqrt_matrix, Δ) SOS = compute_SOS(sqrt_matrix, Δ)
SOS_diff = EOI(Δ, κ) - SOS SOS_diff = EOI(Δ, λ) - SOS
eoi_SOS_L₁_dist = norm(SOS_diff,1) eoi_SOS_L₁_dist = norm(SOS_diff,1)
info(logger, "κ = $κ") info(logger, "λ = $λ")
ɛ_dist = GroupAlgebras.ɛ(SOS_diff) ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))") info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L₁_dist))") info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L₁_dist))")
distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist distance_to_cone = λ - 2^3*eoi_SOS_L₁_dist
return distance_to_cone return distance_to_cone
end end
function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; function check_distance_to_positive_cone(Δ::GroupAlgebraElement, λ, P;
tol=1e-7, rational=false) tol=1e-7, rational=false)
isapprox(eigvals(A), abs(eigvals(A)), atol=tol) || isapprox(eigvals(P), abs(eigvals(P)), atol=tol) ||
warn("The solution matrix doesn't seem to be positive definite!") warn("The solution matrix doesn't seem to be positive definite!")
@assert A == Symmetric(A) @assert P == Symmetric(P)
A_sqrt = real(sqrtm(A)) Q = real(sqrtm(P))
info(logger, "------------------------------------------------------------") info(logger, "------------------------------------------------------------")
info(logger, "") info(logger, "")
info(logger, "Checking in floating-point arithmetic...") info(logger, "Checking in floating-point arithmetic...")
t = @timed fp_distance = distance_to_cone(κ, A_sqrt, Δ) t = @timed fp_distance = distance_to_cone(λ, Q, Δ)
info(logger, timed_msg(t)) info(logger, timed_msg(t))
info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))") info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))")
info(logger, "------------------------------------------------------------") info(logger, "------------------------------------------------------------")
@ -133,17 +133,17 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A;
# return fpdistance # return fpdistance
# end # end
info(logger, "Projecting columns of rationalized A_sqrt to the augmentation ideal...") info(logger, "Projecting columns of rationalized Q to the augmentation ideal...")
δ = eps(κ) δ = eps(λ)
A_sqrt_ℚ = ℚ(A_sqrt, δ) Q_ℚ = ℚ(Q, δ)
t = @timed A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ) t = @timed Q_ℚω = correct_to_augmentation_ideal(Q_ℚ)
info(logger, timed_msg(t)) info(logger, timed_msg(t))
κ_ℚ = ℚ(κ, δ) λ_ℚ = ℚ(λ, δ)
Δ_ℚ = ℚ(Δ, δ) Δ_ℚ = ℚ(Δ, δ)
info(logger, "Checking in interval arithmetic") info(logger, "Checking in interval arithmetic")
A_sqrt_ℚ_augⁱⁿᵗ = Float64.(A_sqrt_ℚ_aug) ± δ Q_ℚωⁱⁿᵗ = Float64.(Q_ℚω) ± δ
t = @timed Interval_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_augⁱⁿᵗ, Δ_ℚ) t = @timed Interval_dist_to_Σ² = distance_to_cone(λ_ℚ, Q_ℚωⁱⁿᵗ, Δ_ℚ)
info(logger, timed_msg(t)) info(logger, timed_msg(t))
info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_Σ²)") info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_Σ²)")
info(logger, "------------------------------------------------------------") info(logger, "------------------------------------------------------------")
@ -152,7 +152,7 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A;
return Interval_dist_to_Σ² return Interval_dist_to_Σ²
else else
info(logger, "Checking Projected SOS decomposition in exact rational arithmetic...") info(logger, "Checking Projected SOS decomposition in exact rational arithmetic...")
t = @timed ℚ_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_aug, Δ_ℚ) t = @timed ℚ_dist_to_Σ² = distance_to_cone(λ_ℚ, Q_ℚω, Δ_ℚ)
info(logger, timed_msg(t)) info(logger, timed_msg(t))
@assert isa(ℚ_dist_to_Σ², Rational) @assert isa(ℚ_dist_to_Σ², Rational)
info(logger, "Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))") info(logger, "Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))")

20
src/sdps.jl
View File

@ -52,21 +52,21 @@ function create_SDP_problem(matrix_constraints, Δ::GroupAlgebraElement; upper_b
Δ² = Δ*Δ Δ² = Δ*Δ
@assert length(Δ) == length(matrix_constraints) @assert length(Δ) == length(matrix_constraints)
m = JuMP.Model(); m = JuMP.Model();
JuMP.@variable(m, A[1:N, 1:N], SDP) JuMP.@variable(m, P[1:N, 1:N], SDP)
JuMP.@SDconstraint(m, A >= 0) JuMP.@SDconstraint(m, P >= 0)
JuMP.@constraint(m, sum(A[i] for i in eachindex(A)) == 0) JuMP.@constraint(m, sum(P[i] for i in eachindex(P)) == 0)
if upper_bound < Inf if upper_bound < Inf
JuMP.@variable(m, 0.0 <= κ <= upper_bound) JuMP.@variable(m, 0.0 <= λ <= upper_bound)
else else
JuMP.@variable(m, κ >= 0) JuMP.@variable(m, λ >= 0)
end end
for (pairs, δ², δ) in zip(matrix_constraints, Δ².coefficients, Δ.coefficients) for (pairs, δ², δ) in zip(matrix_constraints, Δ².coefficients, Δ.coefficients)
JuMP.@constraint(m, sum(A[i,j] for (i,j) in pairs) == δ² - κ*δ) JuMP.@constraint(m, sum(P[i,j] for (i,j) in pairs) == δ² - λ*δ)
end end
JuMP.@objective(m, Max, κ) JuMP.@objective(m, Max, λ)
return m return m
end end
@ -95,7 +95,7 @@ function solve_SDP(SDP_problem, solver)
end end
info(logger, solution_status) info(logger, solution_status)
κ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :κ)) λ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :λ))
A = JuMP.getvalue(JuMP.getvariable(SDP_problem, :A)) P = JuMP.getvalue(JuMP.getvariable(SDP_problem, :P))
return κ, A return λ, P
end end