From d6435c1d4486815e9b5b4acc095d6ccf4143554a Mon Sep 17 00:00:00 2001 From: kalmar Date: Mon, 13 Mar 2017 14:49:55 +0100 Subject: [PATCH] Bootstrap Julia Package --- .codecov.yml | 1 + .gitignore | 7 +- .travis.yml | 19 +++ README.md | 1 + REQUIRE | 1 + appveyor.yml | 34 ++++ property(T).jl | 366 ------------------------------------------- src/GroupAlgebras.jl | 133 ++++++++++++++++ src/Property(T).jl | 135 ++++++++++++++++ src/checksolution.jl | 142 +++++++++++++++++ src/sdps.jl | 87 ++++++++++ test/runtests.jl | 5 + 12 files changed, 561 insertions(+), 370 deletions(-) create mode 100644 .codecov.yml create mode 100644 .travis.yml create mode 100644 README.md create mode 100644 REQUIRE create mode 100644 appveyor.yml delete mode 100644 property(T).jl create mode 100644 src/GroupAlgebras.jl create mode 100644 src/Property(T).jl create mode 100644 src/checksolution.jl create mode 100644 src/sdps.jl create mode 100644 test/runtests.jl diff --git a/.codecov.yml b/.codecov.yml new file mode 100644 index 0000000..69cb760 --- /dev/null +++ b/.codecov.yml @@ -0,0 +1 @@ +comment: false diff --git a/.gitignore b/.gitignore index 03781b8..8c960ec 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,3 @@ -.*~ -.ipynb_checkpoints -*.ipynb -*.gws +*.jl.cov +*.jl.*.cov +*.jl.mem diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..69cdfee --- /dev/null +++ b/.travis.yml @@ -0,0 +1,19 @@ +# Documentation: http://docs.travis-ci.com/user/languages/julia/ +language: julia +os: + - linux + - osx +julia: + - release + - nightly +notifications: + email: false +# uncomment the following lines to override the default test script +#script: +# - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi +# - julia -e 'Pkg.clone(pwd()); Pkg.build("Property(T)"); Pkg.test("Property(T)"; coverage=true)' +after_success: + # push coverage results to Coveralls + - julia -e 'cd(Pkg.dir("Property(T)")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())' + # push coverage results to Codecov + - julia -e 'cd(Pkg.dir("Property(T)")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())' diff --git a/README.md b/README.md new file mode 100644 index 0000000..6f370e7 --- /dev/null +++ b/README.md @@ -0,0 +1 @@ +# Property(T) diff --git a/REQUIRE b/REQUIRE new file mode 100644 index 0000000..94237c0 --- /dev/null +++ b/REQUIRE @@ -0,0 +1 @@ +julia 0.5 diff --git a/appveyor.yml b/appveyor.yml new file mode 100644 index 0000000..475c68c --- /dev/null +++ b/appveyor.yml @@ -0,0 +1,34 @@ +environment: + matrix: + - JULIAVERSION: "julialang/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe" + - JULIAVERSION: "julialang/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe" + - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe" + - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe" + +branches: + only: + - master + - /release-.*/ + +notifications: + - provider: Email + on_build_success: false + on_build_failure: false + on_build_status_changed: false + +install: +# Download most recent Julia Windows binary + - ps: (new-object net.webclient).DownloadFile( + $("http://s3.amazonaws.com/"+$env:JULIAVERSION), + "C:\projects\julia-binary.exe") +# Run installer silently, output to C:\projects\julia + - C:\projects\julia-binary.exe /S /D=C:\projects\julia + +build_script: +# Need to convert from shallow to complete for Pkg.clone to work + - IF EXIST .git\shallow (git fetch --unshallow) + - C:\projects\julia\bin\julia -e "versioninfo(); + Pkg.clone(pwd(), \"Property(T)\"); Pkg.build(\"Property(T)\")" + +test_script: + - C:\projects\julia\bin\julia -e "Pkg.test(\"Property(T)\")" diff --git a/property(T).jl b/property(T).jl deleted file mode 100644 index bf9c5ab..0000000 --- a/property(T).jl +++ /dev/null @@ -1,366 +0,0 @@ -using JuMP -import MathProgBase: AbstractMathProgSolver -import Base: rationalize -using GroupAlgebras - -using ProgressMeter -using ValidatedNumerics - -function create_product_matrix(basis, limit) - product_matrix = zeros(Int, (limit,limit)) - basis_dict = Dict{Array, Int}(x => i - for (i,x) in enumerate(basis)) - for i in 1:limit - x_inv::eltype(basis) = inv(basis[i]) - for j in 1:limit - w = x_inv*basis[j] - product_matrix[i,j] = basis_dict[w] - # index = findfirst(basis, w) - # index ≠ 0 || throw(ArgumentError("Product is not supported on basis: $w")) - # product_matrix[i,j] = index - end - end - return product_matrix -end - -function constraints_from_pm(pm, total_length=maximum(pm)) - n = size(pm,1) - constraints = constraints = [Array{Int,1}[] for x in 1:total_length] - for j in 1:n - Threads.@threads for i in 1:n - idx = pm[i,j] - push!(constraints[idx], [i,j]) - end - end - return constraints -end - -function splaplacian_coeff(S, basis, n=length(basis)) - result = spzeros(n) - result[1] = float(length(S)) - for s in S - ind = findfirst(basis, s) - result[ind] += -1.0 - end - return result -end - -function laplacian_coeff(S, basis) - return full(splaplacian_coeff(S,basis)) -end - - -function create_SDP_problem(matrix_constraints, Δ::GroupAlgebraElement; upper_bound=Inf) - N = size(Δ.product_matrix,1) - const Δ² = Δ*Δ - @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, κ >= 0.0) - if upper_bound < Inf - JuMP.@constraint(m, κ <= upper_bound) - end - JuMP.@objective(m, Max, κ) - - for (pairs, δ², δ) in zip(matrix_constraints, Δ².coefficients, Δ.coefficients) - JuMP.@constraint(m, sum(A[i,j] for (i,j) in pairs) == δ² - κ*δ) - end - return m -end - -function solve_SDP(SDP_problem, solver) - @show SDP_problem - @show solver - - JuMP.setsolver(SDP_problem, solver); - # @time MathProgBase.writeproblem(SDP_problem, "/tmp/SDP_problem") - solution_status = JuMP.solve(SDP_problem); - - if solution_status != :Optimal - warn("The solver did not solve the problem successfully!") - end - @show solution_status - - κ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :κ)) - A = JuMP.getvalue(JuMP.getvariable(SDP_problem, :A)) - return κ, A -end - -function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, κ::T) - return Δ*Δ - κ*Δ -end - -function square_as_elt(vector, elt) - zzz = zeros(elt.coefficients) - zzz[1:length(vector)] = vector -# new_base_elt = GroupAlgebraElement(zzz, elt.product_matrix) -# return (new_base_elt*new_base_elt).coefficients - return GroupAlgebras.algebra_multiplication(zzz, zzz, elt.product_matrix) -end - -function compute_SOS{T<:Number}(sqrt_matrix::Array{T,2}, - elt::GroupAlgebraElement{T}) - n = size(sqrt_matrix,2) - result = zeros(T, length(elt.coefficients)) - p = Progress(n, 1, "Checking SOS decomposition...", 50) - for i in 1:n - result .+= square_as_elt(sqrt_matrix[:,i], elt) - next!(p) - end - return GroupAlgebraElement{T}(result, elt.product_matrix) -end - -function correct_to_augmentation_ideal{T<:Rational}(sqrt_matrix::Array{T,2}) - sqrt_corrected = similar(sqrt_matrix) - l = size(sqrt_matrix,2) - for i in 1:l - col = view(sqrt_matrix,:,i) - sqrt_corrected[:,i] = col - sum(col)//l - # @assert sum(sqrt_corrected[:,i]) == 0 - end - return sqrt_corrected -end - -function check_solution{T<:Number}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T}; verbose=true, augmented=false) - result = compute_SOS(sqrt_matrix, Δ) - if augmented - epsilon = GroupAlgebras.ɛ(result) - if isa(epsilon, Interval) - @assert 0 in epsilon - elseif isa(epsilon, Rational) - @assert epsilon == 0//1 - else - warn("Does checking for augmentation has meaning for $(typeof(epsilon))?") - end - end - SOS_diff = EOI(Δ, κ) - result - - eoi_SOS_L₁_dist = norm(SOS_diff,1) - - if verbose - @show κ - if augmented - println("ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) = ", GroupAlgebras.ɛ(SOS_diff)) - else - ɛ_dist = GroupAlgebras.ɛ(SOS_diff) - if typeof(ɛ_dist) <: Interval - ɛ_dist = ɛ_dist.lo - end - @printf("ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ≈ %.10f\n", ɛ_dist) - end - - L₁_dist = eoi_SOS_L₁_dist - if typeof(L₁_dist) <: Interval - L₁_dist = L₁_dist.lo - end - @printf("‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ≈ %.10f\n", L₁_dist) - end - - distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist - return distance_to_cone -end - -import ValidatedNumerics.± -function (±)(X::AbstractArray, tol::Real) - r{T}(x::T) = ( x==zero(T) ? @interval(x) : x ± tol) - return r.(X) -end - -(±)(X::GroupAlgebraElement, tol::Real) = GroupAlgebraElement(X.coefficients ± tol, X.product_matrix) - -function Base.rationalize{T<:Integer, S<:Real}(::Type{T}, - X::AbstractArray{S}; tol::Real=eps(eltype(X))) - r(x) = rationalize(T, x, tol=tol) - return r.(X) -end - -ℚ(x, tol::Real) = rationalize(BigInt, x, tol=tol) - - -function ℚ_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; - tol=1e-7, verbose=true, rational=false) - - isapprox(eigvals(A), abs(eigvals(A)), atol=tol) || - warn("The solution matrix doesn't seem to be positive definite!") - @assert A == Symmetric(A) - A_sqrt = real(sqrtm(A)) - - # println("") - # println("Checking in floating-point arithmetic...") - # @time fp_distance = check_solution(κ, A_sqrt, Δ, verbose=verbose) - # println("Floating point distance (to positive cone) ≈ $(Float64(trunc(fp_distance,8)))") - # println("-------------------------------------------------------------") - # println("") - # - # if fp_distance ≤ 0 - # return fp_distance - # end - - println("Checking in interval arithmetic...") - A_sqrtᴵ = A_sqrt ± tol - κᴵ = κ ± tol - Δᴵ = Δ ± tol - @time Interval_distance = check_solution(κᴵ, A_sqrtᴵ, Δᴵ, verbose=verbose) - # @assert isa(ℚ_distance, Rational) - println("The actual distance (to positive cone) is contained in $Interval_distance") - println("-------------------------------------------------------------") - println("") - - if Interval_distance.lo ≤ 0 - return Interval_distance.lo - end - - println("Projecting columns of A_sqrt to the augmentation ideal...") - A_sqrt_ℚ = ℚ(A_sqrt, tol) - A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ) - κ_ℚ = ℚ(κ, tol) - Δ_ℚ = ℚ(Δ, tol) - - A_sqrt_ℚ_augᴵ = A_sqrt_ℚ_aug ± tol - κᴵ = κ_ℚ ± tol - Δᴵ = Δ_ℚ ± tol - @time Interval_dist_to_Σ² = check_solution(κᴵ, A_sqrt_ℚ_augᴵ, Δᴵ, verbose=verbose, augmented=true) - println("The Augmentation-projected actual distance (to positive cone) is contained in $Interval_dist_to_Σ²") - println("-------------------------------------------------------------") - println("") - - if Interval_dist_to_Σ².lo ≤ 0 || !rational - return Interval_dist_to_Σ².lo - else - - println("Checking Projected SOS decomposition in exact rational arithmetic...") - @time ℚ_dist_to_Σ² = check_solution(κ_ℚ, A_sqrt_ℚ_aug, Δ_ℚ, verbose=verbose, augmented=true) - @assert isa(ℚ_dist_to_Σ², Rational) - println("Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))") - println("-------------------------------------------------------------") - return ℚ_dist_to_Σ² - end -end - -function pmΔfilenames(name::String) - if !isdir(name) - mkdir(name) - end - prefix = name - pm_filename = joinpath(prefix, "product_matrix.jld") - Δ_coeff_filename = joinpath(prefix, "delta.coefficients.jld") - return pm_filename, Δ_coeff_filename -end - -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 -end - -function ΔandSDPconstraints(name::String) - pm_fname, Δ_fname = pmΔfilenames(name) - f₁ = isfile(pm_fname) - f₂ = isfile(Δ_fname) - if f₁ && f₂ - println("Loading precomputed pm, Δ, sdp_constraints...") - product_matrix = load(pm_fname, "pm") - L = load(Δ_fname, "Δ")[:, 1] - Δ = GroupAlgebraElement(L, Array{Int,2}(product_matrix)) - sdp_constraints = constraints_from_pm(product_matrix) - else - throw(ArgumentError("You need to precompute pm and Δ to load it!")) - end - return Δ, sdp_constraints -end - -function ΔandSDPconstraints(name::String, ID, generating_func::Function) - pm_fname, Δ_fname = pmΔfilenames(name) - Δ, sdp_constraints = ΔandSDPconstraints(ID, generating_func()) - save(pm_fname, "pm", Δ.product_matrix) - save(Δ_fname, "Δ", Δ.coefficients) - return Δ, sdp_constraints -end - -function κandA(name::String) - κ_fname, SDP_fname = κSDPfilenames(name) - f₁ = isfile(κ_fname) - f₂ = isfile(SDP_fname) - if f₁ && f₂ - println("Loading precomputed κ, A...") - κ = load(κ_fname, "κ") - A = load(SDP_fname, "A") - else - throw(ArgumentError("You need to precompute κ and SDP matrix A to load it!")) - end - return κ, A -end - -function κandA(name::String, sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver; upper_bound=Inf) - println("Creating SDP problem...") - @time SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound) - println("Solving SDP problem maximizing κ...") - κ, A = solve_SDP(SDP_problem, solver) - κ_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 - -function check_property_T(name::String, ID, generate_B₄::Function; - verbose=true, tol=1e-6, upper_bound=Inf) - - # solver = MosekSolver(INTPNT_CO_TOL_REL_GAP=tol, QUIET=!verbose) - solver = SCSSolver(eps=tol, max_iters=100000, verbose=verbose) - - @show name - @show verbose - @show tol - - - Δ, sdp_constraints = try - ΔandSDPconstraints(name) - catch err - if isa(err, ArgumentError) - ΔandSDPconstraints(name, ID, generate_B₄) - else - throw(err) - end - end - println("|S| = $(countnz(Δ.coefficients) -1)") - @show length(Δ) - @show 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) - end - end - - @show κ - @show sum(A) - @show maximum(A) - @show minimum(A) - - if κ > 0 - - true_kappa = ℚ_distance_to_positive_cone(Δ, κ, A, tol=tol, verbose=verbose, rational=true) - true_kappa = Float64(trunc(true_kappa,12)) - if true_kappa > 0 - println("κ($name, S) ≥ $true_kappa: Group HAS property (T)!") - else - println("κ($name, S) ≥ $true_kappa: Group may NOT HAVE property (T)!") - end - else - println("κ($name, S) ≥ $κ < 0: Tells us nothing about property (T)") - end -end diff --git a/src/GroupAlgebras.jl b/src/GroupAlgebras.jl new file mode 100644 index 0000000..12ebdd2 --- /dev/null +++ b/src/GroupAlgebras.jl @@ -0,0 +1,133 @@ +module GroupAlgebras + +import Base: convert, show, isequal, == +import Base: +, -, *, // +import Base: size, length, norm, rationalize + +export GroupAlgebraElement + + +immutable GroupAlgebraElement{T<:Number} + coefficients::AbstractVector{T} + product_matrix::Array{Int,2} + # basis::Array{Any,1} + + function GroupAlgebraElement(coefficients::AbstractVector, + product_matrix::Array{Int,2}) + + size(product_matrix, 1) == size(product_matrix, 2) || + throw(ArgumentError("Product matrix has to be square")) + new(coefficients, product_matrix) + end +end + +# GroupAlgebraElement(c,pm,b) = GroupAlgebraElement(c,pm) +GroupAlgebraElement{T}(c::AbstractVector{T},pm) = GroupAlgebraElement{T}(c,pm) + +convert{T<:Number}(::Type{T}, X::GroupAlgebraElement) = + GroupAlgebraElement(convert(AbstractVector{T}, X.coefficients), X.product_matrix) + +show{T}(io::IO, X::GroupAlgebraElement{T}) = print(io, + "Element of Group Algebra over $T of length $(length(X)):\n $(X.coefficients)") + + +function isequal{T, S}(X::GroupAlgebraElement{T}, Y::GroupAlgebraElement{S}) + if T != S + warn("Comparing elements with different coefficients Rings!") + end + X.product_matrix == Y.product_matrix || return false + X.coefficients == Y.coefficients || return false + return true +end + +(==)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = isequal(X,Y) + +function add{T<:Number}(X::GroupAlgebraElement{T}, Y::GroupAlgebraElement{T}) + X.product_matrix == Y.product_matrix || throw(ArgumentError( + "Elements don't seem to belong to the same Group Algebra!")) + return GroupAlgebraElement(X.coefficients+Y.coefficients, X.product_matrix) +end + +function add{T<:Number, S<:Number}(X::GroupAlgebraElement{T}, + Y::GroupAlgebraElement{S}) + warn("Adding elements with different base rings!") + return GroupAlgebraElement(+(promote(X.coefficients, Y.coefficients)...), + X.product_matrix) +end + +(+)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = add(X,Y) +(-)(X::GroupAlgebraElement) = GroupAlgebraElement(-X.coefficients, X.product_matrix) +(-)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = add(X,-Y) + +function algebra_multiplication{T<:Number}(X::AbstractVector{T}, Y::AbstractVector{T}, pm::Array{Int,2}) + result = zeros(X) + for (j,y) in enumerate(Y) + if y != zero(T) + for (i, index) in enumerate(pm[:,j]) + if X[i] != zero(T) + index == 0 && throw(ArgumentError("The product don't seem to belong to the span of basis!")) + result[index] += X[i]*y + end + end + end + end + return result +end + +function group_star_multiplication{T<:Number}(X::GroupAlgebraElement{T}, + Y::GroupAlgebraElement{T}) + X.product_matrix == Y.product_matrix || ArgumentError( + "Elements don't seem to belong to the same Group Algebra!") + result = algebra_multiplication(X.coefficients, Y.coefficients, X.product_matrix) + return GroupAlgebraElement(result, X.product_matrix) +end + +function group_star_multiplication{T<:Number, S<:Number}( + X::GroupAlgebraElement{T}, + Y::GroupAlgebraElement{S}) + S == T || warn("Multiplying elements with different base rings!") + return group_star_multiplication(promote(X,Y)...) +end + +(*){T<:Number, S<:Number}(X::GroupAlgebraElement{T}, + Y::GroupAlgebraElement{S}) = group_star_multiplication(X,Y); + +(*){T<:Number}(a::T, X::GroupAlgebraElement{T}) = GroupAlgebraElement( + a*X.coefficients, X.product_matrix) + +function scalar_multiplication{T<:Number, S<:Number}(a::T, + X::GroupAlgebraElement{S}) + promote_type(T,S) == S || warn("Scalar and coefficients are in different rings! Promoting result to $(promote_type(T,S))") + return GroupAlgebraElement(a*X.coefficients, X.product_matrix) +end + +(*){T<:Number}(a::T,X::GroupAlgebraElement) = scalar_multiplication(a, X) + +//{T<:Rational, S<:Rational}(X::GroupAlgebraElement{T}, a::S) = + GroupAlgebraElement(X.coefficients//a, X.product_matrix) + +//{T<:Rational, S<:Integer}(X::GroupAlgebraElement{T}, a::S) = + X//convert(T,a) + +length(X::GroupAlgebraElement) = length(X.coefficients) +size(X::GroupAlgebraElement) = size(X.coefficients) + +function norm(X::GroupAlgebraElement, p=2) + if p == 1 + return sum(abs(X.coefficients)) + elseif p == Inf + return max(abs(X.coefficients)) + else + return norm(X.coefficients, p) + end +end + +ɛ(X::GroupAlgebraElement) = sum(X.coefficients) + +function rationalize{T<:Integer, S<:Number}( + ::Type{T}, X::GroupAlgebraElement{S}; tol=eps(S)) + v = rationalize(T, X.coefficients, tol=tol) + return GroupAlgebraElement(v, X.product_matrix) +end + +end diff --git a/src/Property(T).jl b/src/Property(T).jl new file mode 100644 index 0000000..cebb14f --- /dev/null +++ b/src/Property(T).jl @@ -0,0 +1,135 @@ +module Property(T) + +using GroupAlgebras +import SCS.SCSSolver + +include("sdps.jl") +include("checksolution.jl") + +function pmΔfilenames(name::String) + if !isdir(name) + mkdir(name) + end + prefix = name + pm_filename = joinpath(prefix, "product_matrix.jld") + Δ_coeff_filename = joinpath(prefix, "delta.coefficients.jld") + return pm_filename, Δ_coeff_filename +end + +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 +end + +function ΔandSDPconstraints(name::String) + pm_fname, Δ_fname = pmΔfilenames(name) + f₁ = isfile(pm_fname) + f₂ = isfile(Δ_fname) + if f₁ && f₂ + println("Loading precomputed pm, Δ, sdp_constraints...") + product_matrix = load(pm_fname, "pm") + L = load(Δ_fname, "Δ")[:, 1] + Δ = GroupAlgebraElement(L, Array{Int,2}(product_matrix)) + sdp_constraints = constraints_from_pm(product_matrix) + else + throw(ArgumentError("You need to precompute pm and Δ to load it!")) + end + return Δ, sdp_constraints +end + +function ΔandSDPconstraints(name::String, ID, generating_func::Function) + pm_fname, Δ_fname = pmΔfilenames(name) + Δ, sdp_constraints = ΔandSDPconstraints(ID, generating_func()) + save(pm_fname, "pm", Δ.product_matrix) + save(Δ_fname, "Δ", Δ.coefficients) + return Δ, sdp_constraints +end + +function κandA(name::String) + κ_fname, SDP_fname = κSDPfilenames(name) + f₁ = isfile(κ_fname) + f₂ = isfile(SDP_fname) + if f₁ && f₂ + println("Loading precomputed κ, A...") + κ = load(κ_fname, "κ") + A = load(SDP_fname, "A") + else + throw(ArgumentError("You need to precompute κ and SDP matrix A to load it!")) + end + return κ, A +end + +function κandA(name::String, sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver; upper_bound=Inf) + println("Creating SDP problem...") + @time SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound) + println("Solving SDP problem maximizing κ...") + κ, A = solve_SDP(SDP_problem, solver) + κ_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 + +function check_property_T(name::String, ID, generate_B₄::Function; + verbose=true, tol=1e-6, upper_bound=Inf) + + # solver = MosekSolver(INTPNT_CO_TOL_REL_GAP=tol, QUIET=!verbose) + solver = SCSSolver(eps=tol, max_iters=100000, verbose=verbose) + + @show name + @show verbose + @show tol + + + Δ, sdp_constraints = try + ΔandSDPconstraints(name) + catch err + if isa(err, ArgumentError) + ΔandSDPconstraints(name, ID, generate_B₄) + else + throw(err) + end + end + println("|S| = $(countnz(Δ.coefficients) -1)") + @show length(Δ) + @show 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) + end + end + + @show κ + @show sum(A) + @show maximum(A) + @show minimum(A) + + if κ > 0 + + true_kappa = ℚ_distance_to_positive_cone(Δ, κ, A, tol=tol, verbose=verbose, rational=true) + true_kappa = Float64(trunc(true_kappa,12)) + if true_kappa > 0 + println("κ($name, S) ≥ $true_kappa: Group HAS property (T)!") + else + println("κ($name, S) ≥ $true_kappa: Group may NOT HAVE property (T)!") + end + else + println("κ($name, S) ≥ $κ < 0: Tells us nothing about property (T)") + end +end + +end # module Property(T) diff --git a/src/checksolution.jl b/src/checksolution.jl new file mode 100644 index 0000000..a97c9a1 --- /dev/null +++ b/src/checksolution.jl @@ -0,0 +1,142 @@ +using ProgressMeter +using ValidatedNumerics +import Base: rationalize + +function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, κ::T) + return Δ*Δ - κ*Δ +end + +function square_as_elt(vector, elt) + zzz = zeros(elt.coefficients) + zzz[1:length(vector)] = vector +# new_base_elt = GroupAlgebraElement(zzz, elt.product_matrix) +# return (new_base_elt*new_base_elt).coefficients + return GroupAlgebras.algebra_multiplication(zzz, zzz, elt.product_matrix) +end + +function compute_SOS{T<:Number}(sqrt_matrix::Array{T,2}, + elt::GroupAlgebraElement{T}) + n = size(sqrt_matrix,2) + result = zeros(T, length(elt.coefficients)) + p = Progress(n, 1, "Checking SOS decomposition...", 50) + for i in 1:n + result .+= square_as_elt(sqrt_matrix[:,i], elt) + next!(p) + end + return GroupAlgebraElement{T}(result, elt.product_matrix) +end + +function correct_to_augmentation_ideal{T<:Rational}(sqrt_matrix::Array{T,2}) + sqrt_corrected = similar(sqrt_matrix) + l = size(sqrt_matrix,2) + for i in 1:l + col = view(sqrt_matrix,:,i) + sqrt_corrected[:,i] = col - sum(col)//l + # @assert sum(sqrt_corrected[:,i]) == 0 + end + return sqrt_corrected +end + +function check_solution{T<:Number}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T}; verbose=true, augmented=false) + result = compute_SOS(sqrt_matrix, Δ) + if augmented + epsilon = GroupAlgebras.ɛ(result) + if isa(epsilon, Interval) + @assert 0 in epsilon + elseif isa(epsilon, Rational) + @assert epsilon == 0//1 + else + warn("Does checking for augmentation has meaning for $(typeof(epsilon))?") + end + end + SOS_diff = EOI(Δ, κ) - result + + eoi_SOS_L₁_dist = norm(SOS_diff,1) + + if verbose + @show κ + ɛ_dist = GroupAlgebras.ɛ(SOS_diff) + @printf("ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ≈ %.10f\n", ɛ_dist) + @printf("‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ≈ %.10f\n", eoi_SOS_L₁_dist) + end + + distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist + return distance_to_cone +end + +import ValidatedNumerics.± +function (±)(X::AbstractArray, tol::Real) + r{T}(x::T) = ( x==zero(T) ? @interval(x) : x ± tol) + return r.(X) +end + +(±)(X::GroupAlgebraElement, tol::Real) = GroupAlgebraElement(X.coefficients ± tol, X.product_matrix) + +function Base.rationalize{T<:Integer, S<:Real}(::Type{T}, + X::AbstractArray{S}; tol::Real=eps(eltype(X))) + r(x) = rationalize(T, x, tol=tol) + return r.(X) +end + +ℚ(x, tol::Real) = rationalize(BigInt, x, tol=tol) + +function ℚ_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A; + tol=1e-7, verbose=true, rational=false) + + isapprox(eigvals(A), abs(eigvals(A)), atol=tol) || + warn("The solution matrix doesn't seem to be positive definite!") + @assert A == Symmetric(A) + A_sqrt = real(sqrtm(A)) + + # println("") + # println("Checking in floating-point arithmetic...") + # @time fp_distance = check_solution(κ, A_sqrt, Δ, verbose=verbose) + # println("Floating point distance (to positive cone) ≈ $(Float64(trunc(fp_distance,8)))") + # println("-------------------------------------------------------------") + # println("") + # + # if fp_distance ≤ 0 + # return fp_distance + # end + + println("Checking in interval arithmetic...") + A_sqrtᴵ = A_sqrt ± tol + κᴵ = κ ± tol + Δᴵ = Δ ± tol + @time Interval_distance = check_solution(κᴵ, A_sqrtᴵ, Δᴵ, verbose=verbose) + # @assert isa(ℚ_distance, Rational) + println("The actual distance (to positive cone) is contained in $Interval_distance") + println("-------------------------------------------------------------") + println("") + + if Interval_distance.lo ≤ 0 + return Interval_distance.lo + end + + println("Projecting columns of A_sqrt to the augmentation ideal...") + A_sqrt_ℚ = ℚ(A_sqrt, tol) + A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ) + κ_ℚ = ℚ(κ, tol) + Δ_ℚ = ℚ(Δ, tol) + + A_sqrt_ℚ_augᴵ = A_sqrt_ℚ_aug ± tol + κᴵ = κ_ℚ ± tol + Δᴵ = Δ_ℚ ± tol + @time Interval_dist_to_Σ² = check_solution(κᴵ, A_sqrt_ℚ_augᴵ, Δᴵ, verbose=verbose, augmented=true) + println("The Augmentation-projected actual distance (to positive cone) is contained in $Interval_dist_to_Σ²") + println("-------------------------------------------------------------") + println("") + + if Interval_dist_to_Σ².lo ≤ 0 || !rational + return Interval_dist_to_Σ².lo + else + + println("Checking Projected SOS decomposition in exact rational arithmetic...") + @time ℚ_dist_to_Σ² = check_solution(κ_ℚ, A_sqrt_ℚ_aug, Δ_ℚ, verbose=verbose, augmented=true) + @assert isa(ℚ_dist_to_Σ², Rational) + println("Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))") + println("-------------------------------------------------------------") + return ℚ_dist_to_Σ² + end +end + diff --git a/src/sdps.jl b/src/sdps.jl new file mode 100644 index 0000000..873d41d --- /dev/null +++ b/src/sdps.jl @@ -0,0 +1,87 @@ +using JuMP +import MathProgBase: AbstractMathProgSolver + +function create_product_matrix(basis, limit) + product_matrix = zeros(Int, (limit,limit)) + basis_dict = Dict{Array, Int}(x => i + for (i,x) in enumerate(basis)) + for i in 1:limit + x_inv::eltype(basis) = inv(basis[i]) + for j in 1:limit + w = x_inv*basis[j] + product_matrix[i,j] = basis_dict[w] + # index = findfirst(basis, w) + # index ≠ 0 || throw(ArgumentError("Product is not supported on basis: $w")) + # product_matrix[i,j] = index + end + end + return product_matrix +end + +function constraints_from_pm(pm, total_length) + n = size(pm,1) + constraints = constraints = [Array{Int,1}[] for x in 1:total_length] + for j in 1:n + Threads.@threads for i in 1:n + idx = pm[i,j] + push!(constraints[idx], [i,j]) + end + end + return constraints +end + +constraints_from_pm(pm) = constraints_from_pm(pm, maximum(pm)) + +function splaplacian_coeff(S, basis, n=length(basis)) + result = spzeros(n) + result[1] = float(length(S)) + for s in S + ind = findfirst(basis, s) + result[ind] += -1.0 + end + return result +end + +function laplacian_coeff(S, basis) + return full(splaplacian_coeff(S,basis)) +end + + +function create_SDP_problem(matrix_constraints, Δ::GroupAlgebraElement; upper_bound=Inf) + N = size(Δ.product_matrix,1) + const Δ² = Δ*Δ + @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, κ >= 0.0) + if upper_bound < Inf + JuMP.@constraint(m, κ <= upper_bound) + end + JuMP.@objective(m, Max, κ) + + for (pairs, δ², δ) in zip(matrix_constraints, Δ².coefficients, Δ.coefficients) + JuMP.@constraint(m, sum(A[i,j] for (i,j) in pairs) == δ² - κ*δ) + end + return m +end + +function solve_SDP(SDP_problem, solver) + @show SDP_problem + @show solver + + JuMP.setsolver(SDP_problem, solver); + # @time MathProgBase.writeproblem(SDP_problem, "/tmp/SDP_problem") + solution_status = JuMP.solve(SDP_problem); + + if solution_status != :Optimal + warn("The solver did not solve the problem successfully!") + end + @show solution_status + + κ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :κ)) + A = JuMP.getvalue(JuMP.getvariable(SDP_problem, :A)) + return κ, A +end + diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..3b3d048 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,5 @@ +using Property(T) +using Base.Test + +# write your own tests here +@test 1 == 2