From a180025dc029a33d575f18a778e79da616aad93c Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 15:21:45 +0200 Subject: [PATCH 01/17] log time of projecting A onto augmentation ideal --- src/checksolution.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/checksolution.jl b/src/checksolution.jl index d56b53f..7cdd41e 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -129,7 +129,8 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; info(logger, "Projecting columns of rationalized A_sqrt to the augmentation ideal...") δ = eps(κ) A_sqrt_ℚ = ℚ(A_sqrt, δ) - A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ) + t = @timed A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ) + info(logger, timed_msg(t)) κ_ℚ = ℚ(κ, δ) Δ_ℚ = ℚ(Δ, δ) From f487c4635c55ab6a191d17fb242286fca5e120b6 Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 16:07:08 +0200 Subject: [PATCH 02/17] Cleanup of the check_propertyT function --- src/PropertyT.jl | 93 +++++++++++++++++++++++++++--------------------- 1 file changed, 52 insertions(+), 41 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 2426719..8f5b592 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -46,13 +46,21 @@ function ΔandSDPconstraints(name::String) return Δ, sdp_constraints end -function ΔandSDPconstraints(name::String, generating_set::Function) - pm_fname, Δ_fname = pmΔfilenames(name) - S, ID = generating_set() - Δ, sdp_constraints = Main.ΔandSDPconstraints(ID, S) - save(pm_fname, "pm", Δ.product_matrix) - save(Δ_fname, "Δ", Δ.coefficients) - return Δ, sdp_constraints +function ΔandSDPconstraints(name::String, generating_set::Function, radius::Int) + try + return ΔandSDPconstraints(name) + catch err + if isa(err, ArgumentError) + pm_fname, Δ_fname = pmΔfilenames(name) + S, Id = generating_set() + Δ, sdp_constraints = Main.ΔandSDPconstraints(Id, S, radius) + save(pm_fname, "pm", Δ.product_matrix) + save(Δ_fname, "Δ", Δ.coefficients) + return Δ, sdp_constraints + else + error(logger, err) + end + end end function κandA(name::String) @@ -78,14 +86,43 @@ function timed_msg(t) return "took: $elapsed s, allocated: $bytes_alloc bytes ($(gc_diff.poolalloc) allocations)." end -function κandA(name::String, sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver; upper_bound=Inf) - if isfile("$name/solver.log") - rm("$name/solver.log") + +function κandA(name::String, opts...) + try + return κandA(name) + catch err + if isa(err, ArgumentError) + if isfile("$name/solver.log") + rm("$name/solver.log") + end + + add_handler(solver_logger, DefaultHandler("./$name/solver.log", DefaultFormatter("{date}| {msg}")), "solver_log") + + info(logger, "Creating SDP problem...") + + return compute_κandA(opts...) + + remove_handler(solver_logger, "solver_log") + + κ_fname, A_fname = κSDPfilenames(name) + + if κ > 0 + save(κ_fname, "κ", κ) + save(A_fname, "A", A) + else + throw(ErrorException("Solver $solver did not produce a valid solution!: κ = $κ")) + end + return κ, A + + else + # throw(err) + error(logger, err) + end end +end - add_handler(solver_logger, DefaultHandler("./$name/solver.log", DefaultFormatter("{date}| {msg}")), "solver_log") +function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver; upper_bound=Inf) - info(logger, "Creating SDP problem...") t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound) info(logger, timed_msg(t)) @@ -98,16 +135,6 @@ function κandA(name::String, sdp_constraints, Δ::GroupAlgebraElement, solver:: warn(solver_logger, y) end end - - remove_handler(solver_logger, "solver_log") - - κ_fname, A_fname = κSDPfilenames(name) - if κ > 0 - save(κ_fname, "κ", κ) - save(A_fname, "A", A) - else - throw(ErrorException("Solver $solver did not produce a valid solution!: κ = $κ")) - end return κ, A end @@ -124,30 +151,14 @@ function check_property_T(name::String, generating_set::Function, info(logger, "Precision: $tol") info(logger, "Upper bound: $upper_bound") - Δ, sdp_constraints = try - ΔandSDPconstraints(name) - catch err - if isa(err, ArgumentError) - ΔandSDPconstraints(name, generating_set) - else - error(logger, err) - end - end + Δ, sdp_constraints = ΔandSDPconstraints(name, generating_set, radius) + S = countnz(Δ.coefficients) - 1 info(logger, "|S| = $S") info(logger, "length(Δ) = $(length(Δ))") info(logger, "size(Δ.product_matrix) = $(size(Δ.product_matrix))") - κ, A = try - κandA(name) - catch err - if isa(err, ArgumentError) - κandA(name, sdp_constraints, Δ, solver; upper_bound=upper_bound) - else - # throw(err) - error(logger, err) - end - end + κ, A = κandA(name, sdp_constraints, Δ, solver; upper_bound=upper_bound) info(logger, "κ = $κ") info(logger, "sum(A) = $(sum(A))") From 5c36eccc7bce595dff42dcdf9e0514fe8ae51c1b Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 16:12:00 +0200 Subject: [PATCH 03/17] cosmetics --- src/PropertyT.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 8f5b592..15898ca 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -169,9 +169,10 @@ function check_property_T(name::String, generating_set::Function, spectral_gap = check_distance_to_positive_cone(Δ, κ, A, tol=tol, rational=false) if spectral_gap > 0 Kazhdan_κ = sqrt(2*spectral_gap/S) - Kazhdan_κ = Float64(trunc(Kazhdan_κ,12)) + Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12)) info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!") else + spectral_gap = Float64(trunc(spectral_gap, 12)) info(logger, "λ($name, S) ≥ $spectral_gap: Group may NOT HAVE property (T)!") end else From e08e787c1a34c5c5a7fad16ab29b8ea57f1869ed Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:32:36 +0200 Subject: [PATCH 04/17] incorporate radius parameter --- src/PropertyT.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 15898ca..a5998a1 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -121,7 +121,7 @@ function κandA(name::String, opts...) end end -function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver; upper_bound=Inf) +function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver, upper_bound=Inf) t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound) info(logger, timed_msg(t)) @@ -139,7 +139,7 @@ function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::Abstra end function check_property_T(name::String, generating_set::Function, - solver, upper_bound, tol=1e-6) + solver, upper_bound, tol, radius) if !isdir(name) mkdir(name) @@ -158,7 +158,7 @@ function check_property_T(name::String, generating_set::Function, info(logger, "length(Δ) = $(length(Δ))") info(logger, "size(Δ.product_matrix) = $(size(Δ.product_matrix))") - κ, A = κandA(name, sdp_constraints, Δ, solver; upper_bound=upper_bound) + κ, A = κandA(name, sdp_constraints, Δ, solver, upper_bound) info(logger, "κ = $κ") info(logger, "sum(A) = $(sum(A))") From ada830f5d551fa23ce7ff3b2aa13d964eef519b6 Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:33:17 +0200 Subject: [PATCH 05/17] fix: don't return after solving, save kappa and A first --- src/PropertyT.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index a5998a1..702dab1 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -100,7 +100,7 @@ function κandA(name::String, opts...) info(logger, "Creating SDP problem...") - return compute_κandA(opts...) + κ, A = compute_κandA(opts...) remove_handler(solver_logger, "solver_log") From 837988c381e5c7fbc4571e9bd95acd0968a444c9 Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:34:28 +0200 Subject: [PATCH 06/17] fix: don't throw exception in kandA k <=0 is handled in the outside loop --- src/PropertyT.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 702dab1..9585359 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -109,8 +109,6 @@ function κandA(name::String, opts...) if κ > 0 save(κ_fname, "κ", κ) save(A_fname, "A", A) - else - throw(ErrorException("Solver $solver did not produce a valid solution!: κ = $κ")) end return κ, A From 69e74190af1efcbd1b2d2c66df31dd2a3f0b8106 Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:34:56 +0200 Subject: [PATCH 07/17] tad better try...catch syntax --- src/PropertyT.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 9585359..09d091a 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -124,13 +124,11 @@ function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::Abstra t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound) info(logger, timed_msg(t)) - κ = 0.0 - A = nothing while κ == 0.0 - try - κ, A = solve_SDP(SDP_problem, solver) + κ, A = try + solve_SDP(SDP_problem, solver) catch y - warn(solver_logger, y) + warn(logger, y) end end return κ, A From 5aa16314ec50f5b8031bddc778be87195cde5744 Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:35:30 +0200 Subject: [PATCH 08/17] do not force @bigintervals when we use Floats64 in a moment --- src/checksolution.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/checksolution.jl b/src/checksolution.jl index 7cdd41e..ef8680e 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -46,7 +46,7 @@ end import ValidatedNumerics.± function (±){T<:Number}(X::AbstractArray{T}, tol::Real) - r{T}(x::T) = (x == zero(T)? @biginterval(0) : x ± tol) + r{T}(x::T) = (x == zero(T)? @interval(0) : x ± tol) return r.(X) end From 98cdeedc7aad62eaae5289c8d06c2b821eb9a5e3 Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:36:18 +0200 Subject: [PATCH 09/17] make A_sqrt_Q_aug \pm \delta a matrix of Float64 intervals --- src/checksolution.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/checksolution.jl b/src/checksolution.jl index ef8680e..3577e96 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -135,8 +135,8 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; Δ_ℚ = ℚ(Δ, δ) info(logger, "Checking in interval arithmetic") - A_sqrt_ℚ_augᴵ = A_sqrt_ℚ_aug ± δ - t = @timed Interval_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_augᴵ, Δ_ℚ) + A_sqrt_ℚ_augⁱⁿᵗ = Float64.(A_sqrt_ℚ_aug) ± δ + t = @timed Interval_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_augⁱⁿᵗ, Δ_ℚ) info(logger, timed_msg(t)) info(logger, "The Augmentation-projected actual distance (to positive cone) ≥ $(@sprintf("%.10f", Interval_dist_to_Σ².lo))") info(logger, "------------------------------------------------------------") From 2b1bb7b765c94d41c5ea5b107296ec93f458ae3c Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:38:18 +0200 Subject: [PATCH 10/17] break after checking in floats return negatives (commented out) --- src/checksolution.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/checksolution.jl b/src/checksolution.jl index 3577e96..00fecf2 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -126,6 +126,10 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))") info(logger, "------------------------------------------------------------") + # if fp_distance ≤ 0 + # return fpdistance + # end + info(logger, "Projecting columns of rationalized A_sqrt to the augmentation ideal...") δ = eps(κ) A_sqrt_ℚ = ℚ(A_sqrt, δ) From 450d43ed14e11862d6fe705c717f47e10a5cf087 Mon Sep 17 00:00:00 2001 From: kalmar Date: Fri, 31 Mar 2017 22:38:51 +0200 Subject: [PATCH 11/17] ValidatedNumerics parameters -- to be checked --- src/checksolution.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/checksolution.jl b/src/checksolution.jl index 00fecf2..abf50aa 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -2,7 +2,9 @@ import Base: rationalize using ValidatedNumerics ValidatedNumerics.setrounding(Interval, :correct) +# ValidatedNumerics.setrounding(Interval, :fast) #which is slower?? ValidatedNumerics.setformat(:standard) +# setprecision(Interval, 53) # slightly faster than 256 function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, κ::T) return Δ*Δ - κ*Δ From a827f8c42528c066ed295bf6597ff9feb2faaf74 Mon Sep 17 00:00:00 2001 From: kalmar Date: Sat, 1 Apr 2017 08:27:54 +0200 Subject: [PATCH 12/17] Revert "fix: don't throw exception in kandA" This reverts commit 837988c381e5c7fbc4571e9bd95acd0968a444c9. --- src/PropertyT.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 09d091a..3b8f4ab 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -109,6 +109,8 @@ function κandA(name::String, opts...) if κ > 0 save(κ_fname, "κ", κ) save(A_fname, "A", A) + else + throw(ErrorException("Solver $solver did not produce a valid solution!: κ = $κ")) end return κ, A From b20cd02d5bc50f3fc8de45d85e3dbf148f7bbba9 Mon Sep 17 00:00:00 2001 From: kalmar Date: Sat, 1 Apr 2017 08:32:34 +0200 Subject: [PATCH 13/17] Revert "tad better try...catch syntax" This reverts commit 69e74190af1efcbd1b2d2c66df31dd2a3f0b8106. --- src/PropertyT.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 3b8f4ab..702dab1 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -126,11 +126,13 @@ function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::Abstra t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound) info(logger, timed_msg(t)) + κ = 0.0 + A = nothing while κ == 0.0 - κ, A = try - solve_SDP(SDP_problem, solver) + try + κ, A = solve_SDP(SDP_problem, solver) catch y - warn(logger, y) + warn(solver_logger, y) end end return κ, A From dbed75103eb6c2243ef2c277f91aeb18241851c4 Mon Sep 17 00:00:00 2001 From: kalmar Date: Sat, 1 Apr 2017 14:21:01 +0200 Subject: [PATCH 14/17] allow check_distance_to_positive_cone return Interval --- src/PropertyT.jl | 3 +++ src/checksolution.jl | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 702dab1..e9b8b30 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -167,6 +167,9 @@ function check_property_T(name::String, generating_set::Function, if κ > 0 spectral_gap = check_distance_to_positive_cone(Δ, κ, A, tol=tol, rational=false) + if isa(spectral_gap, Interval) + spectral_gap = spectral_gap.lo + end if spectral_gap > 0 Kazhdan_κ = sqrt(2*spectral_gap/S) Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12)) diff --git a/src/checksolution.jl b/src/checksolution.jl index abf50aa..4d0a9cc 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -144,11 +144,11 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; A_sqrt_ℚ_augⁱⁿᵗ = Float64.(A_sqrt_ℚ_aug) ± δ t = @timed Interval_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_augⁱⁿᵗ, Δ_ℚ) info(logger, timed_msg(t)) - info(logger, "The Augmentation-projected actual distance (to positive cone) ≥ $(@sprintf("%.10f", Interval_dist_to_Σ².lo))") + info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_Σ²)") info(logger, "------------------------------------------------------------") if Interval_dist_to_Σ².lo ≤ 0 || !rational - return Interval_dist_to_Σ².lo + return Interval_dist_to_Σ² else info(logger, "Checking Projected SOS decomposition in exact rational arithmetic...") t = @timed ℚ_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_aug, Δ_ℚ) From 41e53a598a9a8fc0dade43f19f5a54e1201410a2 Mon Sep 17 00:00:00 2001 From: kalmar Date: Sat, 1 Apr 2017 14:22:01 +0200 Subject: [PATCH 15/17] Finally move all intense interval arithmetic to Interval{Float64} --- src/checksolution.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/checksolution.jl b/src/checksolution.jl index 4d0a9cc..a1f40d7 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -83,11 +83,12 @@ end function distance_to_cone{T<:Rational, S<:Interval}(κ::T, sqrt_matrix::Array{S,2}, Δ::GroupAlgebraElement{T}) SOS = compute_SOS(sqrt_matrix, Δ) info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(GroupAlgebras.ɛ(SOS))") - - SOS_diff = EOI(Δ, κ) - SOS + κⁱⁿᵗ = @interval(κ) + Δⁱⁿᵗ = GroupAlgebraElement([@interval(c) for c in Δ.coefficients], Δ.product_matrix) + SOS_diff = EOI(Δⁱⁿᵗ, κⁱⁿᵗ) - SOS eoi_SOS_L₁_dist = norm(SOS_diff,1) - info(logger, "κ = $κ (≈$(@sprintf("%.10f",float(κ))))") + info(logger, "κ = $κ (≈≥$(@sprintf("%.10f",float(κ))))") ɛ_dist = GroupAlgebras.ɛ(SOS_diff) info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)") From 028979bfd14a8ef487b7a64ef11ec2da523e0a09 Mon Sep 17 00:00:00 2001 From: kalmar Date: Sat, 1 Apr 2017 14:22:30 +0200 Subject: [PATCH 16/17] better printing --- src/PropertyT.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index e9b8b30..bf325c4 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -138,6 +138,8 @@ function compute_κandA(sdp_constraints, Δ::GroupAlgebraElement, solver::Abstra return κ, A end +Kazhdan_from_sgap(λ,N) = sqrt(2*λ/N) + function check_property_T(name::String, generating_set::Function, solver, upper_bound, tol, radius) @@ -171,9 +173,10 @@ function check_property_T(name::String, generating_set::Function, spectral_gap = spectral_gap.lo end if spectral_gap > 0 - Kazhdan_κ = sqrt(2*spectral_gap/S) - Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12)) - info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!") + @show spectral_gap + Kazhdan_const = Kazhdan_from_sgap(spectral_gap, S) + Kazhdan_const = Float64(trunc(Kazhdan_const, 12)) + info(logger, "κ($name, S) ≥ $Kazhdan_const: Group HAS property (T)!") else spectral_gap = Float64(trunc(spectral_gap, 12)) info(logger, "λ($name, S) ≥ $spectral_gap: Group may NOT HAVE property (T)!") From 5876998cbaf2ed807de3b76b5c9abb531708b685 Mon Sep 17 00:00:00 2001 From: kalmar Date: Sat, 1 Apr 2017 15:21:57 +0200 Subject: [PATCH 17/17] Change names: kappa -> lambda, A -> P --- src/PropertyT.jl | 83 ++++++++++++++++++++++---------------------- src/checksolution.jl | 68 ++++++++++++++++++------------------ src/sdps.jl | 20 +++++------ 3 files changed, 86 insertions(+), 85 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index bf325c4..9b62c7d 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -20,14 +20,14 @@ function pmΔfilenames(name::String) return pm_filename, Δ_coeff_filename end -function κSDPfilenames(name::String) +function λSDPfilenames(name::String) if !isdir(name) mkdir(name) end prefix = name - κ_filename = joinpath(prefix, "kappa.jld") - SDP_filename = joinpath(prefix, "SDPmatrixA.jld") - return κ_filename, SDP_filename + λ_filename = joinpath(prefix, "lambda.jld") + SDP_filename = joinpath(prefix, "SDPmatrix.jld") + return λ_filename, SDP_filename end function ΔandSDPconstraints(name::String) @@ -63,18 +63,19 @@ function ΔandSDPconstraints(name::String, generating_set::Function, radius::Int end end -function κandA(name::String) - κ_fname, SDP_fname = κSDPfilenames(name) - f₁ = isfile(κ_fname) +function λandP(name::String) + λ_fname, SDP_fname = λSDPfilenames(name) + f₁ = isfile(λ_fname) f₂ = isfile(SDP_fname) + if f₁ && f₂ - info(logger, "Loading precomputed κ, A...") - κ = load(κ_fname, "κ") - A = load(SDP_fname, "A") + info(logger, "Loading precomputed λ, P...") + λ = load(λ_fname, "λ") + P = load(SDP_fname, "P") 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 - return κ, A + return λ, P end function timed_msg(t) @@ -87,32 +88,32 @@ function timed_msg(t) end -function κandA(name::String, opts...) +function λandP(name::String, opts...) try - return κandA(name) + return λandP(name) catch err if isa(err, ArgumentError) - if isfile("$name/solver.log") - rm("$name/solver.log") + if isfile(joinpath(name, "solver.log")) + rm(joinpath(name, "solver.log")) 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...") - κ, A = compute_κandA(opts...) + λ, P = compute_λandP(opts...) remove_handler(solver_logger, "solver_log") - κ_fname, A_fname = κSDPfilenames(name) + λ_fname, P_fname = λSDPfilenames(name) - if κ > 0 - save(κ_fname, "κ", κ) - save(A_fname, "A", A) + if λ > 0 + save(λ_fname, "λ", λ) + save(P_fname, "P", P) else - throw(ErrorException("Solver $solver did not produce a valid solution!: κ = $κ")) + throw(ErrorException("Solver $solver did not produce a valid solution!: λ = $λ")) end - return κ, A + return λ, P else # throw(err) @@ -121,21 +122,21 @@ function κandA(name::String, opts...) 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) info(logger, timed_msg(t)) - κ = 0.0 - A = nothing - while κ == 0.0 + λ = 0.0 + P = nothing + while λ == 0.0 try - κ, A = solve_SDP(SDP_problem, solver) + λ, P = solve_SDP(SDP_problem, solver) catch y warn(solver_logger, y) end end - return κ, A + return λ, P end 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, "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, "sum(A) = $(sum(A))") - info(logger, "maximum(A) = $(maximum(A))") - info(logger, "minimum(A) = $(minimum(A))") + info(logger, "λ = $λ") + info(logger, "sum(P) = $(sum(P))") + info(logger, "maximum(P) = $(maximum(P))") + info(logger, "minimum(P) = $(minimum(P))") - if κ > 0 - spectral_gap = check_distance_to_positive_cone(Δ, κ, A, tol=tol, rational=false) + if λ > 0 + spectral_gap = check_distance_to_positive_cone(Δ, λ, P, tol=tol, rational=false) if isa(spectral_gap, Interval) spectral_gap = spectral_gap.lo end if spectral_gap > 0 @show spectral_gap - Kazhdan_const = Kazhdan_from_sgap(spectral_gap, S) - Kazhdan_const = Float64(trunc(Kazhdan_const, 12)) - info(logger, "κ($name, S) ≥ $Kazhdan_const: Group HAS property (T)!") + Kazhdan_κ = Kazhdan_from_sgap(spectral_gap, S) + Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12)) + info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!") else spectral_gap = Float64(trunc(spectral_gap, 12)) info(logger, "λ($name, S) ≥ $spectral_gap: Group may NOT HAVE property (T)!") end 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 diff --git a/src/checksolution.jl b/src/checksolution.jl index a1f40d7..8d72f8e 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -6,8 +6,8 @@ ValidatedNumerics.setrounding(Interval, :correct) ValidatedNumerics.setformat(:standard) # setprecision(Interval, 53) # slightly faster than 256 -function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, κ::T) - return Δ*Δ - κ*Δ +function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, λ::T) + return Δ*Δ - λ*Δ end function algebra_square(vector, elt) @@ -62,69 +62,69 @@ end ℚ(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_diff = EOI(Δ, κ) - SOS + SOS_diff = EOI(Δ, λ) - SOS eoi_SOS_L₁_dist = norm(SOS_diff,1) - info(logger, "κ = $κ (≈$(@sprintf("%.10f", float(κ)))") + info(logger, "λ = $λ (≈$(@sprintf("%.10f", float(λ)))") ɛ_dist = GroupAlgebras.ɛ(SOS_diff) if ɛ_dist ≠ 0//1 warn(logger, "The SOS is not in the augmentation ideal, number below are meaningless!") end - info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) = $ɛ_dist") - info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ = $(@sprintf("%.10f", float(eoi_SOS_L₁_dist)))") + info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) = $ɛ_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 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, Δ) info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(GroupAlgebras.ɛ(SOS))") - κⁱⁿᵗ = @interval(κ) + λⁱⁿᵗ = @interval(λ) Δⁱⁿᵗ = 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) - info(logger, "κ = $κ (≈≥$(@sprintf("%.10f",float(κ))))") + info(logger, "λ = $λ (≈≥$(@sprintf("%.10f",float(λ))))") ɛ_dist = GroupAlgebras.ɛ(SOS_diff) - info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)") - info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L₁_dist)") + info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_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 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_diff = EOI(Δ, κ) - SOS + SOS_diff = EOI(Δ, λ) - SOS eoi_SOS_L₁_dist = norm(SOS_diff,1) - info(logger, "κ = $κ") + info(logger, "λ = $λ") ɛ_dist = GroupAlgebras.ɛ(SOS_diff) - info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))") - info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L₁_dist))") + info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_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 end -function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; +function check_distance_to_positive_cone(Δ::GroupAlgebraElement, λ, P; 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!") - @assert A == Symmetric(A) - A_sqrt = real(sqrtm(A)) + @assert P == Symmetric(P) + Q = real(sqrtm(P)) info(logger, "------------------------------------------------------------") info(logger, "") 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, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))") info(logger, "------------------------------------------------------------") @@ -133,17 +133,17 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; # return fpdistance # end - info(logger, "Projecting columns of rationalized A_sqrt to the augmentation ideal...") - δ = eps(κ) - A_sqrt_ℚ = ℚ(A_sqrt, δ) - t = @timed A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ) + info(logger, "Projecting columns of rationalized Q to the augmentation ideal...") + δ = eps(λ) + Q_ℚ = ℚ(Q, δ) + t = @timed Q_ℚω = correct_to_augmentation_ideal(Q_ℚ) info(logger, timed_msg(t)) - κ_ℚ = ℚ(κ, δ) + λ_ℚ = ℚ(λ, δ) Δ_ℚ = ℚ(Δ, δ) info(logger, "Checking in interval arithmetic") - A_sqrt_ℚ_augⁱⁿᵗ = Float64.(A_sqrt_ℚ_aug) ± δ - t = @timed Interval_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_augⁱⁿᵗ, Δ_ℚ) + Q_ℚωⁱⁿᵗ = Float64.(Q_ℚω) ± δ + t = @timed Interval_dist_to_Σ² = distance_to_cone(λ_ℚ, Q_ℚωⁱⁿᵗ, Δ_ℚ) info(logger, timed_msg(t)) info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_Σ²)") info(logger, "------------------------------------------------------------") @@ -152,7 +152,7 @@ function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; return Interval_dist_to_Σ² else 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)) @assert isa(ℚ_dist_to_Σ², Rational) info(logger, "Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))") diff --git a/src/sdps.jl b/src/sdps.jl index 4eeb980..b8199c6 100644 --- a/src/sdps.jl +++ b/src/sdps.jl @@ -52,21 +52,21 @@ function create_SDP_problem(matrix_constraints, Δ::GroupAlgebraElement; upper_b Δ² = Δ*Δ @assert length(Δ) == length(matrix_constraints) m = JuMP.Model(); - JuMP.@variable(m, A[1:N, 1:N], SDP) - JuMP.@SDconstraint(m, A >= 0) - JuMP.@constraint(m, sum(A[i] for i in eachindex(A)) == 0) + JuMP.@variable(m, P[1:N, 1:N], SDP) + JuMP.@SDconstraint(m, P >= 0) + JuMP.@constraint(m, sum(P[i] for i in eachindex(P)) == 0) if upper_bound < Inf - JuMP.@variable(m, 0.0 <= κ <= upper_bound) + JuMP.@variable(m, 0.0 <= λ <= upper_bound) else - JuMP.@variable(m, κ >= 0) + JuMP.@variable(m, λ >= 0) end 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 - JuMP.@objective(m, Max, κ) + JuMP.@objective(m, Max, λ) return m end @@ -95,7 +95,7 @@ function solve_SDP(SDP_problem, solver) end info(logger, solution_status) - κ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :κ)) - A = JuMP.getvalue(JuMP.getvariable(SDP_problem, :A)) - return κ, A + λ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :λ)) + P = JuMP.getvalue(JuMP.getvariable(SDP_problem, :P)) + return λ, P end