Merge branch 'master' into enh/julia-v0.6

src/CheckSolution.jl

@@ -3,8 +3,7 @@ import Base: rationalize
 using IntervalArithmetic
 
 IntervalArithmetic.setrounding(Interval, :tight)
-IntervalArithmetic.setformat(sigfigs=10)
+IntervalArithmetic.setformat(sigfigs=12)
-IntervalArithmetic.setprecision(Interval, 53) # slightly faster than 256
 
 import IntervalArithmetic.±
 
@@ -15,131 +14,97 @@ 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)
     zzz = zeros(eltype(vect), l)
-    zzz[1:length(vect)] .= vect
+    return GroupRings.mul!(zzz, vect, vect, pm)
-    return GroupRings.mul!(similar(zzz), zzz, zzz, pm)
 end
 
-function compute_SOS(sqrt_matrix, elt::GroupRingElem)
+function compute_SOS(Q::AbstractArray, pm::Array{Int,2}, l::Int)
-   n = size(sqrt_matrix,2)
-   l = length(elt.coeffs)
-   pm = parent(elt).pm
 
-   result = zeros(eltype(sqrt_matrix), l)
+   # result = zeros(eltype(Q), l)
-   for i in 1:n
+   # r = similar(result)
-      result .+= groupring_square(view(sqrt_matrix,:,i), l, pm)
+   # for i in 1:size(Q,2)
-   end
+   #    print(" $i")
-
+   #    result += GroupRings.mul!(r, view(Q,:,i), view(Q,:,i), pm)
-   # @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))
+   @everywhere groupring_square = PropertyT.groupring_square
+
+   result = @parallel (+) for i in 1:size(Q,2)
+      groupring_square(Q[:,i], l, pm)
+   end
+
+   println("")
+
+   return result
 end
 
-function correct_to_augmentation_ideal{T<:Rational}(sqrt_matrix::Array{T,2})
+function compute_SOS(Q::AbstractArray, RG::GroupRing, l::Int)
-   l = size(sqrt_matrix, 2)
+   result = compute_SOS(Q, RG.pm, l)
-   sqrt_corrected = Array{Interval{Float64}}(l,l)
+   return GroupRingElem(result, RG)
+end
+
+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
+
+   ɛ_dist = GroupRings.augmentation(SOS_diff)
+   info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
+
+   eoi_SOS_L1_dist = norm(SOS_diff,1)
+   info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)")
+
+   dist = 2^(wlen-1)*eoi_SOS_L1_dist
+   return dist
+end
+
+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
+
+   ɛ_dist = GroupRings.augmentation(SOS_diff)
+   info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
+
+   eoi_SOS_L1_dist = norm(SOS_diff,1)
+   info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L1_dist))")
+
+   dist = 2^(wlen-1)*eoi_SOS_L1_dist
+   return dist
+end
+
+function augIdproj{T, I<:AbstractInterval}(S::Type{I}, Q::AbstractArray{T,2})
+   l = size(Q, 2)
+   R = zeros(S, (l,l))
    Threads.@threads for j in 1:l
-      col = sum(view(sqrt_matrix, :,j))//l
+      col = sum(view(Q, :,j))/l
       for i in 1:l
-         sqrt_corrected[i,j] = (Float64(sqrt_matrix[i,j]) - Float64(col)) ± eps(0.0)
+         R[i,j] = Q[i,j] - col ± eps(0.0)
       end
    end
-   return sqrt_corrected
+   return R
 end
 
-function distance_to_cone{T<:Rational}(λ::T, sqrt_matrix::Array{T,2}, Δ::GroupRingElem{T}, wlen)
+function augIdproj{T}(Q::AbstractArray{T,2}, logger)
-    SOS = compute_SOS(sqrt_matrix, Δ)
+   info(logger, "Projecting columns of Q to the augmentation ideal...")
-
+   @logtime logger Q = augIdproj(Interval{T}, Q)
-    SOS_diff = EOI(Δ, λ) - SOS
-    eoi_SOS_L1_dist = norm(SOS_diff,1)
-
-    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)))")
-
-    distance_to_cone = λ - 2^(wlen-1)*eoi_SOS_L1_dist
-    return distance_to_cone
-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)
-
-    info(logger, "λ = $λ (≈≥$(@sprintf("%.10f",float(λ))))")
-    ɛ_dist = GroupRings.augmentation(SOS_diff)
-
-    info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
-    info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L1_dist)")
-
-    distance_to_cone = λ - 2^(wlen-1)*eoi_SOS_L1_dist
-    return distance_to_cone
-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)
-   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_int = 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;
+function distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen::Int)
-    tol=1e-14, rational=false)
-
     info(logger, "------------------------------------------------------------")
-    info(logger, "")
+    info(logger, "λ = $λ")
     info(logger, "Checking in floating-point arithmetic...")
-    t = @timed fp_distance = distance_to_cone(λ, Q, Δ, wlen)
+    Δ²_λΔ = EOI(Δ, λ)
-    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, "------------------------------------------------------------")
 
@@ -148,26 +113,16 @@ function check_distance_to_positive_cone(Δ::GroupRingElem, λ, Q, wlen;
     end
 
     info(logger, "")
-    Q_ℚω_int = rationalize_and_project(Q, tol, logger)
+    Q = augIdproj(Q, logger)
-    λ_ℚ = ℚ(λ, tol)
-    Δ_ℚ = ℚ(Δ, tol)
 
     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)
+    @logtime logger 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
-        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
 end

src/Orbit-wise.jl

@@ -8,44 +8,50 @@ immutable Settings
    N::Int
    G::Group
    S::Vector
-   AutS::Group
+   autS::Group
    radius::Int
    solver::SCSSolver
    upper_bound::Float64
    tol::Float64
 end
 
-immutable OrbitData
+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{T<:AbstractArray{Float64, 2}, LapType <:AbstractVector{Float64}}
    name::String
-   Us::Vector
+   Us::Vector{T}
    Ps::Vector{Array{JuMP.Variable,2}}
-   cnstr::Vector
+   cnstr::Vector{SparseMatrixCSC{Float64, Int}}
-   laplacian::Vector
+   laplacian::LapType
-   laplacianSq::Vector
+   laplacianSq::LapType
    dims::Vector{Int}
 end
 
-function OrbitData(name::String)
+function OrbitData(sett::Settings)
-   splap = load(joinpath(name, "delta.jld"), "Δ");
+   splap = load(joinpath(prepath(sett), "delta.jld"), "Δ");
-   pm = load(joinpath(name, "pm.jld"), "pm");
+   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);
 
    # 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"), "Uπs")
+   Uπs = sparsify!.(Uπs, sett.tol, check=true, verbose=true)
    #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);
+   m, P = init_model(size(Uπs,1), [size(U,2) for U in 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);
 
@@ -89,19 +95,19 @@ 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)
 
-   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"])
+   files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld", "preps.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
@@ -118,9 +124,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)
+   l = endof(data.Us)
-    lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l))
+   lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l))
-    return lhs
+   return lhs
 end
 
 function constrLHS(m::JuMP.Model, cnstr, Us, Ust, dims, vars, eps=100*eps(1.0))
@@ -154,36 +160,32 @@ 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();
+   m = JuMP.Model();
-    l = size(Uπs,1)
+   P = Vector{Array{JuMP.Variable,2}}(n)
-    P = Vector{Array{JuMP.Variable,2}}(l)
 
-    for k in 1:l
+   for (k,s) in enumerate(sizes)
-        s = size(Uπs[k],2)
+      P[k] = JuMP.@variable(m, [i=1:s, j=1:s])
-        P[k] = JuMP.@variable(m, [i=1:s, j=1:s])
+      JuMP.@SDconstraint(m, P[k] >= 0.0)
-        JuMP.@SDconstraint(m, P[k] >= 0.0)
+   end
-    end
 
-    JuMP.@variable(m, λ >= 0.0)
+   JuMP.@variable(m, λ >= 0.0)
-    JuMP.@objective(m, Max, λ)
+   JuMP.@objective(m, Max, λ)
-    return m, P
+   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);
+   @logtime logger 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)
    info(logger, "Adding $t constraints ... ")
-   t = @timed addconstraints!(SDP_problem, orb_data)
+   @logtime logger addconstraints!(SDP_problem, orb_data)
-   info(logger, PropertyT.timed_msg(t))
 
    return SDP_problem, orb_data
 end
@@ -201,28 +203,36 @@ 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)
+   preps = load_preps(joinpath(prepath(sett), "preps.jld"), sett.autS)
-   info(logger, PropertyT.timed_msg(t))
 
-   t = @timed recP = reconstruct_sol(preps, data.Us, Ps, data.dims)
+   @logtime logger 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
 
+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)
 
-   fnames = PropertyT.λSDPfilenames(sett.name)
+   if all(isfile.(λSDPfilenames(fullpath(sett))))
-
+      λ, P = PropertyT.λandP(fullpath(sett))
-   if all(isfile.(fnames))
-      λ, 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)
@@ -234,17 +244,16 @@ function check_property_T(sett::Settings)
    info(logger, "minimum(P) = $(minimum(P))")
 
    if λ > 0
-      pm_fname = joinpath(sett.name, "pm.jld")
+      pm_fname, Δ_fname = pmΔfilenames(prepath(sett))
       RG = GroupRing(sett.G, load(pm_fname, "pm"))
-      Δ_fname = joinpath(sett.name, "delta.jld")
       Δ = 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)
-      Q = real(sqrtm(Symmetric(P)))
+      @logtime logger 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

src/OrbitDecomposition.jl

@@ -10,7 +10,7 @@ mutable struct FFEltsIter{T<:Generic.FinField}
     all::Int
     field::T
 
-    function FFEltsIter{T}(F::T) where {T} 
+    function FFEltsIter{T}(F::T) where {T}
         return new(Int(characteristic(F)^degree(F)), F)
     end
 end
@@ -79,89 +79,61 @@ 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
 
 ###############################################################################
 #
-#  Matrix- and C*-representations
+#  Matrix-, Permutation- and C*-representations
 #
 ###############################################################################
 
-function matrix_repr(g::GroupElem, E, E_dict)
+function matrix_repr(p::perm)
-   rep_matrix = spzeros(Int, length(E), length(E))
+    N = parent(p).n
+    return sparse(1:N, p.d, [1.0 for _ in 1:N])
+end
+
+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)
-      j = E_dict[g(elt)]
+      p[i] = E_dict[g(elt)]
-      rep_matrix[i,j] = 1
    end
-   return rep_matrix
+   return p
 end
 
-function matrix_reps{T<:GroupElem}(G::Group, S::Vector{T}, AutS::Group, radius::Int)
+function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
-   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)
    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)
-   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)
-
-   elts = collect(elements(AutS))
-   l = length(elts)
    preps = Vector{Generic.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 Dict(elts[i]=>preps[i] for i in 1:l)
-
-   return preps_dict
 end
 
-function perm_repr(g::GroupElem, E, E_dict)
+function perm_reps(S::Vector, autS::Group, radius::Int)
-   l = length(E)
+   E, _ = Groups.generate_balls(S, radius=radius)
-   p = Vector{Int}(l)
+   return perm_reps(autS, E)
-   for (i,elt) in enumerate(E)
-      j = E_dict[g(elt)]
-      p[i] = j
-   end
-   return p
 end
 
 function reconstruct_sol{T<:GroupElem}(preps::Dict{T, Generic.perm},
@@ -189,22 +161,17 @@ function reconstruct_sol{T<:GroupElem}(preps::Dict{T, Generic.perm},
 end
 
 function Cstar_repr(x::GroupRingElem{T}, mreps::Dict) where {T}
-   nzindx = [i for i in eachindex(x.coeffs) if x[i] != zero(T)]
+   return sum(x[i].*mreps[parent(x).basis[i]] for i in findn(x.coeffs))
-   RG = parent(x)
-   res = sum(Float64(x[i]).*mreps[RG.basis[i]] for i in nzindx)
-
-   return res
 end
 
-function orthSVD(M::AbstractMatrix)
+function orthSVD{T}(M::AbstractMatrix{T})
-    M = full(M)
+   M = full(M)
-    fact = svdfact(M)
+   fact = svdfact(M)
-    singv = fact[:S]
+   M_rank = sum(fact[:S] .> maximum(size(M))*eps(T))
-    M_rank = sum(singv .> maximum(size(M))*eps(eltype(singv)))
+   return fact[:U][:,1:M_rank]
-    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)")
@@ -212,44 +179,46 @@ 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 E_2R, 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_rdict = GroupRings.reverse_dict(E_2R)
 
    info(logger, "Product matrix")
-   @time pm = GroupRings.create_pm(E4, E_dict, sizes[radius], twisted=true)
+   @logtime logger pm = GroupRings.create_pm(E_2R, E_rdict, sizes[radius], twisted=true)
-   RG = GroupRing(G, E4, E_dict, pm)
+   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)")
+   info(logger, "Decomposing E into orbits of $(autS)")
-   @time orbs = orbit_decomposition(AutS, E4, E_dict)
+   @logtime logger orbs = orbit_decomposition(autS, E_2R, E_rdict)
-   @assert sum(length(o) for o in orbs) == length(E4)
+   @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")
-   @time 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")
-   @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, reps) 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")
-   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]
 
    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

src/Projections.jl

@@ -4,7 +4,7 @@
 #
 ###############################################################################
 
-abstract type AbstractCharacter <: Function end
+abstract type AbstractCharacter end
 
 struct PermCharacter <: AbstractCharacter
    p::Generic.Partition
@@ -20,6 +20,8 @@ function (chi::PermCharacter)(g::Generic.perm)
    return Int(Nemo.Generic.MN1inner(R, p, 1, Nemo.Generic._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))
@@ -47,18 +49,17 @@ end
 #
 ###############################################################################
 
-function central_projection(RG::GroupRing, chi::AbstractCharacter,
+function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int})
-    T::Type=Rational{Int})
+   result = RG(T)
-    result = RG(T)
+   result.coeffs = full(result.coeffs)
-    result.coeffs = full(result.coeffs)
+   dim = chi(RG.group())
-    dim = chi(RG.group())
+   ord = Int(order(RG.group))
-    ord = Int(order(RG.group))
 
-    for g in RG.basis
+   for g in RG.basis
-        result[g] = convert(T, (dim//ord)*chi(g))
+      result[g] = convert(T, (dim//ord)*chi(g))
-    end
+   end
 
-    return result
+   return result
 end
 
 function idempotents(RG::GroupRing{Generic.PermGroup}, T::Type=Rational{Int})
@@ -90,12 +91,10 @@ function idempotents(RG::GroupRing{Generic.PermGroup}, T::Type=Rational{Int})
     return unique(idems)
 end
 
-function rankOne_projection(chi::PropertyT.PermCharacter, idems::Vector{S}) where {S<:GroupRingElem}
+function rankOne_projection(chi::PropertyT.PermCharacter, idems::Vector{T}) where {T<:GroupRingElem}
 
    RG = parent(first(idems))
 
-   T = eltype(first(idems))
-
    ids = [[one(RG, T)]; idems]
 
    for (i,j,k) in Base.product(ids, ids, ids)
@@ -112,7 +111,7 @@ function rankOne_projection(chi::PropertyT.PermCharacter, idems::Vector{S}) wher
    throw("Couldn't find rank-one projection for $chi")
 end
 
-function minimalprojections(G::Generic.PermGroup, T::Type=Rational{Int})
+function rankOne_projections(G::Generic.PermGroup, T::Type=Rational{Int})
    if G.n == 1
       return [one(GroupRing(G), T)]
    elseif G.n < 8
@@ -128,7 +127,7 @@ function minimalprojections(G::Generic.PermGroup, 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
@@ -136,15 +135,11 @@ function minimalprojections(G::Generic.PermGroup, T::Type=Rational{Int})
    return min_projs
 end
 
-function rankOne_projections(G::Generic.PermGroup, T::Type=Rational{Int})
-   return minimalprojections(G, T)
-end
-
 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)
@@ -163,7 +158,7 @@ function rankOne_projections(BN::WreathProduct, T::Type=Rational{Int})
       last_emb = g->BN(Nemo.Generic.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 +182,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}()
+   result = Vector{T}()
-    seen = Set{T}()
+   seen = Set{T}()
-    for x in X
+   for x in X
-        for y in Y
+      for y in Y
-            z = op(x,y)
+         z = op(x,y)
-            if !in(z, seen)
+         if !in(z, seen)
-                push!(seen, z)
+            push!(seen, z)
-                push!(result, z)
+            push!(result, z)
-            end
+         end
-        end
+      end
-    end
+   end
-    return result
+   return result
 end
 
 doc"""

src/PropertyT.jl

@@ -26,36 +26,69 @@ 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
 
-function pmΔfilenames(name::String)
+function pmΔfilenames(prefix::String)
-    if !isdir(name)
+   isdir(prefix) || mkdir(prefix)
-        mkdir(name)
+   pm_filename = joinpath(prefix, "pm.jld")
-    end
+   Δ_coeff_filename = joinpath(prefix, "delta.jld")
-    prefix = name
+   return pm_filename, Δ_coeff_filename
-    pm_filename = joinpath(prefix, "pm.jld")
-    Δ_coeff_filename = joinpath(prefix, "delta.jld")
-    return pm_filename, Δ_coeff_filename
 end
 
-function λSDPfilenames(name::String)
+function λSDPfilenames(prefix::String)
-    if !isdir(name)
+   isdir(prefix) || mkdir(prefix)
-        mkdir(name)
+   λ_filename = joinpath(prefix, "lambda.jld")
-    end
+   SDP_filename = joinpath(prefix, "SDPmatrix.jld")
-    prefix = name
+   return λ_filename, SDP_filename
-    λ_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)
+    sdp_constraints = constraints(product_matrix)
 
     RG = GroupRing(G, product_matrix)
     Δ = GroupRingElem(load(Δ_fname, "Δ")[:, 1], RG)
@@ -73,46 +106,21 @@ 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, "Generating balls of sizes $sizes")
-    info(logger, "Generated balls of sizes $sizes")
+    @logtime logger E_R, sizes = Groups.generate_balls(S, Id, radius=2*radius)
 
     info(logger, "Creating product matrix...")
-    t = @timed 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, timed_msg(t))
 
     info(logger, "Creating sdp_constratints...")
-    t = @timed sdp_constraints = PropertyT.constraints_from_pm(pm)
+    @logtime logger sdp_constraints = PropertyT.constraints(pm)
-    info(logger, timed_msg(t))
 
-    RG = GroupRing(parent(Id), B, pm)
+    RG = GroupRing(parent(Id), E_R, pm)
 
-    Δ = splaplacian(RG, S, Id)
+    Δ = splaplacian(RG, S)
     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]
-    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)
@@ -148,7 +156,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
 
@@ -188,15 +196,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...")
-
+      SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound)
-      t = @timed SDP_problem, λ, P = create_SDP_problem(Δ, sdp_constraints, upper_bound=upper_bound)
-      info(logger, timed_msg(t))
-
       JuMP.setsolver(SDP_problem, solver)
 
       λ, P = λandP(name, SDP_problem, λ, P)
@@ -208,13 +211,16 @@ 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!")
      #  @assert P == Symmetric(P)
-      Q = real(sqrtm(Symmetric(P)))
+      @logtime logger 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

src/SDPs.jl

@@ -1,30 +1,30 @@
 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
 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
@@ -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)
+    @logtime logger solution_status = JuMP.solve(SDP_problem)
-    info(logger, timed_msg(t))
     Base.Libc.flush_cstdio()
 
     redirect_stdout(o)