merge two sparsify functions

src/Orbit-wise.jl

@@ -53,13 +53,18 @@ end
 
 include("OrbitDecomposition.jl")
 
-function sparsify{T}(U::AbstractArray{T}, eps=eps(T), check=true)
+dens(M::SparseMatrixCSC) = length(M.nzval)/length(M)
+dens(M::AbstractArray) = sum(abs.(M) .!= 0)/length(M)
+
+function sparsify{T}(U::AbstractArray{T}, check=true)
     W = deepcopy(U)
-    W[abs.(W) .< eps] = zero(T)
+    W[W .< eps(T)] .= zero(T)
     if check && rank(W) != rank(U)
-        warn("Sparsification would decrease the rank!")
+        info("Sparsification would decrease the rank!")
         W = U
-    end
+    else
+      info("Sparsified:", rpad(dens(U), 10), "\t", rpad(dens(W),10))
+   end
     W = sparse(W)
     dropzeros!(W)
     return W
@@ -79,14 +84,14 @@ function init_orbit_data(logger, sett::Settings; radius=2)
 end
 
 function transform(U::AbstractArray, V::AbstractArray; sparse=false)
-   w = U'*V*U
    if sparse
-      w = sparsify(w)
+      return sparsify(U'*V*U)
+   else
+      return U'*V*U
    end
-   return w
 end
 
-A(data::OrbitData, π, t) = data.dims[π]*transform(data.Us[π], data.cnstr[t], sparse=true)
+A(data::OrbitData, π, t) = data.dims[π]*transform(data.Us[π], data.cnstr[t])
 
 function constrLHS(m::JuMP.Model, data::OrbitData, t)
     l = endof(data.Us)

src/OrbitDecomposition.jl

@@ -140,19 +140,6 @@ function Cstar_repr(x::GroupRingElem, mreps::Dict)
     return res
 end
 
-dens(M::SparseMatrixCSC) = length(M.nzval)/length(M)
-dens(M::AbstractArray) = sum(abs.(M) .!= 0)/length(M)
-
-function sparsify2(M::AbstractArray)
-   println("Density before sparsification: \t$(dens(M))")
-   M = deepcopy(M)
-   M[M .< eps(eltype(M))] .= zero(eltype(M))
-   M = sparse(M)
-   dropzeros!(M)
-   println("Density after sparsification: \t$(dens(M))")
-   return M
-end
-
 function orthSVD(M::AbstractMatrix)
     M = full(M)
     fact = svdfact(M)
@@ -205,6 +192,6 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S
    info(logger, "dimensions = $dimensions")
    @assert dot(multiplicities, dimensions) == sizes[radius]
 
-   save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, "spUπs", sparsify2.(Uπs), "dims", dimensions)
+   save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, "spUπs", sparsify.(Uπs), "dims", dimensions)
    return 0
 end