mirror of
https://github.com/kalmarek/PropertyT.jl.git
synced 2024-11-26 00:55:27 +01:00
fix: deprecations for julia-0.7/1.0
This commit is contained in:
parent
0ee12f76a5
commit
10d319f48b
@ -5,11 +5,13 @@ version = "0.1.0"
|
|||||||
|
|
||||||
[deps]
|
[deps]
|
||||||
AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
|
AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
|
||||||
|
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
|
||||||
GroupRings = "0befed6a-bd73-11e8-1e41-a1190947c2f5"
|
GroupRings = "0befed6a-bd73-11e8-1e41-a1190947c2f5"
|
||||||
Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
|
Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
|
||||||
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
|
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
|
||||||
JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
|
JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
|
||||||
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
|
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
|
||||||
|
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
|
||||||
Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
|
Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
|
||||||
MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
|
MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
|
||||||
Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
|
Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
|
||||||
|
@ -1,3 +1,5 @@
|
|||||||
|
using Printf
|
||||||
|
|
||||||
###############################################################################
|
###############################################################################
|
||||||
#
|
#
|
||||||
# Settings and filenames
|
# Settings and filenames
|
||||||
@ -84,14 +86,15 @@ end
|
|||||||
function computeλandP(sett::Naive, Δ::GroupRingElem;
|
function computeλandP(sett::Naive, Δ::GroupRingElem;
|
||||||
solverlog=tempname()*".log")
|
solverlog=tempname()*".log")
|
||||||
|
|
||||||
info("Creating SDP problem...")
|
@info("Creating SDP problem...")
|
||||||
info(Base.repr(SDP_problem))
|
|
||||||
SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, upper_bound=sett.upper_bound)
|
SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, upper_bound=sett.upper_bound)
|
||||||
JuMP.setsolver(SDP_problem, sett.solver)
|
JuMP.setsolver(SDP_problem, sett.solver)
|
||||||
|
@info(Base.repr(SDP_problem))
|
||||||
|
|
||||||
ws = warmstart(sett)
|
ws = warmstart(sett)
|
||||||
@time status, (λ, P, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws)
|
@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, "λ", λ)
|
save(filename(sett, :warmstart), "warmstart", ws, "P", P, "λ", λ)
|
||||||
|
|
||||||
return λ, P
|
return λ, P
|
||||||
@ -108,17 +111,17 @@ function computeλandP(sett::Symmetrized, Δ::GroupRingElem;
|
|||||||
orbit_data = load(filename(sett, :OrbitData), "OrbitData")
|
orbit_data = load(filename(sett, :OrbitData), "OrbitData")
|
||||||
orbit_data = decimate(orbit_data)
|
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)
|
SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, orbit_data, upper_bound=sett.upper_bound)
|
||||||
JuMP.setsolver(SDP_problem, sett.solver)
|
JuMP.setsolver(SDP_problem, sett.solver)
|
||||||
|
@info(Base.repr(SDP_problem))
|
||||||
|
|
||||||
ws = warmstart(sett)
|
ws = warmstart(sett)
|
||||||
@time status, (λ, Ps, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws)
|
@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, "λ", λ)
|
save(filename(sett, :warmstart), "warmstart", ws, "Ps", Ps, "λ", λ)
|
||||||
|
|
||||||
info("Reconstructing P...")
|
info("Reconstructing P...")
|
||||||
@time P = reconstruct(Ps, orbit_data)
|
@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("Floating point distance (to positive cone) ≈")
|
||||||
@info("$(@sprintf("%.10f", distance))")
|
@info("$(@sprintf("%.10f", distance))")
|
||||||
@info("")
|
|
||||||
|
|
||||||
if distance ≤ 0
|
if distance ≤ 0
|
||||||
return distance
|
return distance
|
||||||
@ -154,15 +156,14 @@ function distance_to_positive_cone(Δ::GroupRingElem, λ, Q; R::Int=2)
|
|||||||
|
|
||||||
@info("------------------------------------------------------------")
|
@info("------------------------------------------------------------")
|
||||||
@info("Checking in interval arithmetic...")
|
@info("Checking in interval arithmetic...")
|
||||||
@info("λ ∈ $λ")
|
|
||||||
|
|
||||||
λ = @interval(λ)
|
λ = @interval(λ)
|
||||||
|
@info("λ ∈ $λ")
|
||||||
eoi = Δ^2 - λ*Δ
|
eoi = Δ^2 - λ*Δ
|
||||||
|
|
||||||
@info("Projecting columns of Q to the augmentation ideal...")
|
@info("Projecting columns of Q to the augmentation ideal...")
|
||||||
@time Q, check = augIdproj(Interval, Q)
|
@time Q, check = augIdproj(Interval, Q)
|
||||||
@info("Checking that sum of every column contains 0.0... ")
|
@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!")
|
check || @warn("The following numbers are meaningless!")
|
||||||
|
|
||||||
@time residual = eoi - compute_SOS(parent(eoi), Q)
|
@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("Interval distance (to positive cone) ∈")
|
||||||
@info("$(distance)")
|
@info("$(distance)")
|
||||||
@info("")
|
@info("------------------------------------------------------------")
|
||||||
|
|
||||||
return distance.lo
|
return distance.lo
|
||||||
end
|
end
|
||||||
@ -192,11 +193,11 @@ function interpret_results(sett::Settings, sgap::Number)
|
|||||||
if sgap > 0
|
if sgap > 0
|
||||||
Kazhdan_κ = Kazhdan(sgap, length(sett.S))
|
Kazhdan_κ = Kazhdan(sgap, length(sett.S))
|
||||||
if Kazhdan_κ > 0
|
if Kazhdan_κ > 0
|
||||||
info("κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
@info("κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
||||||
return true
|
return true
|
||||||
end
|
end
|
||||||
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
|
return false
|
||||||
end
|
end
|
||||||
|
|
||||||
@ -222,20 +223,20 @@ function check_property_T(sett::Settings)
|
|||||||
save(filename(sett, :solution), "λ", λ, "P", P)
|
save(filename(sett, :solution), "λ", λ, "P", P)
|
||||||
|
|
||||||
if λ < 0
|
if λ < 0
|
||||||
warn("Solver did not produce a valid solution!")
|
@warn("Solver did not produce a valid solution!")
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
info("λ = $λ")
|
@info("λ = $λ")
|
||||||
info("sum(P) = $(sum(P))")
|
@info("sum(P) = $(sum(P))")
|
||||||
info("maximum(P) = $(maximum(P))")
|
@info("maximum(P) = $(maximum(P))")
|
||||||
info("minimum(P) = $(minimum(P))")
|
@info("minimum(P) = $(minimum(P))")
|
||||||
|
|
||||||
isapprox(eigvals(P), abs.(eigvals(P))) ||
|
isapprox(eigvals(P), abs.(eigvals(P))) ||
|
||||||
@warn("The solution matrix doesn't seem to be positive definite!")
|
@warn("The solution matrix doesn't seem to be positive definite!")
|
||||||
|
|
||||||
@time Q = real(sqrtm((P+P')/2))
|
@time Q = real(sqrt( (P.+ P')./2 ))
|
||||||
sgap = distance_to_positive_cone(Δ, λ, Q, wlen=2*sett.radius)
|
sgap = distance_to_positive_cone(Δ, λ, Q, R=2*sett.radius)
|
||||||
|
|
||||||
return interpret_results(sett, sgap)
|
return interpret_results(sett, sgap)
|
||||||
end
|
end
|
||||||
|
@ -2,6 +2,11 @@ __precompile__()
|
|||||||
module PropertyT
|
module PropertyT
|
||||||
|
|
||||||
using AbstractAlgebra
|
using AbstractAlgebra
|
||||||
|
using LinearAlgebra
|
||||||
|
using SparseArrays
|
||||||
|
using Markdown
|
||||||
|
using Dates
|
||||||
|
|
||||||
using Groups
|
using Groups
|
||||||
using GroupRings
|
using GroupRings
|
||||||
|
|
||||||
@ -17,5 +22,6 @@ include("RGprojections.jl")
|
|||||||
include("orbitdata.jl")
|
include("orbitdata.jl")
|
||||||
include("sos_sdps.jl")
|
include("sos_sdps.jl")
|
||||||
include("checksolution.jl")
|
include("checksolution.jl")
|
||||||
|
include("1712.07167.jl")
|
||||||
|
|
||||||
end # module Property(T)
|
end # module Property(T)
|
||||||
|
@ -3,6 +3,7 @@ module Projections
|
|||||||
using AbstractAlgebra
|
using AbstractAlgebra
|
||||||
using Groups
|
using Groups
|
||||||
using GroupRings
|
using GroupRings
|
||||||
|
using Markdown
|
||||||
|
|
||||||
export PermCharacter, DirectProdCharacter, rankOne_projections
|
export PermCharacter, DirectProdCharacter, rankOne_projections
|
||||||
|
|
||||||
@ -200,13 +201,13 @@ end
|
|||||||
#
|
#
|
||||||
##############################################################################
|
##############################################################################
|
||||||
|
|
||||||
doc"""
|
@doc doc"""
|
||||||
products(X::Vector{GroupElem}, Y::Vector{GroupElem}, op=*)
|
products(X::Vector{GroupElem}, Y::Vector{GroupElem}, op=*)
|
||||||
> Returns a vector of all possible products (or `op(x,y)`), where $x\in X$ and
|
> 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
|
> $y\in Y$ are group elements. You may specify which operation is used when
|
||||||
> forming 'products' by adding `op` (which is `*` by default).
|
> 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}()
|
result = Vector{T}()
|
||||||
seen = Set{T}()
|
seen = Set{T}()
|
||||||
for x in X
|
for x in X
|
||||||
@ -221,7 +222,7 @@ function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=*
|
|||||||
return result
|
return result
|
||||||
end
|
end
|
||||||
|
|
||||||
doc"""
|
@doc doc"""
|
||||||
generateGroup(gens::Vector{GroupElem}, r=2, Id=parent(first(gens))(), op=*)
|
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
|
> Produces all elements of a group generated by elements in `gens` in ball of
|
||||||
> radius `r` (word-length metric induced by `gens`).
|
> 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
|
> The identity element `Id` and binary operation function `op` can be supplied
|
||||||
> to e.g. take advantage of additive group structure.
|
> 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
|
n = 0
|
||||||
R = 1
|
R = 1
|
||||||
elts = gens
|
elts = gens
|
||||||
|
@ -33,12 +33,12 @@ function Laplacian(S::Vector{E}, radius) where E<:AbstractAlgebra.GroupElem
|
|||||||
end
|
end
|
||||||
|
|
||||||
function Laplacian(S, Id, radius)
|
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)
|
@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)
|
@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)
|
RG = GroupRing(parent(Id), E_R, pm)
|
||||||
Δ = spLaplacian(RG, S)
|
Δ = spLaplacian(RG, S)
|
||||||
@ -52,7 +52,7 @@ end
|
|||||||
|
|
||||||
function loadGRElem(fname::String, G::Group)
|
function loadGRElem(fname::String, G::Group)
|
||||||
if isfile(fname)
|
if isfile(fname)
|
||||||
info("Loading precomputed Δ...")
|
@info("Loading precomputed Δ...")
|
||||||
coeffs, pm = load(fname, "coeffs", "pm")
|
coeffs, pm = load(fname, "coeffs", "pm")
|
||||||
RG = GroupRing(G, pm)
|
RG = GroupRing(G, pm)
|
||||||
Δ = GroupRingElem(coeffs, RG)
|
Δ = GroupRingElem(coeffs, RG)
|
||||||
|
@ -12,25 +12,25 @@ struct OrbitData{T<:AbstractArray{Float64, 2}, GEl<:GroupElem, P<:perm}
|
|||||||
end
|
end
|
||||||
|
|
||||||
function OrbitData(RG::GroupRing, autS::Group, verbose=true)
|
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)
|
@time orbs = orbit_decomposition(autS, RG.basis, RG.basis_dict)
|
||||||
@assert sum(length(o) for o in orbs) == length(RG.basis)
|
@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))
|
@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 preps = perm_reps(autS, RG.basis[1:size(RG.pm,1)], RG.basis_dict)
|
||||||
@time mreps = matrix_reps(preps)
|
@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]
|
@time Uπs = [orthSVD(matrix_repr(p, mreps)) for p in autS_mps]
|
||||||
|
|
||||||
multiplicities = size.(Uπs,2)
|
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]
|
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)
|
@assert dot(multiplicities, dimensions) == size(RG.pm,1)
|
||||||
|
|
||||||
return OrbitData(orbs, preps, Uπs, dimensions)
|
return OrbitData(orbs, preps, Uπs, dimensions)
|
||||||
@ -43,32 +43,32 @@ function decimate(od::OrbitData)
|
|||||||
#dimensions of the corresponding πs:
|
#dimensions of the corresponding πs:
|
||||||
dims = od.dims[nzros]
|
dims = od.dims[nzros]
|
||||||
|
|
||||||
return OrbitData(od.orbits, od.preps, full.(Us), dims);
|
return OrbitData(od.orbits, od.preps, Array{Float64}.(Us), dims);
|
||||||
end
|
end
|
||||||
|
|
||||||
function orthSVD(M::AbstractMatrix{T}) where {T<:AbstractFloat}
|
function orthSVD(M::AbstractMatrix{T}) where {T<:AbstractFloat}
|
||||||
M = full(M)
|
M = Matrix(M)
|
||||||
fact = svdfact(M)
|
fact = svd(M)
|
||||||
M_rank = sum(fact[:S] .> maximum(size(M))*eps(T))
|
M_rank = sum(fact.S .> maximum(size(M))*eps(T))
|
||||||
return fact[:U][:,1:M_rank]
|
return fact.U[:,1:M_rank]
|
||||||
end
|
end
|
||||||
|
|
||||||
function orbit_decomposition(G::Group, E::Vector, rdict=GroupRings.reverse_dict(E))
|
function orbit_decomposition(G::Group, E::Vector, rdict=GroupRings.reverse_dict(E))
|
||||||
|
|
||||||
elts = collect(elements(G))
|
elts = collect(elements(G))
|
||||||
|
|
||||||
tovisit = trues(E);
|
tovisit = trues(size(E));
|
||||||
orbits = Vector{Vector{Int}}()
|
orbits = Vector{Vector{Int}}()
|
||||||
|
|
||||||
orbit = zeros(Int, length(elts))
|
orbit = zeros(Int, length(elts))
|
||||||
|
|
||||||
for i in 1:endof(E)
|
for i in eachindex(E)
|
||||||
if tovisit[i]
|
if tovisit[i]
|
||||||
g = E[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)]
|
orbit[j] = rdict[elts[j](g)]
|
||||||
end
|
end
|
||||||
tovisit[orbit] = false
|
tovisit[orbit] .= false
|
||||||
push!(orbits, unique(orbit))
|
push!(orbits, unique(orbit))
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
@ -82,9 +82,9 @@ end
|
|||||||
###############################################################################
|
###############################################################################
|
||||||
|
|
||||||
dens(M::SparseMatrixCSC) = nnz(M)/length(M)
|
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)
|
densM = dens(M)
|
||||||
for i in eachindex(M.nzval)
|
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)
|
dropzeros!(M)
|
||||||
|
|
||||||
if verbose
|
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
|
end
|
||||||
|
|
||||||
return M
|
return M
|
||||||
end
|
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)
|
densM = dens(M)
|
||||||
if verbose
|
if verbose
|
||||||
info("Sparsifying $(size(M))-matrix... ")
|
@info("Sparsifying $(size(M))-matrix... ")
|
||||||
end
|
end
|
||||||
|
|
||||||
for n in eachindex(M)
|
for n in eachindex(M)
|
||||||
@ -114,13 +114,15 @@ function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); verbose=false)
|
|||||||
end
|
end
|
||||||
|
|
||||||
if verbose
|
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
|
end
|
||||||
|
|
||||||
return sparse(M)
|
return sparse(M)
|
||||||
end
|
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)
|
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)
|
for (i,elt) in enumerate(E)
|
||||||
p[i] = E_dict[g(elt)]
|
p[i] = E_dict[g(elt)]
|
||||||
end
|
end
|
||||||
@ -139,7 +141,7 @@ end
|
|||||||
function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
|
function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
|
||||||
elts = collect(elements(G))
|
elts = collect(elements(G))
|
||||||
l = length(elts)
|
l = length(elts)
|
||||||
preps = Vector{perm}(l)
|
preps = Vector{perm}(undef, l)
|
||||||
|
|
||||||
permG = PermutationGroup(length(E))
|
permG = PermutationGroup(length(E))
|
||||||
|
|
||||||
@ -151,13 +153,13 @@ function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
|
|||||||
end
|
end
|
||||||
|
|
||||||
function matrix_repr(x::GroupRingElem, mreps::Dict)
|
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)
|
return sum(x[i].*mreps[parent(x).basis[i]] for i in nzeros)
|
||||||
end
|
end
|
||||||
|
|
||||||
function matrix_reps(preps::Dict{T,perm{I}}) where {T<:GroupElem, I<:Integer}
|
function matrix_reps(preps::Dict{T,perm{I}}) where {T<:GroupElem, I<:Integer}
|
||||||
kk = collect(keys(preps))
|
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)
|
Threads.@threads for i in 1:length(kk)
|
||||||
mreps[i] = AbstractAlgebra.matrix_repr(preps[kk[i]])
|
mreps[i] = AbstractAlgebra.matrix_repr(preps[kk[i]])
|
||||||
end
|
end
|
||||||
|
@ -75,7 +75,7 @@ function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData
|
|||||||
Ns = size.(data.Uπs, 2)
|
Ns = size.(data.Uπs, 2)
|
||||||
m = JuMP.Model();
|
m = JuMP.Model();
|
||||||
|
|
||||||
P = Vector{Matrix{JuMP.Variable}}(length(Ns))
|
P = Vector{Matrix{JuMP.Variable}}(undef, length(Ns))
|
||||||
|
|
||||||
for (k,s) in enumerate(Ns)
|
for (k,s) in enumerate(Ns)
|
||||||
P[k] = JuMP.@variable(m, [i=1:s, j=1:s])
|
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)
|
JuMP.@constraint(m, λ <= upper_bound)
|
||||||
end
|
end
|
||||||
|
|
||||||
info("Adding $(length(data.orbits)) constraints... ")
|
@info("Adding $(length(data.orbits)) constraints... ")
|
||||||
|
|
||||||
@time addconstraints!(m,P,λ,X,orderunit, data)
|
@time addconstraints!(m,P,λ,X,orderunit, data)
|
||||||
|
|
||||||
@ -96,7 +96,7 @@ function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData
|
|||||||
end
|
end
|
||||||
|
|
||||||
function constraintLHS!(M, cnstr, Us, Ust, dims, eps=1000*eps(eltype(first(M))))
|
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)
|
M[π] = PropertyT.sparsify!(dims[π].*Ust[π]*cnstr*Us[π], eps)
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
@ -112,13 +112,13 @@ function addconstraints!(m::JuMP.Model,
|
|||||||
cnstrs = constraints(parent(X).pm)
|
cnstrs = constraints(parent(X).pm)
|
||||||
orb_cnstr = spzeros(Float64, size(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)
|
for (t, orbit) in enumerate(data.orbits)
|
||||||
orbit_constraint!(orb_cnstr, cnstrs, orbit)
|
orbit_constraint!(orb_cnstr, cnstrs, orbit)
|
||||||
constraintLHS!(M, orb_cnstr, data.Uπs, UπsT, data.dims)
|
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]
|
x, u = X_orb[t], orderunit_orb[t]
|
||||||
JuMP.@constraint(m, lhs == x - λ*u)
|
JuMP.@constraint(m, lhs == x - λ*u)
|
||||||
end
|
end
|
||||||
@ -198,8 +198,8 @@ function solve(solverlog::String, model::JuMP.Model, varλ::JuMP.Variable, varP,
|
|||||||
|
|
||||||
isdir(dirname(solverlog)) || mkpath(dirname(solverlog))
|
isdir(dirname(solverlog)) || mkpath(dirname(solverlog))
|
||||||
|
|
||||||
|
Base.flush(Base.stdout)
|
||||||
status, (λ, P, warmstart) = open(solverlog, "a+") do logfile
|
status, (λ, P, warmstart) = open(solverlog, "a+") do logfile
|
||||||
Base.Libc.flush_cstdio()
|
|
||||||
redirect_stdout(logfile) do
|
redirect_stdout(logfile) do
|
||||||
status, (λ, P, warmstart) = PropertyT.solve(model, varλ, varP, warmstart)
|
status, (λ, P, warmstart) = PropertyT.solve(model, varλ, varP, warmstart)
|
||||||
Base.Libc.flush_cstdio()
|
Base.Libc.flush_cstdio()
|
||||||
@ -224,7 +224,7 @@ function fillfrominternal!(m::JuMP.Model, traits)
|
|||||||
m.objBound = NaN
|
m.objBound = NaN
|
||||||
m.objVal = NaN
|
m.objVal = NaN
|
||||||
m.colVal = fill(NaN, numCols)
|
m.colVal = fill(NaN, numCols)
|
||||||
m.linconstrDuals = Array{Float64}(0)
|
m.linconstrDuals = Array{Float64}(undef, 0)
|
||||||
|
|
||||||
discrete = (traits.int || traits.sos)
|
discrete = (traits.int || traits.sos)
|
||||||
|
|
||||||
@ -262,7 +262,7 @@ function fillfrominternal!(m::JuMP.Model, traits)
|
|||||||
@assert length(infray) == numRows
|
@assert length(infray) == numRows
|
||||||
infray
|
infray
|
||||||
catch
|
catch
|
||||||
suppress_warnings || warn("Infeasibility ray (Farkas proof) not available")
|
@warn("Infeasibility ray (Farkas proof) not available")
|
||||||
fill(NaN, numRows)
|
fill(NaN, numRows)
|
||||||
end
|
end
|
||||||
elseif stat == :Unbounded
|
elseif stat == :Unbounded
|
||||||
@ -271,7 +271,7 @@ function fillfrominternal!(m::JuMP.Model, traits)
|
|||||||
@assert length(unbdray) == numCols
|
@assert length(unbdray) == numCols
|
||||||
unbdray
|
unbdray
|
||||||
catch
|
catch
|
||||||
suppress_warnings || warn("Unbounded ray not available")
|
@warn("Unbounded ray not available")
|
||||||
fill(NaN, numCols)
|
fill(NaN, numCols)
|
||||||
end
|
end
|
||||||
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
|
# Do a separate try since getobjval could work while getobjbound does not and vice versa
|
||||||
objBound = MathProgBase.getobjbound(m.internalModel) + m.obj.aff.constant
|
objBound = MathProgBase.getobjbound(m.internalModel) + m.obj.aff.constant
|
||||||
m.objBound = objBound
|
m.objBound = objBound
|
||||||
|
catch
|
||||||
|
@warn("objBound could not be obtained")
|
||||||
end
|
end
|
||||||
|
|
||||||
try
|
try
|
||||||
objVal = MathProgBase.getobjval(m.internalModel) + m.obj.aff.constant
|
objVal = MathProgBase.getobjval(m.internalModel) + m.obj.aff.constant
|
||||||
colVal = MathProgBase.getsolution(m.internalModel)[1:numCols]
|
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
|
# Don't corrupt the answers if one of the above two calls fails
|
||||||
m.objVal = objVal
|
m.objVal = objVal
|
||||||
m.colVal = colVal
|
m.colVal = colVal
|
||||||
|
catch
|
||||||
|
@warn("objVal/colVal could not be obtained")
|
||||||
end
|
end
|
||||||
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
|
return stat
|
||||||
end
|
end
|
||||||
|
Loading…
Reference in New Issue
Block a user