From dbca3945a83367d8df4868553517fa59ad2c43a0 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Tue, 19 Sep 2017 17:37:35 +0200 Subject: [PATCH 01/55] call flush_cstdio() AFTER redirect_stdout to fix SCS log buffering issue --- src/SDPs.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/SDPs.jl b/src/SDPs.jl index ab1cffb..b2a9fa5 100644 --- a/src/SDPs.jl +++ b/src/SDPs.jl @@ -58,11 +58,8 @@ end function solve_SDP(SDP_problem) info(logger, Base.repr(SDP_problem)) - # to change buffering mode of stdout to _IOLBF (line bufferin) - # see https://github.com/JuliaLang/julia/issues/8765 - ccall((:printf, "libc"), Int, (Ptr{UInt8},), "\n"); - o = redirect_stdout(solver_logger.handlers["solver_log"].io) + Base.Libc.flush_cstdio() t = @timed solution_status = JuMP.solve(SDP_problem) info(logger, timed_msg(t)) From 344e11a9744bda67981b4f565e91f90096ad082d Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 9 Oct 2017 18:06:29 +0200 Subject: [PATCH 02/55] overwrite Q to allow gc --- src/CheckSolution.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index b80756f..7bcb74a 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -116,20 +116,20 @@ end function rationalize_and_project{T}(Q::AbstractArray{T}, δ::T, logger) info(logger, "") info(logger, "Rationalizing with accuracy $δ") - t = @timed Q_ℚ = ℚ(Q, δ) + t = @timed Q = ℚ(Q, δ) info(logger, timed_msg(t)) info(logger, "Projecting columns of the rationalized Q to the augmentation ideal...") - t = @timed Q_int = correct_to_augmentation_ideal(Q_ℚ) + t = @timed Q = correct_to_augmentation_ideal(Q) info(logger, timed_msg(t)) info(logger, "Checking that sum of every column contains 0.0... ") - check = all([0.0 in sum(view(Q_int, :, i)) for i in 1:size(Q_int, 2)]) + check = all([0.0 in sum(view(Q, :, i)) for i in 1:size(Q, 2)]) info(logger, (check? "They do." : "FAILED!")) @assert check - return Q_int + return Q end function check_distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen; From d7bc94b64cffc4eee2c74f869dc6b85fb076ae10 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 9 Oct 2017 18:08:22 +0200 Subject: [PATCH 03/55] @parallel compute_SOS --- src/CheckSolution.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index 7bcb74a..0fb7610 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -36,17 +36,17 @@ function compute_SOS(sqrt_matrix, elt::GroupRingElem) l = length(elt.coeffs) pm = parent(elt).pm - result = zeros(eltype(sqrt_matrix), l) - for i in 1:n - result .+= groupring_square(view(sqrt_matrix,:,i), l, pm) - end - - # @everywhere groupring_square = PropertyT.groupring_square - # - # result = @parallel (+) for i in 1:n - # groupring_square(view(sqrt_matrix,:,i), length(elt.coeffs), parent(elt).pm) + # result = zeros(eltype(sqrt_matrix), l) + # for i in 1:n + # result .+= groupring_square(view(sqrt_matrix,:,i), l, pm) # end + @everywhere groupring_square = PropertyT.groupring_square + + result = @parallel (+) for i in 1:n + groupring_square(view(sqrt_matrix,:,i), length(elt.coeffs), parent(elt).pm) + end + return GroupRingElem(result, parent(elt)) end From 6070bbf709a9aba7cd63ad2e427a7c3f2777b38d Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 9 Oct 2017 19:44:54 +0200 Subject: [PATCH 04/55] rename name -> prefix in filenames functions --- src/PropertyT.jl | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 6c5b3f8..d55bc40 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -30,29 +30,23 @@ function exists(fname::String) return isfile(fname) || islink(fname) end -function pmΔfilenames(name::String) - if !isdir(name) - mkdir(name) - end - prefix = name - pm_filename = joinpath(prefix, "pm.jld") - Δ_coeff_filename = joinpath(prefix, "delta.jld") - return pm_filename, Δ_coeff_filename +function pmΔfilenames(prefix::String) + isdir(prefix) || mkdir(prefix) + pm_filename = joinpath(prefix, "pm.jld") + Δ_coeff_filename = joinpath(prefix, "delta.jld") + return pm_filename, Δ_coeff_filename end -function λSDPfilenames(name::String) - if !isdir(name) - mkdir(name) - end - prefix = name - λ_filename = joinpath(prefix, "lambda.jld") - SDP_filename = joinpath(prefix, "SDPmatrix.jld") - return λ_filename, SDP_filename +function λSDPfilenames(prefix::String) + isdir(prefix) || mkdir(prefix) + λ_filename = joinpath(prefix, "lambda.jld") + SDP_filename = joinpath(prefix, "SDPmatrix.jld") + return λ_filename, SDP_filename end -function ΔandSDPconstraints(name::String, G::Group) +function ΔandSDPconstraints(prefix::String, G::Group) info(logger, "Loading precomputed pm, Δ, sdp_constraints...") - pm_fname, Δ_fname = pmΔfilenames(name) + pm_fname, Δ_fname = pmΔfilenames(prefix) product_matrix = load(pm_fname, "pm") sdp_constraints = constraints_from_pm(product_matrix) From 36a1151a739689a68db58a62b67ed0002ceb5c08 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 9 Oct 2017 19:45:28 +0200 Subject: [PATCH 05/55] add functions to produce prepath, fullpath out of Settings --- src/Orbit-wise.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index fca5438..377cf82 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -15,6 +15,11 @@ immutable Settings tol::Float64 end +prefix(s::Settings) = s.name +suffix(s::Settings) = "$(s.upper_bound)" +prepath(s::Settings) = prefix(s) +fullpath(s::Settings) = joinpath(prefix(s), suffix(s)) + immutable OrbitData name::String Us::Vector From b3f1961fbf2e7fb2543f5946b45f7d30247da21d Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 9 Oct 2017 19:46:08 +0200 Subject: [PATCH 06/55] use prepath/fullpath(::Settings) in OrbitData and in create_SDP_problem generation --- src/Orbit-wise.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 377cf82..7ba0ec6 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -30,27 +30,27 @@ immutable OrbitData dims::Vector{Int} end -function OrbitData(name::String) - splap = load(joinpath(name, "delta.jld"), "Δ"); - pm = load(joinpath(name, "pm.jld"), "pm"); +function OrbitData(sett::Settings) + splap = load(joinpath(prepath(sett), "delta.jld"), "Δ"); + pm = load(joinpath(prepath(sett), "pm.jld"), "pm"); cnstr = PropertyT.constraints_from_pm(pm); splap² = similar(splap) splap² = GroupRings.mul!(splap², splap, splap, pm); # Uπs = load(joinpath(name, "U_pis.jld"), "Uπs"); - Uπs = load(joinpath(name, "U_pis.jld"), "spUπs"); + Uπs = load(joinpath(prepath(sett), "U_pis.jld"), "spUπs"); #dimensions of the corresponding πs: - dims = load(joinpath(name, "U_pis.jld"), "dims") + dims = load(joinpath(prepath(sett), "U_pis.jld"), "dims") m, P = init_model(Uπs); - orbits = load(joinpath(name, "orbits.jld"), "orbits"); + orbits = load(joinpath(prepath(sett), "orbits.jld"), "orbits"); n = size(Uπs[1],1) orb_spcnstrm = [orbit_constraint(cnstr[collect(orb)], n) for orb in orbits] orb_splap = orbit_spvector(splap, orbits) orb_splap² = orbit_spvector(splap², orbits) - orbData = OrbitData(name, Uπs, P, orb_spcnstrm, orb_splap, orb_splap², dims); + orbData = OrbitData(fullpath(sett), Uπs, P, orb_spcnstrm, orb_splap, orb_splap², dims); # orbData = OrbitData(name, Uπs, P, orb_spcnstrm, splap, splap², dims); @@ -175,14 +175,14 @@ function init_model(Uπs) return m, P end -function create_SDP_problem(name::String; upper_bound=Inf) +function create_SDP_problem(sett::Settings) info(logger, "Loading orbit data....") - t = @timed SDP_problem, orb_data = OrbitData(name); + t = @timed SDP_problem, orb_data = OrbitData(sett); info(logger, PropertyT.timed_msg(t)) - if upper_bound < Inf + if sett.upper_bound < Inf λ = JuMP.getvariable(SDP_problem, :λ) - JuMP.@constraint(SDP_problem, λ <= upper_bound) + JuMP.@constraint(SDP_problem, λ <= sett.upper_bound) end t = length(orb_data.laplacian) @@ -227,7 +227,7 @@ function check_property_T(sett::Settings) λ, P = PropertyT.λandP(sett.name) else info(logger, "Creating SDP problem...") - SDP_problem, orb_data = create_SDP_problem(sett.name, upper_bound=sett.upper_bound) + SDP_problem, orb_data = create_SDP_problem(sett) JuMP.setsolver(SDP_problem, sett.solver) λ, P = λandP(SDP_problem, orb_data, sett) From bc2a827342c71d6cc5da0227cf437bf928ec2442 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 9 Oct 2017 19:49:49 +0200 Subject: [PATCH 07/55] use prepath(sett)/fullpath(sett) instead of sett.name --- src/Orbit-wise.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 7ba0ec6..fa612bc 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -101,12 +101,12 @@ sparsify{T}(U::AbstractArray{T}, tol=eps(T)) = sparsify!(deepcopy(U), tol) function init_orbit_data(logger, sett::Settings; radius=2) - ex(fname) = isfile(joinpath(sett.name, fname)) + ex(fname) = isfile(joinpath(prepath(sett), fname)) files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld"]) if !all(files_exists) - compute_orbit_data(logger, sett.name, sett.G, sett.S, sett.AutS, radius=radius) + compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.AutS, radius=radius) end return 0 @@ -212,7 +212,7 @@ function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) t = @timed recP = reconstruct_sol(preps, data.Us, Ps, data.dims) info(logger, PropertyT.timed_msg(t)) - fname = PropertyT.λSDPfilenames(data.name)[2] + fname = PropertyT.λSDPfilenames(fullpath(sett))[2] save(fname, "origP", Ps, "P", recP) return λ, recP end @@ -221,10 +221,10 @@ function check_property_T(sett::Settings) init_orbit_data(logger, sett, radius=sett.radius) - fnames = PropertyT.λSDPfilenames(sett.name) + fnames = PropertyT.λSDPfilenames(fullpath(sett)) if all(isfile.(fnames)) - λ, P = PropertyT.λandP(sett.name) + λ, P = PropertyT.λandP(fullpath(sett)) else info(logger, "Creating SDP problem...") SDP_problem, orb_data = create_SDP_problem(sett) @@ -239,9 +239,9 @@ function check_property_T(sett::Settings) info(logger, "minimum(P) = $(minimum(P))") if λ > 0 - pm_fname = joinpath(sett.name, "pm.jld") + pm_fname = joinpath(prepath(sett), "pm.jld") RG = GroupRing(sett.G, load(pm_fname, "pm")) - Δ_fname = joinpath(sett.name, "delta.jld") + Δ_fname = joinpath(prepath(sett), "delta.jld") Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || From 6cbf5e694462ad9e29f3af14e5c1bc1a4ab8d898 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 14:27:26 +0200 Subject: [PATCH 08/55] Allow for AbstractMathProfSolver in Settings --- src/Orbit-wise.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index fa612bc..9d7ef7c 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -10,7 +10,7 @@ immutable Settings S::Vector AutS::Group radius::Int - solver::SCSSolver + solver::AbstractMathProgSolver upper_bound::Float64 tol::Float64 end @@ -166,8 +166,10 @@ function init_model(Uπs) for k in 1:l s = size(Uπs[k],2) - P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) - JuMP.@SDconstraint(m, P[k] >= 0.0) + if s > 0 + P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) + JuMP.@SDconstraint(m, P[k] >= 0.0) + end end JuMP.@variable(m, λ >= 0.0) From 46a915026616bf5be6a048db063e6da73b2a1f9b Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 14:28:01 +0200 Subject: [PATCH 09/55] solver is not defined here --- src/PropertyT.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index d55bc40..35937e8 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -142,7 +142,7 @@ function λandP(name::String, SDP_problem::JuMP.Model, varλ, varP) save(λ_fname, "λ", λ) save(P_fname, "P", P) else - throw(ErrorException("Solver $solver did not produce a valid solution!: λ = $λ")) + throw(ErrorException("Solver did not produce a valid solution!: λ = $λ")) end return λ, P From 32c020fd8816fe0198a013bb242e1aad6c78bf1b Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 14:28:38 +0200 Subject: [PATCH 10/55] make initial P nullable --- src/PropertyT.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 35937e8..101c5da 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -150,7 +150,7 @@ end function compute_λandP(m, varλ, varP) λ = 0.0 - P = nothing + P = Vector{Nullable{Array{Float64,2}}} while λ == 0.0 try solve_SDP(m) From 94c9a75a2187495ce65adda5e63a09a11da0da7d Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 14:29:47 +0200 Subject: [PATCH 11/55] we don't need to pass id element -- it could be derived from S --- src/SDPs.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/SDPs.jl b/src/SDPs.jl index b2a9fa5..7a41479 100644 --- a/src/SDPs.jl +++ b/src/SDPs.jl @@ -13,18 +13,18 @@ function constraints_from_pm(pm, total_length=maximum(pm)) return constraints end -function splaplacian{TT<:Group}(RG::GroupRing{TT}, S, Id=RG.group(), T::Type=Float64) +function splaplacian(RG::GroupRing, S, T::Type=Float64) result = RG(T) - result[Id] = T(length(S)) + result[RG.group()] = T(length(S)) for s in S result[s] -= one(T) end return result end -function splaplacian{TT<:Ring}(RG::GroupRing{TT}, S, Id=one(RG.group), T::Type=Float64) +function splaplacian{TT<:Ring}(RG::GroupRing{TT}, S, T::Type=Float64) result = RG(T) - result[Id] = T(length(S)) + result[one(RG.group)] = T(length(S)) for s in S result[s] -= one(T) end From 01e6625571c658056f2d1937dfc404cda3e58dfc Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:19:53 +0200 Subject: [PATCH 12/55] Revert "Allow for AbstractMathProfSolver in Settings" This reverts commit 6cbf5e694462ad9e29f3af14e5c1bc1a4ab8d898. --- src/Orbit-wise.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 9d7ef7c..fa612bc 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -10,7 +10,7 @@ immutable Settings S::Vector AutS::Group radius::Int - solver::AbstractMathProgSolver + solver::SCSSolver upper_bound::Float64 tol::Float64 end @@ -166,10 +166,8 @@ function init_model(Uπs) for k in 1:l s = size(Uπs[k],2) - if s > 0 - P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) - JuMP.@SDconstraint(m, P[k] >= 0.0) - end + P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) + JuMP.@SDconstraint(m, P[k] >= 0.0) end JuMP.@variable(m, λ >= 0.0) From b96b39e910a3afa5196c523a99d0b41950fbcfa5 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:20:28 +0200 Subject: [PATCH 13/55] Revert "make initial P nullable" This reverts commit 32c020fd8816fe0198a013bb242e1aad6c78bf1b. --- src/PropertyT.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 101c5da..35937e8 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -150,7 +150,7 @@ end function compute_λandP(m, varλ, varP) λ = 0.0 - P = Vector{Nullable{Array{Float64,2}}} + P = nothing while λ == 0.0 try solve_SDP(m) From d8b2f0ab3063ea474d706710240f822d6bee2a7a Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:25:13 +0200 Subject: [PATCH 14/55] increase sigfigs to 12; remove explicit precision setting (it's 53 anyway); --- src/CheckSolution.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index 0fb7610..8ddc3c6 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -3,8 +3,7 @@ import Base: rationalize using IntervalArithmetic IntervalArithmetic.setrounding(Interval, :tight) -IntervalArithmetic.setformat(sigfigs=10) -IntervalArithmetic.setprecision(Interval, 53) # slightly faster than 256 +IntervalArithmetic.setformat(sigfigs=12) import IntervalArithmetic.± From 8defd71b596e2f3c262a46bbc6dd175a21a6c87d Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:29:19 +0200 Subject: [PATCH 15/55] new low-level compute_SOS function --- src/CheckSolution.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index 8ddc3c6..0da910d 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -30,23 +30,22 @@ function groupring_square(vect::AbstractVector, l, pm) return GroupRings.mul!(similar(zzz), zzz, zzz, pm) end -function compute_SOS(sqrt_matrix, elt::GroupRingElem) - n = size(sqrt_matrix,2) - l = length(elt.coeffs) - pm = parent(elt).pm +function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int) + n = size(Q,2) - # result = zeros(eltype(sqrt_matrix), l) + # result = zeros(eltype(Q), l) # for i in 1:n - # result .+= groupring_square(view(sqrt_matrix,:,i), l, pm) + # result .+= groupring_square(view(Q,:,i), l, pm) # end @everywhere groupring_square = PropertyT.groupring_square result = @parallel (+) for i in 1:n - groupring_square(view(sqrt_matrix,:,i), length(elt.coeffs), parent(elt).pm) + print(" $i") + groupring_square(view(Q,:,i), l, pm) end - return GroupRingElem(result, parent(elt)) + return result end function correct_to_augmentation_ideal{T<:Rational}(sqrt_matrix::Array{T,2}) From 0225f38066bfd026dcac1b236bdea38ef7f43a6a Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:30:59 +0200 Subject: [PATCH 16/55] just two simple distance_to_cone functions for Numbers and for Intervals --- src/CheckSolution.jl | 57 ++++++++++++++------------------------------ 1 file changed, 18 insertions(+), 39 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index 0da910d..c7d1dee 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -60,55 +60,34 @@ function correct_to_augmentation_ideal{T<:Rational}(sqrt_matrix::Array{T,2}) return sqrt_corrected end -function distance_to_cone{T<:Rational}(λ::T, sqrt_matrix::Array{T,2}, Δ::GroupRingElem{T}, wlen) - SOS = compute_SOS(sqrt_matrix, Δ) +function distance_to_cone{S<:Interval}(elt::GroupRingElem, Q::AbstractArray{S,2}, wlen::Int) + SOS = compute_SOS(Q, parent(elt), length(elt.coeffs)) + SOS_diff = elt - SOS - SOS_diff = EOI(Δ, λ) - SOS - eoi_SOS_L1_dist = norm(SOS_diff,1) + ɛ_dist = GroupRings.augmentation(SOS_diff) + info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)") - info(logger, "λ = $λ (≈$(@sprintf("%.10f", float(λ)))") - ɛ_dist = GroupRings.augmentation(SOS_diff) - if ɛ_dist ≠ 0//1 - warn(logger, "The SOS is not in the augmentation ideal, numbers below are meaningless!") - end - info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) = $ɛ_dist") - info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ = $(@sprintf("%.10f", float(eoi_SOS_L1_dist)))") + eoi_SOS_L1_dist = norm(SOS_diff,1) + info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)") - distance_to_cone = λ - 2^(wlen-1)*eoi_SOS_L1_dist - return distance_to_cone + dist = 2^(wlen-1)*eoi_SOS_L1_dist + return dist end -function distance_to_cone{T<:Rational, S<:Interval}(λ::T, sqrt_matrix::AbstractArray{S,2}, Δ::GroupRingElem{T}, wlen) - SOS = compute_SOS(sqrt_matrix, Δ) - info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(GroupRings.augmentation(SOS))") - λ_int = @interval(λ) - Δ_int = GroupRingElem([@interval(c) for c in Δ.coeffs], parent(Δ)) - SOS_diff = EOI(Δ_int, λ_int) - SOS - eoi_SOS_L1_dist = norm(SOS_diff,1) +function distance_to_cone{T}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int) + SOS = compute_SOS(Q, parent(elt), length(elt.coeffs)) + SOS_diff = elt - SOS - info(logger, "λ = $λ (≈≥$(@sprintf("%.10f",float(λ))))") - ɛ_dist = GroupRings.augmentation(SOS_diff) + ɛ_dist = GroupRings.augmentation(SOS_diff) + info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))") - info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)") - info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)") + eoi_SOS_L1_dist = norm(SOS_diff,1) + info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L1_dist))") - distance_to_cone = λ - 2^(wlen-1)*eoi_SOS_L1_dist - return distance_to_cone + dist = 2^(wlen-1)*eoi_SOS_L1_dist + return dist end -function distance_to_cone(λ, sqrt_matrix::AbstractArray, Δ::GroupRingElem, wlen) - SOS = compute_SOS(sqrt_matrix, Δ) - - SOS_diff = EOI(Δ, λ) - SOS - eoi_SOS_L1_dist = norm(SOS_diff,1) - - info(logger, "λ = $λ") - ɛ_dist = GroupRings.augmentation(SOS_diff) - info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))") - info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L1_dist))") - - distance_to_cone = λ - 2^(wlen-1)*eoi_SOS_L1_dist - return distance_to_cone end function rationalize_and_project{T}(Q::AbstractArray{T}, δ::T, logger) From 96ce02852f57a2784ca6043cf13fd0f981e66024 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:32:58 +0200 Subject: [PATCH 17/55] replace rationalise&project with interval&project --- src/CheckSolution.jl | 32 +++++++++++++------------------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index c7d1dee..4d0668e 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -48,16 +48,6 @@ function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int) return result end -function correct_to_augmentation_ideal{T<:Rational}(sqrt_matrix::Array{T,2}) - l = size(sqrt_matrix, 2) - sqrt_corrected = Array{Interval{Float64}}(l,l) - Threads.@threads for j in 1:l - col = sum(view(sqrt_matrix, :,j))//l - for i in 1:l - sqrt_corrected[i,j] = (Float64(sqrt_matrix[i,j]) - Float64(col)) ± eps(0.0) - end - end - return sqrt_corrected end function distance_to_cone{S<:Interval}(elt::GroupRingElem, Q::AbstractArray{S,2}, wlen::Int) @@ -88,17 +78,21 @@ function distance_to_cone{T}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::In return dist end +function augIdproj{T, I<:AbstractInterval}(S::Type{I}, Q::AbstractArray{T,2}) + l = size(Q, 2) + R = zeros(Interval, (l,l)) + Threads.@threads for j in 1:l + col = sum(view(Q, :,j))/l + for i in 1:l + R[i,j] = Q[i,j] - col ± eps(0.0) + end + end + return R end -function rationalize_and_project{T}(Q::AbstractArray{T}, δ::T, logger) - info(logger, "") - info(logger, "Rationalizing with accuracy $δ") - t = @timed Q = ℚ(Q, δ) - info(logger, timed_msg(t)) - - info(logger, "Projecting columns of the rationalized Q to the augmentation ideal...") - t = @timed Q = correct_to_augmentation_ideal(Q) - info(logger, timed_msg(t)) +function augIdproj{T}(Q::AbstractArray{T,2}, logger) + info(logger, "Projecting columns of Q to the augmentation ideal...") + @logtime logger Q = augIdproj(Interval, Q) info(logger, "Checking that sum of every column contains 0.0... ") check = all([0.0 in sum(view(Q, :, i)) for i in 1:size(Q, 2)]) From 4ae2c617cb9568c581d62d8ae0b0415ef8f2d974 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:33:18 +0200 Subject: [PATCH 18/55] Higher level compute_SOS --- src/CheckSolution.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index 4d0668e..f4a975a 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -48,6 +48,9 @@ function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int) return result end +function compute_SOS(Q::AbstractArray, RG::GroupRing, l::Int) + result = compute_SOS(Q, RG.pm, l) + return GroupRingElem(result, RG) end function distance_to_cone{S<:Interval}(elt::GroupRingElem, Q::AbstractArray{S,2}, wlen::Int) From 098a7ce948edf8a55aae7d7f2d2292301703815f Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:34:07 +0200 Subject: [PATCH 19/55] check_distance_to... -> distance_to_... and a lot of old code removal --- src/CheckSolution.jl | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index f4a975a..df0bb3f 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -106,13 +106,12 @@ function augIdproj{T}(Q::AbstractArray{T,2}, logger) return Q end -function check_distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen; - tol=1e-14, rational=false) - +function distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen::Int) info(logger, "------------------------------------------------------------") - info(logger, "") + info(logger, "λ = $λ") info(logger, "Checking in floating-point arithmetic...") - t = @timed fp_distance = distance_to_cone(λ, Q, Δ, wlen) + Δ²_λΔ = EOI(Δ, λ) + t = @timed fp_distance = λ - distance_to_cone(Δ²_λΔ, Q, wlen) info(logger, timed_msg(t)) info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))") info(logger, "------------------------------------------------------------") @@ -122,26 +121,17 @@ function check_distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen; end info(logger, "") - Q_ℚω_int = rationalize_and_project(Q, tol, logger) - λ_ℚ = ℚ(λ, tol) - Δ_ℚ = ℚ(Δ, tol) + Q = augIdproj(Q, logger) info(logger, "Checking in interval arithmetic") + λ = @interval(λ) + Δ = GroupRingElem([@interval(c) for c in Δ.coeffs], parent(Δ)) + Δ²_λΔ = EOI(Δ, λ) - t = @timed Interval_dist_to_ΣSq = distance_to_cone(λ_ℚ, Q_ℚω_int, Δ_ℚ, wlen) + t = @timed Interval_dist_to_ΣSq = λ - distance_to_cone(Δ²_λΔ, Q, wlen) info(logger, timed_msg(t)) info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_ΣSq)") info(logger, "------------------------------------------------------------") - if Interval_dist_to_ΣSq.lo ≤ 0 || !rational - return Interval_dist_to_ΣSq - else - info(logger, "Checking Projected SOS decomposition in exact rational arithmetic...") - t = @timed ℚ_dist_to_ΣSq = distance_to_cone(λ_ℚ, Q_ℚω, Δ_ℚ, wlen) - info(logger, timed_msg(t)) - @assert isa(ℚ_dist_to_ΣSq, Rational) - info(logger, "Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_ΣSq,8)))") - info(logger, "------------------------------------------------------------") - return ℚ_dist_to_ΣSq - end + return Interval_dist_to_ΣSq end From f206a1398056a21950642ed701afe2bd75428c5e Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:35:50 +0200 Subject: [PATCH 20/55] update function names --- src/Orbit-wise.jl | 2 +- src/PropertyT.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index fa612bc..5d81173 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -249,7 +249,7 @@ function check_property_T(sett::Settings) # @assert P == Symmetric(P) Q = real(sqrtm(Symmetric(P))) - sgap = PropertyT.check_distance_to_positive_cone(Δ, λ, Q, 2*sett.radius, tol=sett.tol, rational=false) + sgap = distance_to_positive_cone(Δ, λ, Q, 2*sett.radius) if isa(sgap, Interval) sgap = sgap.lo end diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 35937e8..7368f6a 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -208,7 +208,7 @@ function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius) # @assert P == Symmetric(P) Q = real(sqrtm(Symmetric(P))) - sgap = check_distance_to_positive_cone(Δ, λ, Q, 2*radius, tol=tol) + sgap = distance_to_positive_cone(Δ, λ, Q, 2*radius) if isa(sgap, Interval) sgap = sgap.lo end From ca18f75efdd742603db9dcd188221bdf48acee1b Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:36:32 +0200 Subject: [PATCH 21/55] clean the main check_property_T in both cases --- src/Orbit-wise.jl | 7 ++----- src/PropertyT.jl | 10 ++++------ 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 5d81173..670c329 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -221,9 +221,7 @@ function check_property_T(sett::Settings) init_orbit_data(logger, sett, radius=sett.radius) - fnames = PropertyT.λSDPfilenames(fullpath(sett)) - - if all(isfile.(fnames)) + if all(isfile.(λSDPfilenames(fullpath(sett)))) λ, P = PropertyT.λandP(fullpath(sett)) else info(logger, "Creating SDP problem...") @@ -239,9 +237,8 @@ function check_property_T(sett::Settings) info(logger, "minimum(P) = $(minimum(P))") if λ > 0 - pm_fname = joinpath(prepath(sett), "pm.jld") + pm_fname, Δ_fname = pmΔfilenames(prepath(sett)) RG = GroupRing(sett.G, load(pm_fname, "pm")) - Δ_fname = joinpath(prepath(sett), "delta.jld") Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 7368f6a..b27d3cb 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -182,15 +182,10 @@ function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius) info(logger, "|R[G]|.pm = $(size(parent(Δ).pm))") if all(exists.(λSDPfilenames(name))) - # cached λ, P = λandP(name) else - # compute info(logger, "Creating SDP problem...") - - t = @timed SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound) - info(logger, timed_msg(t)) - + SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound) JuMP.setsolver(SDP_problem, solver) λ, P = λandP(name, SDP_problem, λ, P) @@ -202,6 +197,9 @@ function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius) info(logger, "minimum(P) = $(minimum(P))") if λ > 0 + pm_fname, Δ_fname = pmΔfilenames(name) + RG = GroupRing(parent(first(S)), load(pm_fname, "pm")) + Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) isapprox(eigvals(P), abs(eigvals(P)), atol=tol) || warn("The solution matrix doesn't seem to be positive definite!") From 5bb311141d80c3c268e95f848509fcd72e62af1d Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 27 Oct 2017 18:36:59 +0200 Subject: [PATCH 22/55] splaplacian doesn't have Id argument --- src/PropertyT.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index b27d3cb..abecd44 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -80,7 +80,7 @@ function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2) RG = GroupRing(parent(Id), B, pm) - Δ = splaplacian(RG, S, Id) + Δ = splaplacian(RG, S) return Δ, sdp_constraints end From 7fa49756ba24156a556be37b7246a8e5a5071daa Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 2 Nov 2017 13:44:08 +0100 Subject: [PATCH 23/55] @logtime macro and time_string function --- src/PropertyT.jl | 52 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 39 insertions(+), 13 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index abecd44..0fee9f0 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -26,6 +26,45 @@ function setup_logging(name::String) return logger end +macro logtime(logger, ex) + quote + local stats = Base.gc_num() + local elapsedtime = Base.time_ns() + local val = $(esc(ex)) + elapsedtime = Base.time_ns() - elapsedtime + local diff = Base.GC_Diff(Base.gc_num(), stats) + local ts = time_string(elapsedtime, diff.allocd, diff.total_time, + Base.gc_alloc_count(diff)) + esc(info(logger, ts)) + val + end +end + +function time_string(elapsedtime, bytes, gctime, allocs) + str = @sprintf("%10.6f seconds", elapsedtime/1e9) + if bytes != 0 || allocs != 0 + bytes, mb = Base.prettyprint_getunits(bytes, length(Base._mem_units), Int64(1024)) + allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000)) + if ma == 1 + str*= @sprintf(" (%d%s allocation%s: ", allocs, Base._cnt_units[ma], allocs==1 ? "" : "s") + else + str*= @sprintf(" (%.2f%s allocations: ", allocs, Base._cnt_units[ma]) + end + if mb == 1 + str*= @sprintf("%d %s%s", bytes, Base._mem_units[mb], bytes==1 ? "" : "s") + else + str*= @sprintf("%.3f %s", bytes, Base._mem_units[mb]) + end + if gctime > 0 + str*= @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime) + end + str*=")" + elseif gctime > 0 + str*= @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime) + end + return str +end + function exists(fname::String) return isfile(fname) || islink(fname) end @@ -84,19 +123,6 @@ function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2) return Δ, sdp_constraints end -macro logtime(logger, ex) - quote - local stats = Base.gc_num() - local elapsedtime = Base.time_ns() - local val = $(esc(ex)) - elapsedtime = Base.time_ns() - elapsedtime - local diff = Base.GC_Diff(Base.gc_num(), stats) - local ts = time_string(elapsedtime, diff.allocd, diff.total_time, - Base.gc_alloc_count(diff)) - esc(warn($(esc(logger)), ts)) - val - end -end function timed_msg(t) elapsed = t[2] From 3eab5b721a15eb2bddd7db52b748e3491f0b8849 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 2 Nov 2017 13:45:12 +0100 Subject: [PATCH 24/55] replace @timed by @logtime --- src/CheckSolution.jl | 6 ++---- src/Orbit-wise.jl | 12 ++++-------- src/PropertyT.jl | 6 ++---- src/SDPs.jl | 3 +-- 4 files changed, 9 insertions(+), 18 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index df0bb3f..b38f838 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -111,8 +111,7 @@ function distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen::Int) info(logger, "λ = $λ") info(logger, "Checking in floating-point arithmetic...") Δ²_λΔ = EOI(Δ, λ) - t = @timed fp_distance = λ - distance_to_cone(Δ²_λΔ, Q, wlen) - info(logger, timed_msg(t)) + @logtime logger fp_distance = λ - distance_to_cone(Δ²_λΔ, Q, wlen) info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))") info(logger, "------------------------------------------------------------") @@ -128,8 +127,7 @@ function distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen::Int) Δ = GroupRingElem([@interval(c) for c in Δ.coeffs], parent(Δ)) Δ²_λΔ = EOI(Δ, λ) - t = @timed Interval_dist_to_ΣSq = λ - distance_to_cone(Δ²_λΔ, Q, wlen) - info(logger, timed_msg(t)) + @logtime logger Interval_dist_to_ΣSq = λ - distance_to_cone(Δ²_λΔ, Q, wlen) info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_ΣSq)") info(logger, "------------------------------------------------------------") diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 670c329..3846310 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -177,8 +177,7 @@ end function create_SDP_problem(sett::Settings) info(logger, "Loading orbit data....") - t = @timed SDP_problem, orb_data = OrbitData(sett); - info(logger, PropertyT.timed_msg(t)) + @logtime logger SDP_problem, orb_data = OrbitData(sett); if sett.upper_bound < Inf λ = JuMP.getvariable(SDP_problem, :λ) @@ -187,8 +186,7 @@ function create_SDP_problem(sett::Settings) t = length(orb_data.laplacian) info(logger, "Adding $t constraints ... ") - t = @timed addconstraints!(SDP_problem, orb_data) - info(logger, PropertyT.timed_msg(t)) + @logtime logger addconstraints!(SDP_problem, orb_data) return SDP_problem, orb_data end @@ -206,11 +204,9 @@ function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) info(logger, "Reconstructing P...") - t = @timed preps = perm_reps(sett.G, sett.S, sett.AutS, sett.radius) - info(logger, PropertyT.timed_msg(t)) + @logtime logger preps = perm_reps(sett.G, sett.S, sett.AutS, sett.radius) - t = @timed recP = reconstruct_sol(preps, data.Us, Ps, data.dims) - info(logger, PropertyT.timed_msg(t)) + @logtime logger recP = reconstruct_sol(preps, data.Us, Ps, data.dims) fname = PropertyT.λSDPfilenames(fullpath(sett))[2] save(fname, "origP", Ps, "P", recP) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 0fee9f0..4e04287 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -110,12 +110,10 @@ function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2) info(logger, "Generated balls of sizes $sizes") info(logger, "Creating product matrix...") - t = @timed pm = GroupRings.create_pm(B, GroupRings.reverse_dict(B), sizes[radius]; twisted=true) - info(logger, timed_msg(t)) + @logtime logger pm = GroupRings.create_pm(B, GroupRings.reverse_dict(B), sizes[radius]; twisted=true) info(logger, "Creating sdp_constratints...") - t = @timed sdp_constraints = PropertyT.constraints_from_pm(pm) - info(logger, timed_msg(t)) + @logtime logger sdp_constraints = PropertyT.constraints_from_pm(pm) RG = GroupRing(parent(Id), B, pm) diff --git a/src/SDPs.jl b/src/SDPs.jl index 7a41479..972f3b5 100644 --- a/src/SDPs.jl +++ b/src/SDPs.jl @@ -61,8 +61,7 @@ function solve_SDP(SDP_problem) o = redirect_stdout(solver_logger.handlers["solver_log"].io) Base.Libc.flush_cstdio() - t = @timed solution_status = JuMP.solve(SDP_problem) - info(logger, timed_msg(t)) + @logtime logger solution_status = JuMP.solve(SDP_problem) Base.Libc.flush_cstdio() redirect_stdout(o) From 1c67cdbcea8d9554a05f7fa7d349438c97a8dfab Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 2 Nov 2017 13:45:27 +0100 Subject: [PATCH 25/55] remove timed_msg --- src/PropertyT.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 4e04287..cc117c2 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -121,16 +121,6 @@ function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2) return Δ, sdp_constraints end - -function timed_msg(t) - elapsed = t[2] - bytes_alloc = t[3] - gc_time = t[4] - gc_diff = t[5] - - return "took: $elapsed s, allocated: $bytes_alloc bytes ($(gc_diff.poolalloc) allocations)." -end - function λandP(name::String) λ_fname, SDP_fname = λSDPfilenames(name) f₁ = exists(λ_fname) From 79d5c7f933f8d6d1e571ff2dbbf0f228e734e4e8 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 2 Nov 2017 13:48:20 +0100 Subject: [PATCH 26/55] rename parameter S -> T --- src/CheckSolution.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index b38f838..7e1bf80 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -53,7 +53,7 @@ function compute_SOS(Q::AbstractArray, RG::GroupRing, l::Int) return GroupRingElem(result, RG) end -function distance_to_cone{S<:Interval}(elt::GroupRingElem, Q::AbstractArray{S,2}, wlen::Int) +function distance_to_cone{T<:Interval}(elt::GroupRingElem, Q::AbstractArray{T,2}, wlen::Int) SOS = compute_SOS(Q, parent(elt), length(elt.coeffs)) SOS_diff = elt - SOS From a00029007498175621c55446fd77a662264b89d5 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 3 Nov 2017 16:41:12 +0100 Subject: [PATCH 27/55] remove rationalization functions --- src/CheckSolution.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index 7e1bf80..ac99838 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -14,14 +14,6 @@ end (±)(X::GroupRingElem, tol::Real) = GroupRingElem(X.coeffs ± tol, parent(X)) -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) - EOI{T<:Number}(Δ::GroupRingElem{T}, λ::T) = Δ*Δ - λ*Δ function groupring_square(vect::AbstractVector, l, pm) From fba2464a8097280dbb9b7a51c496faa6af82a81c Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 3 Nov 2017 16:42:53 +0100 Subject: [PATCH 28/55] Specify concrete Interval types in augIdproj --- src/CheckSolution.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index ac99838..6e17050 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -75,7 +75,7 @@ end function augIdproj{T, I<:AbstractInterval}(S::Type{I}, Q::AbstractArray{T,2}) l = size(Q, 2) - R = zeros(Interval, (l,l)) + R = zeros(S, (l,l)) Threads.@threads for j in 1:l col = sum(view(Q, :,j))/l for i in 1:l @@ -87,7 +87,7 @@ end function augIdproj{T}(Q::AbstractArray{T,2}, logger) info(logger, "Projecting columns of Q to the augmentation ideal...") - @logtime logger Q = augIdproj(Interval, Q) + @logtime logger Q = augIdproj(Interval{T}, Q) info(logger, "Checking that sum of every column contains 0.0... ") check = all([0.0 in sum(view(Q, :, i)) for i in 1:size(Q, 2)]) From c0dc6937a51444b9a6bbec16edd7d918cb4a3aba Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 3 Nov 2017 16:54:41 +0100 Subject: [PATCH 29/55] cut on allocations in compute_SOS --- src/CheckSolution.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index 6e17050..d9078a6 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -18,25 +18,27 @@ EOI{T<:Number}(Δ::GroupRingElem{T}, λ::T) = Δ*Δ - λ*Δ function groupring_square(vect::AbstractVector, l, pm) zzz = zeros(eltype(vect), l) - zzz[1:length(vect)] .= vect - return GroupRings.mul!(similar(zzz), zzz, zzz, pm) + return GroupRings.mul!(zzz, vect, vect, pm) end function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int) - n = size(Q,2) # result = zeros(eltype(Q), l) - # for i in 1:n - # result .+= groupring_square(view(Q,:,i), l, pm) + # r = similar(result) + # for i in 1:size(Q,2) + # print(" $i") + # result += GroupRings.mul!(r, view(Q,:,i), view(Q,:,i), pm) # end @everywhere groupring_square = PropertyT.groupring_square - result = @parallel (+) for i in 1:n + result = @parallel (+) for i in 1:size(Q,2) print(" $i") groupring_square(view(Q,:,i), l, pm) end + println("") + return result end From c3ea80e6a9f786fcfc80a4bdd1b01cf83fc507d8 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 5 Nov 2017 16:38:51 +0100 Subject: [PATCH 30/55] replace remaining @time with '@logtime logger' --- src/OrbitDecomposition.jl | 16 ++++++++-------- src/PropertyT.jl | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index 2cd28ec..ccfc503 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -216,13 +216,13 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S # TODO: Fix that by multiple dispatch? Id = (isa(G, Nemo.Ring) ? one(G) : G()) - @time E4, sizes = Groups.generate_balls(S, Id, radius=2*radius); + @logtime logger E4, sizes = Groups.generate_balls(S, Id, radius=2*radius); info(logger, "Balls of sizes $sizes.") info(logger, "Reverse dict") - @time E_dict = GroupRings.reverse_dict(E4) + @logtime logger E_dict = GroupRings.reverse_dict(E4) info(logger, "Product matrix") - @time pm = GroupRings.create_pm(E4, E_dict, sizes[radius], twisted=true) + @logtime logger pm = GroupRings.create_pm(E4, E_dict, sizes[radius], twisted=true) RG = GroupRing(G, E4, E_dict, pm) Δ = PropertyT.splaplacian(RG, S) @assert GroupRings.augmentation(Δ) == 0 @@ -230,20 +230,20 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S save(joinpath(name, "pm.jld"), "pm", pm) info(logger, "Decomposing E into orbits of $(AutS)") - @time orbs = orbit_decomposition(AutS, E4, E_dict) + @logtime logger orbs = orbit_decomposition(AutS, E4, E_dict) @assert sum(length(o) for o in orbs) == length(E4) save(joinpath(name, "orbits.jld"), "orbits", orbs) info(logger, "Action matrices") - @time AutS_mreps = matrix_reps(AutS, E4[1:sizes[radius]], E_dict) + @logtime logger AutS_mreps = matrix_reps(AutS, E4[1:sizes[radius]], E_dict) info(logger, "Projections") - @time AutS_mps = rankOne_projections(AutS); + @logtime logger AutS_mps = rankOne_projections(AutS); - @time π_E_projections = [Cstar_repr(p, AutS_mreps) for p in AutS_mps] + @logtime logger π_E_projections = [Cstar_repr(p, AutS_mreps) for p in AutS_mps] info(logger, "Uπs...") - @time Uπs = orthSVD.(π_E_projections) + @logtime logger Uπs = orthSVD.(π_E_projections) multiplicities = size.(Uπs,2) info(logger, "multiplicities = $multiplicities") diff --git a/src/PropertyT.jl b/src/PropertyT.jl index cc117c2..2c20fe7 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -106,8 +106,8 @@ function ΔandSDPconstraints{T<:GroupElem}(name::String, S::Vector{T}, Id::T; ra end function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2) - B, sizes = Groups.generate_balls(S, Id, radius=2*radius) - info(logger, "Generated balls of sizes $sizes") + info(logger, "Generating balls of sizes $sizes") + @logtime logger B, sizes = Groups.generate_balls(S, Id, radius=2*radius) info(logger, "Creating product matrix...") @logtime logger pm = GroupRings.create_pm(B, GroupRings.reverse_dict(B), sizes[radius]; twisted=true) From f7f9ceef070b79637b04c161ec25123a57bce10f Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 5 Nov 2017 19:28:51 +0100 Subject: [PATCH 31/55] pass only the i-th column in @parallel mode this cuts on memory per worker drastically --- src/CheckSolution.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/CheckSolution.jl b/src/CheckSolution.jl index d9078a6..0a7e32e 100644 --- a/src/CheckSolution.jl +++ b/src/CheckSolution.jl @@ -33,8 +33,7 @@ function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int) @everywhere groupring_square = PropertyT.groupring_square result = @parallel (+) for i in 1:size(Q,2) - print(" $i") - groupring_square(view(Q,:,i), l, pm) + groupring_square(Q[:,i], l, pm) end println("") From 36e5979774bbf7f0bacf082c397b4b97301903da Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 5 Nov 2017 19:30:11 +0100 Subject: [PATCH 32/55] @logtime Q = sqrtm(P) --- src/PropertyT.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 2c20fe7..2c2f7a4 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -218,7 +218,7 @@ function check_property_T(name::String, S, Id, solver, upper_bound, tol, radius) isapprox(eigvals(P), abs(eigvals(P)), atol=tol) || warn("The solution matrix doesn't seem to be positive definite!") # @assert P == Symmetric(P) - Q = real(sqrtm(Symmetric(P))) + @logtime logger Q = real(sqrtm(Symmetric(P))) sgap = distance_to_positive_cone(Δ, λ, Q, 2*radius) if isa(sgap, Interval) From bc8c28170e7b7e6609ad62db73572daf63eba1cb Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 5 Nov 2017 20:54:33 +0100 Subject: [PATCH 33/55] missing @logtime in orbit-check_property_T --- src/Orbit-wise.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 3846310..69d7896 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -240,7 +240,7 @@ function check_property_T(sett::Settings) isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || warn("The solution matrix doesn't seem to be positive definite!") # @assert P == Symmetric(P) - Q = real(sqrtm(Symmetric(P))) + @logtime logger Q = real(sqrtm(Symmetric(P))) sgap = distance_to_positive_cone(Δ, λ, Q, 2*sett.radius) if isa(sgap, Interval) From 84813fedc4fc557d9acc0b7ba7063b16bdbf27f0 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 5 Nov 2017 20:55:53 +0100 Subject: [PATCH 34/55] indentation --- src/Orbit-wise.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 69d7896..1cfca2b 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -123,9 +123,9 @@ end A(data::OrbitData, π, t) = data.dims[π].*transform(data.Us[π], data.cnstr[t]) function constrLHS(m::JuMP.Model, data::OrbitData, t) - l = endof(data.Us) - lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l)) - return lhs + l = endof(data.Us) + lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l)) + return lhs end function constrLHS(m::JuMP.Model, cnstr, Us, Ust, dims, vars, eps=100*eps(1.0)) @@ -160,19 +160,19 @@ function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.lapl end function init_model(Uπs) - m = JuMP.Model(); - l = size(Uπs,1) - P = Vector{Array{JuMP.Variable,2}}(l) + m = JuMP.Model(); + l = size(Uπs,1) + P = Vector{Array{JuMP.Variable,2}}(l) - for k in 1:l - s = size(Uπs[k],2) - P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) - JuMP.@SDconstraint(m, P[k] >= 0.0) - end + for k in 1:l + s = size(Uπs[k],2) + P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) + JuMP.@SDconstraint(m, P[k] >= 0.0) + end - JuMP.@variable(m, λ >= 0.0) - JuMP.@objective(m, Max, λ) - return m, P + JuMP.@variable(m, λ >= 0.0) + JuMP.@objective(m, Max, λ) + return m, P end function create_SDP_problem(sett::Settings) @@ -238,7 +238,7 @@ function check_property_T(sett::Settings) Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || - warn("The solution matrix doesn't seem to be positive definite!") + warn("The solution matrix doesn't seem to be positive definite!") # @assert P == Symmetric(P) @logtime logger Q = real(sqrtm(Symmetric(P))) From 12c5c47ca8f30fa9a405b7a5efbdbe24889042cb Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 6 Nov 2017 11:31:15 +0100 Subject: [PATCH 35/55] add Nemo.isone(::GroupElem) --- src/Projections.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Projections.jl b/src/Projections.jl index 06c0b8a..5755f30 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -20,6 +20,8 @@ function (chi::PermCharacter)(g::Nemo.perm) return Int(Nemo.MN1inner(R, p, 1, Nemo._charvalsTable)) end +Nemo.isone(p::GroupElem) = p == parent(p)() + ## NOTE: this works only for Z/2!!!! function (chi::DirectProdCharacter)(g::DirectProductGroupElem) return reduce(*, 1, ((-1)^isone(g.elts[j]) for j in 1:chi.i)) From 31715701cca2f62385835c6fd3e76ac74f937064 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Tue, 7 Nov 2017 09:41:45 +0100 Subject: [PATCH 36/55] de-thread Projections to avoid UndefRefError: access to undefined reference in rankOne_projections(::WreathProduct) --- src/Projections.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Projections.jl b/src/Projections.jl index 5755f30..78c2c3f 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -128,7 +128,7 @@ function minimalprojections(G::PermutationGroup, T::Type=Rational{Int}) chars = [PropertyT.PermCharacter(p) for p in parts] min_projs = Vector{eltype(RGidems)}(l) - Threads.@threads for i in 1:l + for i in 1:l chi = PropertyT.PermCharacter(parts[i]) min_projs[i] = rankOne_projection(chi,RGidems)*central_projection(RG,chi) end From 0c9fb40e654f077b3e22961d2cb4ee2e310859f8 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:25:40 +0100 Subject: [PATCH 37/55] constraints_from_pm -> constraints returns a vector of constraints, i.e. vectors of Tuple{Int, Int} --- src/OrbitDecomposition.jl | 6 +++--- src/PropertyT.jl | 10 +++++----- src/SDPs.jl | 6 +++--- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index ccfc503..16a3ded 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -79,14 +79,14 @@ function orbit_spvector(vect::AbstractVector, orbits) return orb_vector end -function orbit_constraint(constraints::Vector{Vector{Vector{Int64}}}, n) +function orbit_constraint(constraints::Vector{Vector{Tuple{Int,Int}}}, n) result = spzeros(n,n) for cnstr in constraints for p in cnstr - result[p[2], p[1]] += 1.0 + result[p[2], p[1]] += 1.0/length(constraints) end end - return 1/length(constraints)*result + return result end ############################################################################### diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 2c2f7a4..32d4897 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -88,7 +88,7 @@ function ΔandSDPconstraints(prefix::String, G::Group) pm_fname, Δ_fname = pmΔfilenames(prefix) product_matrix = load(pm_fname, "pm") - sdp_constraints = constraints_from_pm(product_matrix) + sdp_constraints = constraints(product_matrix) RG = GroupRing(G, product_matrix) Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG) @@ -107,15 +107,15 @@ end function ΔandSDPconstraints{T<:GroupElem}(S::Vector{T}, Id::T; radius::Int=2) info(logger, "Generating balls of sizes $sizes") - @logtime logger B, sizes = Groups.generate_balls(S, Id, radius=2*radius) + @logtime logger E_R, sizes = Groups.generate_balls(S, Id, radius=2*radius) info(logger, "Creating product matrix...") - @logtime logger pm = GroupRings.create_pm(B, GroupRings.reverse_dict(B), sizes[radius]; twisted=true) + @logtime logger pm = GroupRings.create_pm(E_R, GroupRings.reverse_dict(E_R), sizes[radius]; twisted=true) info(logger, "Creating sdp_constratints...") - @logtime logger sdp_constraints = PropertyT.constraints_from_pm(pm) + @logtime logger sdp_constraints = PropertyT.constraints(pm) - RG = GroupRing(parent(Id), B, pm) + RG = GroupRing(parent(Id), E_R, pm) Δ = splaplacian(RG, S) return Δ, sdp_constraints diff --git a/src/SDPs.jl b/src/SDPs.jl index 972f3b5..a12d166 100644 --- a/src/SDPs.jl +++ b/src/SDPs.jl @@ -1,13 +1,13 @@ using JuMP import MathProgBase: AbstractMathProgSolver -function constraints_from_pm(pm, total_length=maximum(pm)) +function constraints(pm, total_length=maximum(pm)) n = size(pm,1) - constraints = [Array{Int,1}[] for x in 1:total_length] + constraints = [Vector{Tuple{Int,Int}}() for _ in 1:total_length] for j in 1:n for i in 1:n idx = pm[i,j] - push!(constraints[idx], [i,j]) + push!(constraints[idx], (i,j)) end end return constraints From e636c977388e860c10304e232e6b0b1d78b81533 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:28:04 +0100 Subject: [PATCH 38/55] unified marix_reps and perm_reps --- src/OrbitDecomposition.jl | 86 +++++++++++++++++++-------------------- 1 file changed, 42 insertions(+), 44 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index 16a3ded..51f14fa 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -91,77 +91,75 @@ end ############################################################################### # -# Matrix- and C*-representations +# Matrix-, Permutation- and C*-representations # ############################################################################### -function matrix_repr(g::GroupElem, E, E_dict) - rep_matrix = spzeros(Int, length(E), length(E)) +function matrix_repr(p::perm) + N = parent(p).n + return sparse(1:N, p.d, [1.0 for _ in 1:N]) +end + +function matrix_repr(g::GroupElem, E, E_rdict) + repmat = spzeros(Int, length(E), length(E)) for (i,elt) in enumerate(E) - j = E_dict[g(elt)] - rep_matrix[i,j] = 1 + j = E_rdict[g(elt)] + repmat[i,j] = 1 end - return rep_matrix + return repmat end -function matrix_reps{T<:GroupElem}(G::Group, S::Vector{T}, AutS::Group, radius::Int) - Id = (isa(G, Nemo.Ring) ? one(G) : G()) - E2, _ = Groups.generate_balls(S, Id, radius=radius) - Edict = GroupRings.reverse_dict(E2) - - elts = collect(elements(AutS)) - l = length(elts) - mreps = Vector{SparseMatrixCSC{Int, Int}}(l) - - Threads.@threads for i in 1:l - mreps[i] = PropertyT.matrix_repr(elts[i], E2, Edict) - end - - mreps_dict = Dict(elts[i]=>mreps[i] for i in 1:l) - - return mreps_dict -end - -function matrix_reps(G::Group, E2, E_dict) +function matrix_reps(G::Nemo.Group, E, E_rdict=GroupRings.reverse_dict(E)) elts = collect(elements(G)) l = length(elts) mreps = Vector{SparseMatrixCSC{Int, Int}}(l) Threads.@threads for i in 1:l - mreps[i] = matrix_repr(elts[i], E2, E_dict) + mreps[i] = matrix_repr(elts[i], E, E_rdict) end return Dict(elts[i]=>mreps[i] for i in 1:l) end -function perm_reps{T<:GroupElem}(G::Group, S::Vector{T}, AutS::Group, radius::Int) - Id = (isa(G, Nemo.Ring) ? one(G) : G()) - E_R, _ = Groups.generate_balls(S, Id, radius=radius) - Edict = GroupRings.reverse_dict(E_R) +function matrix_reps{T<:GroupElem}(autS::Nemo.Group, S::Vector{T}, radius::Int) + E, _ = Groups.generate_balls(S, radius=radius) + return matrix_reps(autS, E) +end - elts = collect(elements(AutS)) +function matrix_reps{T<:GroupElem}(preps::Dict{T,perm}) + kk = collect(keys(preps)) + mreps = Vector{SparseMatrixCSC{Float64, Int}}(length(kk)) + Threads.@threads for i in 1:length(kk) + mreps[i] = matrix_repr(preps[kk[i]]) + end + return Dict(kk[i] => mreps[i] for i in 1:length(kk)) +end + +function perm_repr(g::GroupElem, E::Vector, E_dict) + p = Vector{Int}(length(E)) + for (i,elt) in enumerate(E) + p[i] = E_dict[g(elt)] + end + return p +end + +function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E)) + elts = collect(elements(G)) l = length(elts) preps = Vector{Nemo.perm}(l) - G = Nemo.PermutationGroup(length(E_R)) + permG = Nemo.PermutationGroup(length(E)) Threads.@threads for i in 1:l - preps[i] = G(perm_repr(elts[i], E_R, Edict)) + preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict)) end - preps_dict = Dict(elts[i]=>preps[i] for i in 1:l) - - return preps_dict + return Dict(elts[i]=>preps[i] for i in 1:l) end -function perm_repr(g::GroupElem, E, E_dict) - l = length(E) - p = Vector{Int}(l) - for (i,elt) in enumerate(E) - j = E_dict[g(elt)] - p[i] = j - end - return p +function perm_reps(S::Vector, AutS::Group, radius::Int) + E, _ = Groups.generate_balls(S, radius=radius) + return perm_reps(AutS, E) end function reconstruct_sol{T<:GroupElem, S<:Nemo.perm}(preps::Dict{T, S}, From 2d02c98947f12f47cfcd28371323a83d763778b5 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:29:12 +0100 Subject: [PATCH 39/55] faster & simpler Cstar_repr --- src/OrbitDecomposition.jl | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index 51f14fa..64ef886 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -187,15 +187,7 @@ function reconstruct_sol{T<:GroupElem, S<:Nemo.perm}(preps::Dict{T, S}, end function Cstar_repr{T}(x::GroupRingElem{T}, mreps::Dict) - res = spzeros(size(mreps[first(keys(mreps))])...) - - for g in parent(x).basis - if x[g] != zero(T) - res .+= Float64(x[g]).*mreps[g] - end - end - - return res + return sum(x[g].*mreps[g] for g in parent(x).basis if x[g] != zero(T)) end function orthSVD(M::AbstractMatrix) From 72261afdebd4885114e333956787e48971b253ee Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:33:00 +0100 Subject: [PATCH 40/55] use preps_reprs instead of matrix_reprs and save them in preps.jld --- src/Orbit-wise.jl | 16 ++++++++++++++-- src/OrbitDecomposition.jl | 4 +++- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 1cfca2b..bbee2b9 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -103,7 +103,7 @@ function init_orbit_data(logger, sett::Settings; radius=2) ex(fname) = isfile(joinpath(prepath(sett), fname)) - files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld"]) + files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld", "preps.jld"]) if !all(files_exists) compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.AutS, radius=radius) @@ -204,7 +204,7 @@ function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) info(logger, "Reconstructing P...") - @logtime logger preps = perm_reps(sett.G, sett.S, sett.AutS, sett.radius) + preps = load_preps(joinpath(prepath(sett), "preps.jld"), sett.AutS) @logtime logger recP = reconstruct_sol(preps, data.Us, Ps, data.dims) @@ -213,6 +213,18 @@ function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) return λ, recP end +function load_preps(fname::String, G::Nemo.Group) + lded_preps = load(fname, "perms_d") + permG = PermutationGroup(length(first(lded_preps))) + @assert length(lded_preps) == order(G) + return Dict(k=>permG(v) for (k,v) in zip(elements(G), lded_preps)) +end + +function save_preps(fname::String, preps) + autS = parent(first(keys(preps))) + JLD.save(fname, "perms_d", [preps[elt].d for elt in elements(autS)]) +end + function check_property_T(sett::Settings) init_orbit_data(logger, sett, radius=sett.radius) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index 64ef886..b5ff66f 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -225,7 +225,9 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S save(joinpath(name, "orbits.jld"), "orbits", orbs) info(logger, "Action matrices") - @logtime logger AutS_mreps = matrix_reps(AutS, E4[1:sizes[radius]], E_dict) + @logtime logger reps = perm_reps(AutS, E_2R[1:sizes[radius]], E_rdict) + save_preps(joinpath(name, "preps.jld"), reps) + reps = matrix_reps(reps) info(logger, "Projections") @logtime logger AutS_mps = rankOne_projections(AutS); From 80deaebb9f20362b02ae9a69fe72bfafc8542270 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:35:26 +0100 Subject: [PATCH 41/55] remove matrix_reprs - we create them through perm_reprs --- src/OrbitDecomposition.jl | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index b5ff66f..a9696ed 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -100,32 +100,6 @@ function matrix_repr(p::perm) return sparse(1:N, p.d, [1.0 for _ in 1:N]) end -function matrix_repr(g::GroupElem, E, E_rdict) - repmat = spzeros(Int, length(E), length(E)) - for (i,elt) in enumerate(E) - j = E_rdict[g(elt)] - repmat[i,j] = 1 - end - return repmat -end - -function matrix_reps(G::Nemo.Group, E, E_rdict=GroupRings.reverse_dict(E)) - elts = collect(elements(G)) - l = length(elts) - mreps = Vector{SparseMatrixCSC{Int, Int}}(l) - - Threads.@threads for i in 1:l - mreps[i] = matrix_repr(elts[i], E, E_rdict) - end - - return Dict(elts[i]=>mreps[i] for i in 1:l) -end - -function matrix_reps{T<:GroupElem}(autS::Nemo.Group, S::Vector{T}, radius::Int) - E, _ = Groups.generate_balls(S, radius=radius) - return matrix_reps(autS, E) -end - function matrix_reps{T<:GroupElem}(preps::Dict{T,perm}) kk = collect(keys(preps)) mreps = Vector{SparseMatrixCSC{Float64, Int}}(length(kk)) From 95e19e6f9d48a11d0c6e921f0aebea735d25fa5e Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:36:34 +0100 Subject: [PATCH 42/55] replace stray constraints_from_pm() -> constraints() --- src/Orbit-wise.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index bbee2b9..19ee896 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -33,7 +33,7 @@ end function OrbitData(sett::Settings) splap = load(joinpath(prepath(sett), "delta.jld"), "Δ"); pm = load(joinpath(prepath(sett), "pm.jld"), "pm"); - cnstr = PropertyT.constraints_from_pm(pm); + cnstr = PropertyT.constraints(pm); splap² = similar(splap) splap² = GroupRings.mul!(splap², splap, splap, pm); From ff14b331995089bb5a838ab19a787a7fd5cc5ed8 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:37:39 +0100 Subject: [PATCH 43/55] AutS_mreps -> reps --- src/OrbitDecomposition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index a9696ed..e903483 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -206,7 +206,7 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S info(logger, "Projections") @logtime logger AutS_mps = rankOne_projections(AutS); - @logtime logger π_E_projections = [Cstar_repr(p, AutS_mreps) for p in AutS_mps] + @logtime logger π_E_projections = [Cstar_repr(p, reps) for p in AutS_mps] info(logger, "Uπs...") @logtime logger Uπs = orthSVD.(π_E_projections) From 3c51a463dc25ed68bdd8611f114ab680345fba1b Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:38:56 +0100 Subject: [PATCH 44/55] rename variables --- src/OrbitDecomposition.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index e903483..d310c9d 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -180,22 +180,22 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S # TODO: Fix that by multiple dispatch? Id = (isa(G, Nemo.Ring) ? one(G) : G()) - @logtime logger E4, sizes = Groups.generate_balls(S, Id, radius=2*radius); + @logtime logger E_2R, sizes = Groups.generate_balls(S, Id, radius=2*radius); info(logger, "Balls of sizes $sizes.") info(logger, "Reverse dict") - @logtime logger E_dict = GroupRings.reverse_dict(E4) + @logtime logger E_rdict = GroupRings.reverse_dict(E_2R) info(logger, "Product matrix") - @logtime logger pm = GroupRings.create_pm(E4, E_dict, sizes[radius], twisted=true) - RG = GroupRing(G, E4, E_dict, pm) + @logtime logger pm = GroupRings.create_pm(E_2R, E_rdict, sizes[radius], twisted=true) + RG = GroupRing(G, E_2R, E_rdict, pm) Δ = PropertyT.splaplacian(RG, S) @assert GroupRings.augmentation(Δ) == 0 save(joinpath(name, "delta.jld"), "Δ", Δ.coeffs) save(joinpath(name, "pm.jld"), "pm", pm) info(logger, "Decomposing E into orbits of $(AutS)") - @logtime logger orbs = orbit_decomposition(AutS, E4, E_dict) - @assert sum(length(o) for o in orbs) == length(E4) + @logtime logger orbs = orbit_decomposition(AutS, E_2R, E_rdict) + @assert sum(length(o) for o in orbs) == length(E_2R) save(joinpath(name, "orbits.jld"), "orbits", orbs) info(logger, "Action matrices") From 4ba2f54bb41b3a1ec6a81da5516e707732205853 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:50:58 +0100 Subject: [PATCH 45/55] improve slightly orthSVD --- src/OrbitDecomposition.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index d310c9d..da1f28d 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -164,12 +164,11 @@ function Cstar_repr{T}(x::GroupRingElem{T}, mreps::Dict) return sum(x[g].*mreps[g] for g in parent(x).basis if x[g] != zero(T)) end -function orthSVD(M::AbstractMatrix) - M = full(M) - fact = svdfact(M) - singv = fact[:S] - M_rank = sum(singv .> maximum(size(M))*eps(eltype(singv))) - return fact[:U][:,1:M_rank] +function orthSVD{T}(M::AbstractMatrix{T}) + M = full(M) + fact = svdfact(M) + M_rank = sum(fact[:S] .> maximum(size(M))*eps(T)) + return fact[:U][:,1:M_rank] end function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S::Vector{T}, AutS; radius=2) From c1dc3850dde534860aac81ab25d5bd64b14f6490 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:56:28 +0100 Subject: [PATCH 46/55] rename AutS -> autS --- src/Orbit-wise.jl | 6 +++--- src/OrbitDecomposition.jl | 18 +++++++++--------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index 19ee896..b4d3e9a 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -8,7 +8,7 @@ immutable Settings N::Int G::Group S::Vector - AutS::Group + autS::Group radius::Int solver::SCSSolver upper_bound::Float64 @@ -106,7 +106,7 @@ function init_orbit_data(logger, sett::Settings; radius=2) files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld", "preps.jld"]) if !all(files_exists) - compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.AutS, radius=radius) + compute_orbit_data(logger, prepath(sett), sett.G, sett.S, sett.autS, radius=radius) end return 0 @@ -204,7 +204,7 @@ function λandP(m::JuMP.Model, data::OrbitData, sett::Settings) info(logger, "Reconstructing P...") - preps = load_preps(joinpath(prepath(sett), "preps.jld"), sett.AutS) + preps = load_preps(joinpath(prepath(sett), "preps.jld"), sett.autS) @logtime logger recP = reconstruct_sol(preps, data.Us, Ps, data.dims) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index da1f28d..7e444e1 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -131,9 +131,9 @@ function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E)) return Dict(elts[i]=>preps[i] for i in 1:l) end -function perm_reps(S::Vector, AutS::Group, radius::Int) +function perm_reps(S::Vector, autS::Group, radius::Int) E, _ = Groups.generate_balls(S, radius=radius) - return perm_reps(AutS, E) + return perm_reps(autS, E) end function reconstruct_sol{T<:GroupElem, S<:Nemo.perm}(preps::Dict{T, S}, @@ -171,7 +171,7 @@ function orthSVD{T}(M::AbstractMatrix{T}) return fact[:U][:,1:M_rank] end -function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S::Vector{T}, AutS; radius=2) +function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S::Vector{T}, autS::Nemo.Group; radius=2) isdir(name) || mkdir(name) info(logger, "Generating ball of radius $(2*radius)") @@ -192,27 +192,27 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S save(joinpath(name, "delta.jld"), "Δ", Δ.coeffs) save(joinpath(name, "pm.jld"), "pm", pm) - info(logger, "Decomposing E into orbits of $(AutS)") - @logtime logger orbs = orbit_decomposition(AutS, E_2R, E_rdict) + info(logger, "Decomposing E into orbits of $(autS)") + @logtime logger orbs = orbit_decomposition(autS, E_2R, E_rdict) @assert sum(length(o) for o in orbs) == length(E_2R) save(joinpath(name, "orbits.jld"), "orbits", orbs) info(logger, "Action matrices") - @logtime logger reps = perm_reps(AutS, E_2R[1:sizes[radius]], E_rdict) + @logtime logger reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict) save_preps(joinpath(name, "preps.jld"), reps) reps = matrix_reps(reps) info(logger, "Projections") - @logtime logger AutS_mps = rankOne_projections(AutS); + @logtime logger autS_mps = rankOne_projections(autS); - @logtime logger π_E_projections = [Cstar_repr(p, reps) for p in AutS_mps] + @logtime logger π_E_projections = [Cstar_repr(p, reps) for p in autS_mps] info(logger, "Uπs...") @logtime logger Uπs = orthSVD.(π_E_projections) multiplicities = size.(Uπs,2) info(logger, "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]; info(logger, "dimensions = $dimensions") @assert dot(multiplicities, dimensions) == sizes[radius] From 2aa82db1ac4a2ce5b2a16e20631d3155c16b1f4a Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 09:58:37 +0100 Subject: [PATCH 47/55] specify exact types for OrbitData --- src/Orbit-wise.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index b4d3e9a..bf17026 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -20,13 +20,13 @@ suffix(s::Settings) = "$(s.upper_bound)" prepath(s::Settings) = prefix(s) fullpath(s::Settings) = joinpath(prefix(s), suffix(s)) -immutable OrbitData +immutable OrbitData{T<:AbstractArray{Float64, 2}, LapType <:AbstractVector{Float64}} name::String - Us::Vector + Us::Vector{T} Ps::Vector{Array{JuMP.Variable,2}} - cnstr::Vector - laplacian::Vector - laplacianSq::Vector + cnstr::Vector{SparseMatrixCSC{Float64, Int}} + laplacian::LapType + laplacianSq::LapType dims::Vector{Int} end From fe1e9ee37c18e6efcb6354409732be9d871d68bf Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 10:00:34 +0100 Subject: [PATCH 48/55] let sparsify actually return sparse matrix --- src/Orbit-wise.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index bf17026..a387703 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -94,10 +94,10 @@ function sparsify!{T}(M::AbstractArray{T}, eps=eps(T); check=false, verbose=fals info(logger, "Sparsified density:", rpad(densM, 20), " → ", rpad(dens(M),20)) end - return M + return sparse(M) end -sparsify{T}(U::AbstractArray{T}, tol=eps(T)) = sparsify!(deepcopy(U), tol) +sparsify{T}(U::AbstractArray{T}, tol=eps(T); check=true, verbose=false) = sparsify!(deepcopy(U), tol, check=check, verbose=verbose) function init_orbit_data(logger, sett::Settings; radius=2) From bcbf8a03c8db8f4d1e3428db3b9a18a9cf2497bd Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 10:01:45 +0100 Subject: [PATCH 49/55] =?UTF-8?q?store=20the=20dense=20U=CF=80s,=20sparsif?= =?UTF-8?q?y=20at=20read-time?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Orbit-wise.jl | 3 ++- src/OrbitDecomposition.jl | 1 - 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index a387703..c5175af 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -38,7 +38,8 @@ function OrbitData(sett::Settings) splap² = GroupRings.mul!(splap², splap, splap, pm); # Uπs = load(joinpath(name, "U_pis.jld"), "Uπs"); - Uπs = load(joinpath(prepath(sett), "U_pis.jld"), "spUπs"); + Uπs = load(joinpath(prepath(sett), "U_pis.jld"), "Uπs") + Uπs = sparsify!.(Uπs, sett.tol, check=true, verbose=true) #dimensions of the corresponding πs: dims = load(joinpath(prepath(sett), "U_pis.jld"), "dims") diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index 7e444e1..17c1c0e 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -218,7 +218,6 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, - "spUπs", sparsify!.(deepcopy(Uπs), check=true, verbose=true), "dims", dimensions) return 0 end From ea6a6722be1f5cd587e5804d348ca87720e4002a Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 10:02:38 +0100 Subject: [PATCH 50/55] =?UTF-8?q?init=20model=20on=20number=20and=20size?= =?UTF-8?q?=20of=20SD=20blocks,=20not=20on=20U=CF=80s?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Orbit-wise.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl index c5175af..d77e556 100644 --- a/src/Orbit-wise.jl +++ b/src/Orbit-wise.jl @@ -43,7 +43,7 @@ function OrbitData(sett::Settings) #dimensions of the corresponding πs: dims = load(joinpath(prepath(sett), "U_pis.jld"), "dims") - m, P = init_model(Uπs); + m, P = init_model(size(Uπs,1), [size(U,2) for U in Uπs]); orbits = load(joinpath(prepath(sett), "orbits.jld"), "orbits"); n = size(Uπs[1],1) @@ -160,13 +160,11 @@ function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.lapl println("") end -function init_model(Uπs) +function init_model(n, sizes) m = JuMP.Model(); - l = size(Uπs,1) - P = Vector{Array{JuMP.Variable,2}}(l) + P = Vector{Array{JuMP.Variable,2}}(n) - for k in 1:l - s = size(Uπs[k],2) + for (k,s) in enumerate(sizes) P[k] = JuMP.@variable(m, [i=1:s, j=1:s]) JuMP.@SDconstraint(m, P[k] >= 0.0) end From 6f0087b8c70f2a0f160de60fb21faf9fb8e43bd3 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 11:10:30 +0100 Subject: [PATCH 51/55] indentation --- src/Projections.jl | 92 +++++++++++++++++++++++----------------------- 1 file changed, 45 insertions(+), 47 deletions(-) diff --git a/src/Projections.jl b/src/Projections.jl index 78c2c3f..9382e02 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -49,47 +49,45 @@ end # ############################################################################### -function central_projection(RG::GroupRing, chi::AbstractCharacter, - T::Type=Rational{Int}) - result = RG(T) - result.coeffs = full(result.coeffs) - dim = chi(RG.group()) - ord = Int(order(RG.group)) +function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int}) + result = RG(T) + result.coeffs = full(result.coeffs) + dim = chi(RG.group()) + ord = Int(order(RG.group)) - for g in RG.basis - result[g] = convert(T, (dim//ord)*chi(g)) - end + for g in RG.basis + result[g] = convert(T, (dim//ord)*chi(g)) + end - return result + return result end function idempotents(RG::GroupRing{PermGroup}, T::Type=Rational{Int}) - if RG.group.n == 1 - return GroupRingElem{T}[one(RG,T)] - elseif RG.group.n == 2 - Id = one(RG,T) - transp = convert(T, RG(RG.group([2,1]))) - return GroupRingElem{T}[1//2*(Id + transp), 1//2*(Id - transp)] + if RG.group.n == 1 + return GroupRingElem{T}[one(RG,T)] + elseif RG.group.n == 2 + Id = one(RG,T) + transp = convert(T, RG(RG.group([2,1]))) + return GroupRingElem{T}[1//2*(Id + transp), 1//2*(Id - transp)] + end + projs = Vector{Vector{perm}}() + for l in 2:RG.group.n + u = RG.group([circshift([i for i in 1:l], -1); [i for i in l+1:RG.group.n]]) + i = 0 + while (l-1)*i <= RG.group.n + v = RG.group(circshift(collect(1:RG.group.n), i)) + k = inv(v)*u*v + push!(projs, generateGroup([k], RG.group.n)) + i += 1 + end + end - end - projs = Vector{Vector{perm}}() - for l in 2:RG.group.n - u = RG.group([circshift([i for i in 1:l], -1); [i for i in l+1:RG.group.n]]) - i = 0 - while (l-1)*i <= RG.group.n - v = RG.group(circshift(collect(1:RG.group.n), i)) - k = inv(v)*u*v - push!(projs, generateGroup([k], RG.group.n)) - i += 1 - end - end + idems = Vector{GroupRingElem{T}}() + for p in projs + append!(idems, [RG(p, T), RG(p, T, alt=true)]) + end - idems = Vector{GroupRingElem{T}}() - for p in projs - append!(idems, [RG(p, T), RG(p, T, alt=true)]) - end - - return unique(idems) + return unique(idems) end function rankOne_projection{S}(chi::PropertyT.PermCharacter, idems::Vector{GroupRingElem{S}}) @@ -163,7 +161,7 @@ function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int}) last_emb = g->BN(Nemo.emb!(BN.P(), g, range[i+1:end])) Sk_first = [RBN(p, first_emb) for p in SNprojs_nc[i]] - Sk_last = [RBN(p, last_emb ) for p in SNprojs_nc[N-i]] + Sk_last = [RBN(p, last_emb) for p in SNprojs_nc[N-i]] append!(all_projs, [Qs[i]*p1*p2 for (p1,p2) in Base.product(Sk_first,Sk_last)]) @@ -187,18 +185,18 @@ doc""" > forming 'products' by adding `op` (which is `*` by default). """ function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) - result = Vector{T}() - seen = Set{T}() - for x in X - for y in Y - z = op(x,y) - if !in(z, seen) - push!(seen, z) - push!(result, z) - end - end - end - return result + result = Vector{T}() + seen = Set{T}() + for x in X + for y in Y + z = op(x,y) + if !in(z, seen) + push!(seen, z) + push!(result, z) + end + end + end + return result end doc""" From d7652c4e0270ef86efe8dd569aac4cd19bad5463 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 15:38:51 +0100 Subject: [PATCH 52/55] make AbstractCharacter just an abstract type --- src/Projections.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Projections.jl b/src/Projections.jl index 9382e02..93b4ee0 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -4,7 +4,7 @@ # ############################################################################### -abstract AbstractCharacter <: Function +abstract AbstractCharacter immutable PermCharacter <: AbstractCharacter p::Partition From e746aaeeff382d2112c742950d066fb964d8f53b Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 15:39:51 +0100 Subject: [PATCH 53/55] rankOne_projections MUST be computed over exact field (e.g. Rational{Int}) --- src/Projections.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Projections.jl b/src/Projections.jl index 93b4ee0..ab1b852 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -142,7 +142,7 @@ function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int}) N = BN.P.n # projections as elements of the group rings RSₙ - SNprojs_nc = [rankOne_projections(PermutationGroup(i), T) for i in 1:N] + SNprojs_nc = [rankOne_projections(PermutationGroup(i)) for i in 1:N] # embedding into group ring of BN RBN = GroupRing(BN) From 38ae1e0656b8a65d10b1449cbc2c7f7ef3b0d27d Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 15:40:22 +0100 Subject: [PATCH 54/55] remove minimalprojections in favour of rankOne_projections --- src/Projections.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/Projections.jl b/src/Projections.jl index ab1b852..1a07cab 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -110,7 +110,7 @@ function rankOne_projection{S}(chi::PropertyT.PermCharacter, idems::Vector{Group throw("Couldn't find rank-one projection for $chi") end -function minimalprojections(G::PermutationGroup, T::Type=Rational{Int}) +function rankOne_projections(G::PermutationGroup, T::Type=Rational{Int}) if G.n == 1 return [one(GroupRing(G), T)] elseif G.n < 8 @@ -134,9 +134,6 @@ function minimalprojections(G::PermutationGroup, T::Type=Rational{Int}) return min_projs end -function rankOne_projections(G::PermutationGroup, T::Type=Rational{Int}) - return minimalprojections(G, T) -end function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int}) From a35122f8bb5c0a66eb012896cb1e78e35982859b Mon Sep 17 00:00:00 2001 From: kalmarek Date: Wed, 8 Nov 2017 15:40:38 +0100 Subject: [PATCH 55/55] cosmetic --- src/OrbitDecomposition.jl | 1 + src/Projections.jl | 9 ++++----- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl index 17c1c0e..ce582cf 100644 --- a/src/OrbitDecomposition.jl +++ b/src/OrbitDecomposition.jl @@ -195,6 +195,7 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S info(logger, "Decomposing E into orbits of $(autS)") @logtime logger orbs = orbit_decomposition(autS, E_2R, E_rdict) @assert sum(length(o) for o in orbs) == length(E_2R) + info(logger, "E consists of $(length(orbs)) orbits!") save(joinpath(name, "orbits.jld"), "orbits", orbs) info(logger, "Action matrices") diff --git a/src/Projections.jl b/src/Projections.jl index 1a07cab..8990ae6 100644 --- a/src/Projections.jl +++ b/src/Projections.jl @@ -90,19 +90,19 @@ function idempotents(RG::GroupRing{PermGroup}, T::Type=Rational{Int}) return unique(idems) end -function rankOne_projection{S}(chi::PropertyT.PermCharacter, idems::Vector{GroupRingElem{S}}) +function rankOne_projection{T}(chi::PropertyT.PermCharacter, idems::Vector{GroupRingElem{T}}) RG = parent(first(idems)) - ids = [[one(RG, S)]; idems] + ids = [[one(RG, T)]; idems] for (i,j,k) in Base.product(ids, ids, ids) - if chi(i) == zero(S) || chi(j) == zero(S) || chi(k) == zero(S) + if chi(i) == zero(T) || chi(j) == zero(T) || chi(k) == zero(T) continue end elt = i*j*k elt^2 == elt || continue - if chi(elt) == one(S) + if chi(elt) == one(T) return elt # return (i,j,k) end @@ -134,7 +134,6 @@ function rankOne_projections(G::PermutationGroup, T::Type=Rational{Int}) return min_projs end - function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int}) N = BN.P.n