diff --git a/Project.toml b/Project.toml index 90f874a..4c85d12 100644 --- a/Project.toml +++ b/Project.toml @@ -5,11 +5,13 @@ version = "0.1.0" [deps] AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" GroupRings = "0befed6a-bd73-11e8-1e41-a1190947c2f5" Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73" Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a" diff --git a/src/1712.07167.jl b/src/1712.07167.jl index 5b8ec0b..e02a2fd 100644 --- a/src/1712.07167.jl +++ b/src/1712.07167.jl @@ -1,3 +1,5 @@ +using Printf + ############################################################################### # # Settings and filenames @@ -84,14 +86,15 @@ end function computeλandP(sett::Naive, Δ::GroupRingElem; solverlog=tempname()*".log") - info("Creating SDP problem...") - info(Base.repr(SDP_problem)) + @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) - @show status + @info("Solver's status: $status") + save(filename(sett, :warmstart), "warmstart", ws, "P", P, "λ", λ) return λ, P @@ -108,17 +111,17 @@ function computeλandP(sett::Symmetrized, Δ::GroupRingElem; orbit_data = load(filename(sett, :OrbitData), "OrbitData") orbit_data = decimate(orbit_data) - info("Creating SDP problem...") + @info("Creating SDP problem...") - info(Base.repr(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) - @show status + @info("Solver's status: $status") + save(filename(sett, :warmstart), "warmstart", ws, "Ps", Ps, "λ", λ) - info("Reconstructing P...") @time P = reconstruct(Ps, orbit_data) @@ -146,7 +149,6 @@ function distance_to_positive_cone(Δ::GroupRingElem, λ, Q; R::Int=2) @info("Floating point distance (to positive cone) ≈") @info("$(@sprintf("%.10f", distance))") - @info("") if distance ≤ 0 return distance @@ -154,15 +156,14 @@ function distance_to_positive_cone(Δ::GroupRingElem, λ, Q; R::Int=2) @info("------------------------------------------------------------") @info("Checking in interval arithmetic...") - @info("λ ∈ $λ") - λ = @interval(λ) + @info("λ ∈ $λ") 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!")) + @info((check ? "They do." : "FAILED!")) check || @warn("The following numbers are meaningless!") @time residual = eoi - compute_SOS(parent(eoi), Q) @@ -174,7 +175,7 @@ function distance_to_positive_cone(Δ::GroupRingElem, λ, Q; R::Int=2) @info("Interval distance (to positive cone) ∈") @info("$(distance)") - @info("") + @info("------------------------------------------------------------") return distance.lo end @@ -192,11 +193,11 @@ 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)!") + @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)") + @info("λ($(sett.name), S) ≥ $sgap < 0: Tells us nothing about property (T)") return false end @@ -222,20 +223,20 @@ function check_property_T(sett::Settings) save(filename(sett, :solution), "λ", λ, "P", P) if λ < 0 - warn("Solver did not produce a valid solution!") + @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))") + @info("λ = $λ") + @info("sum(P) = $(sum(P))") + @info("maximum(P) = $(maximum(P))") + @info("minimum(P) = $(minimum(P))") isapprox(eigvals(P), abs.(eigvals(P))) || @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) + @time Q = real(sqrt( (P.+ P')./2 )) + sgap = distance_to_positive_cone(Δ, λ, Q, R=2*sett.radius) return interpret_results(sett, sgap) end diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 156dd4c..5413a44 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -2,6 +2,11 @@ __precompile__() module PropertyT using AbstractAlgebra +using LinearAlgebra +using SparseArrays +using Markdown +using Dates + using Groups using GroupRings @@ -17,5 +22,6 @@ include("RGprojections.jl") include("orbitdata.jl") include("sos_sdps.jl") include("checksolution.jl") +include("1712.07167.jl") end # module Property(T) diff --git a/src/RGprojections.jl b/src/RGprojections.jl index 3a18463..7da54a1 100644 --- a/src/RGprojections.jl +++ b/src/RGprojections.jl @@ -3,6 +3,7 @@ module Projections using AbstractAlgebra using Groups using GroupRings +using Markdown export PermCharacter, DirectProdCharacter, rankOne_projections @@ -200,13 +201,13 @@ end # ############################################################################## -doc""" +@doc doc""" products(X::Vector{GroupElem}, Y::Vector{GroupElem}, op=*) > Returns a vector of all possible products (or `op(x,y)`), where $x\in X$ and > $y\in Y$ are group elements. You may specify which operation is used when > forming 'products' by adding `op` (which is `*` by default). """ -function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) +function products(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) where {T<:GroupElem} result = Vector{T}() seen = Set{T}() for x in X @@ -221,7 +222,7 @@ function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=* return result end -doc""" +@doc doc""" generateGroup(gens::Vector{GroupElem}, r=2, Id=parent(first(gens))(), op=*) > Produces all elements of a group generated by elements in `gens` in ball of > radius `r` (word-length metric induced by `gens`). @@ -230,7 +231,7 @@ doc""" > The identity element `Id` and binary operation function `op` can be supplied > to e.g. take advantage of additive group structure. """ -function generateGroup{T<:GroupElem}(gens::Vector{T}, r=2, Id::T=parent(first(gens))(), op=*) +function generateGroup(gens::Vector{T}, r=2, Id::T=parent(first(gens))(), op=*) where {T<:GroupElem} n = 0 R = 1 elts = gens diff --git a/src/laplacians.jl b/src/laplacians.jl index 207d1ea..eacb8aa 100644 --- a/src/laplacians.jl +++ b/src/laplacians.jl @@ -33,12 +33,12 @@ function Laplacian(S::Vector{E}, radius) where E<:AbstractAlgebra.GroupElem end function Laplacian(S, Id, radius) - info("Generating metric ball of radius $(2radius)...") + @info("Generating metric ball of radius $(2radius)...") @time E_R, sizes = Groups.generate_balls(S, Id, radius=2radius) - info("Generated balls of sizes $sizes.") + @info("Generated balls of sizes $sizes.") - info("Creating product matrix...") @time pm = GroupRings.create_pm(E_R, GroupRings.reverse_dict(E_R), sizes[radius]; twisted=true) + @info("Creating product matrix...") RG = GroupRing(parent(Id), E_R, pm) Δ = spLaplacian(RG, S) @@ -52,7 +52,7 @@ end function loadGRElem(fname::String, G::Group) if isfile(fname) - info("Loading precomputed Δ...") + @info("Loading precomputed Δ...") coeffs, pm = load(fname, "coeffs", "pm") RG = GroupRing(G, pm) Δ = GroupRingElem(coeffs, RG) diff --git a/src/orbitdata.jl b/src/orbitdata.jl index d18922e..8b2da49 100644 --- a/src/orbitdata.jl +++ b/src/orbitdata.jl @@ -12,25 +12,25 @@ struct OrbitData{T<:AbstractArray{Float64, 2}, GEl<:GroupElem, P<:perm} end function OrbitData(RG::GroupRing, autS::Group, verbose=true) - verbose && info("Decomposing basis of RG into orbits of $(autS)") + verbose && @info("Decomposing basis of RG into orbits of $(autS)") @time orbs = orbit_decomposition(autS, RG.basis, RG.basis_dict) @assert sum(length(o) for o in orbs) == length(RG.basis) - verbose && info("The action has $(length(orbs)) orbits") + verbose && @info("The action has $(length(orbs)) orbits") - verbose && info("Projections in the Group Ring of AutS") @time autS_mps = Projections.rankOne_projections(GroupRing(autS)) + verbose && @info("Projections in the Group Ring of AutS = $autS") - verbose && info("AutS-action matrix representatives") + verbose && @info("AutS-action matrix representatives") @time preps = perm_reps(autS, RG.basis[1:size(RG.pm,1)], RG.basis_dict) @time mreps = matrix_reps(preps) - verbose && info("Projection matrices Uπs") + verbose && @info("Projection matrices Uπs") @time Uπs = [orthSVD(matrix_repr(p, mreps)) for p in autS_mps] multiplicities = size.(Uπs,2) - verbose && info("multiplicities = $multiplicities") + verbose && @info("multiplicities = $multiplicities") dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps] - verbose && info("dimensions = $dimensions") + verbose && @info("dimensions = $dimensions") @assert dot(multiplicities, dimensions) == size(RG.pm,1) return OrbitData(orbs, preps, Uπs, dimensions) @@ -43,32 +43,32 @@ function decimate(od::OrbitData) #dimensions of the corresponding πs: dims = od.dims[nzros] - return OrbitData(od.orbits, od.preps, full.(Us), dims); + return OrbitData(od.orbits, od.preps, Array{Float64}.(Us), dims); end function orthSVD(M::AbstractMatrix{T}) where {T<:AbstractFloat} - M = full(M) - fact = svdfact(M) - M_rank = sum(fact[:S] .> maximum(size(M))*eps(T)) - return fact[:U][:,1:M_rank] + M = Matrix(M) + fact = svd(M) + M_rank = sum(fact.S .> maximum(size(M))*eps(T)) + return fact.U[:,1:M_rank] end function orbit_decomposition(G::Group, E::Vector, rdict=GroupRings.reverse_dict(E)) elts = collect(elements(G)) - tovisit = trues(E); + tovisit = trues(size(E)); orbits = Vector{Vector{Int}}() orbit = zeros(Int, length(elts)) - for i in 1:endof(E) + for i in eachindex(E) if tovisit[i] g = E[i] - Threads.@threads for j in 1:length(elts) + Threads.@threads for j in eachindex(elts) orbit[j] = rdict[elts[j](g)] end - tovisit[orbit] = false + tovisit[orbit] .= false push!(orbits, unique(orbit)) end end @@ -82,9 +82,9 @@ end ############################################################################### dens(M::SparseMatrixCSC) = nnz(M)/length(M) -dens(M::AbstractArray) = countnz(M)/length(M) +dens(M::AbstractArray) = count(!iszero, M)/length(M) -function sparsify!{Tv,Ti}(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv); verbose=false) +function sparsify!(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv); verbose=false) where {Tv,Ti} densM = dens(M) for i in eachindex(M.nzval) @@ -95,16 +95,16 @@ function sparsify!{Tv,Ti}(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv); verbose=false) dropzeros!(M) if verbose - info("Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M), 20), " ($(nnz(M)) non-zeros)") + @info("Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M), 20), " ($(nnz(M)) non-zeros)") end return M end -function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); verbose=false) +function sparsify!(M::AbstractArray{T}, eps=eps(T); verbose=false) where T densM = dens(M) if verbose - info("Sparsifying $(size(M))-matrix... ") + @info("Sparsifying $(size(M))-matrix... ") end for n in eachindex(M) @@ -114,13 +114,15 @@ function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); verbose=false) end if verbose - info("$(rpad(densM, 20)) → $(rpad(dens(M),20))), ($(countnz(M)) non-zeros)") + @info("$(rpad(densM, 20)) → $(rpad(dens(M),20))), ($(count(!iszero, M)) non-zeros)") end return sparse(M) end -sparsify{T}(U::AbstractArray{T}, tol=eps(T); verbose=false) = sparsify!(deepcopy(U), tol, verbose=verbose) +function sparsify(U::AbstractArray{T}, tol=eps(T); verbose=false) where T + return sparsify!(deepcopy(U), tol, verbose=verbose) +end ############################################################################### # @@ -129,7 +131,7 @@ sparsify{T}(U::AbstractArray{T}, tol=eps(T); verbose=false) = sparsify!(deepcopy ############################################################################### function perm_repr(g::GroupElem, E::Vector, E_dict) - p = Vector{Int}(length(E)) + p = Vector{Int}(undef, length(E)) for (i,elt) in enumerate(E) p[i] = E_dict[g(elt)] end @@ -139,7 +141,7 @@ end function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E)) elts = collect(elements(G)) l = length(elts) - preps = Vector{perm}(l) + preps = Vector{perm}(undef, l) permG = PermutationGroup(length(E)) @@ -151,13 +153,13 @@ function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E)) end function matrix_repr(x::GroupRingElem, mreps::Dict) - nzeros = findn(x.coeffs) + nzeros = findall(!iszero, x.coeffs) return sum(x[i].*mreps[parent(x).basis[i]] for i in nzeros) end function matrix_reps(preps::Dict{T,perm{I}}) where {T<:GroupElem, I<:Integer} kk = collect(keys(preps)) - mreps = Vector{SparseMatrixCSC{Float64, Int}}(length(kk)) + mreps = Vector{SparseMatrixCSC{Float64, Int}}(undef, length(kk)) Threads.@threads for i in 1:length(kk) mreps[i] = AbstractAlgebra.matrix_repr(preps[kk[i]]) end diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index a304ef2..fecd929 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -75,7 +75,7 @@ function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData Ns = size.(data.Uπs, 2) m = JuMP.Model(); - P = Vector{Matrix{JuMP.Variable}}(length(Ns)) + P = Vector{Matrix{JuMP.Variable}}(undef, length(Ns)) for (k,s) in enumerate(Ns) P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) @@ -87,7 +87,7 @@ function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData JuMP.@constraint(m, λ <= upper_bound) end - info("Adding $(length(data.orbits)) constraints... ") + @info("Adding $(length(data.orbits)) constraints... ") @time addconstraints!(m,P,λ,X,orderunit, data) @@ -96,7 +96,7 @@ function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData end function constraintLHS!(M, cnstr, Us, Ust, dims, eps=1000*eps(eltype(first(M)))) - for π in 1:endof(Us) + for π in eachindex(Us) M[π] = PropertyT.sparsify!(dims[π].*Ust[π]*cnstr*Us[π], eps) end end @@ -112,13 +112,13 @@ function addconstraints!(m::JuMP.Model, cnstrs = constraints(parent(X).pm) orb_cnstr = spzeros(Float64, size(parent(X).pm)...) - M = [Array{Float64}(n,n) for n in size.(UπsT,1)] + M = [Array{Float64}(undef, n,n) for n in size.(UπsT,1)] for (t, orbit) in enumerate(data.orbits) orbit_constraint!(orb_cnstr, cnstrs, orbit) constraintLHS!(M, orb_cnstr, data.Uπs, UπsT, data.dims) - lhs = @expression(m, sum(vecdot(M[π], P[π]) for π in 1:endof(data.Uπs))) + lhs = @expression(m, sum(dot(M[π], P[π]) for π in eachindex(data.Uπs))) x, u = X_orb[t], orderunit_orb[t] JuMP.@constraint(m, lhs == x - λ*u) end @@ -198,8 +198,8 @@ function solve(solverlog::String, model::JuMP.Model, varλ::JuMP.Variable, varP, isdir(dirname(solverlog)) || mkpath(dirname(solverlog)) + Base.flush(Base.stdout) status, (λ, P, warmstart) = open(solverlog, "a+") do logfile - Base.Libc.flush_cstdio() redirect_stdout(logfile) do status, (λ, P, warmstart) = PropertyT.solve(model, varλ, varP, warmstart) Base.Libc.flush_cstdio() @@ -224,7 +224,7 @@ function fillfrominternal!(m::JuMP.Model, traits) m.objBound = NaN m.objVal = NaN m.colVal = fill(NaN, numCols) - m.linconstrDuals = Array{Float64}(0) + m.linconstrDuals = Array{Float64}(undef, 0) discrete = (traits.int || traits.sos) @@ -262,7 +262,7 @@ function fillfrominternal!(m::JuMP.Model, traits) @assert length(infray) == numRows infray catch - suppress_warnings || warn("Infeasibility ray (Farkas proof) not available") + @warn("Infeasibility ray (Farkas proof) not available") fill(NaN, numRows) end elseif stat == :Unbounded @@ -271,7 +271,7 @@ function fillfrominternal!(m::JuMP.Model, traits) @assert length(unbdray) == numCols unbdray catch - suppress_warnings || warn("Unbounded ray not available") + @warn("Unbounded ray not available") fill(NaN, numCols) end end @@ -292,7 +292,10 @@ function fillfrominternal!(m::JuMP.Model, traits) # Do a separate try since getobjval could work while getobjbound does not and vice versa objBound = MathProgBase.getobjbound(m.internalModel) + m.obj.aff.constant m.objBound = objBound + catch + @warn("objBound could not be obtained") end + try objVal = MathProgBase.getobjval(m.internalModel) + m.obj.aff.constant colVal = MathProgBase.getsolution(m.internalModel)[1:numCols] @@ -304,8 +307,15 @@ function fillfrominternal!(m::JuMP.Model, traits) # Don't corrupt the answers if one of the above two calls fails m.objVal = objVal m.colVal = colVal + catch + @warn("objVal/colVal could not be obtained") end end + + if traits.conic && m.objSense == :Max + m.objBound = -1 * (m.objBound - m.obj.aff.constant) + m.obj.aff.constant + m.objVal = -1 * (m.objVal - m.obj.aff.constant) + m.obj.aff.constant + end return stat end