eps-sparsify for both sparse and dense matrices

src/Orbit-wise.jl

@@ -57,20 +57,42 @@ include("OrbitDecomposition.jl")
 dens(M::SparseMatrixCSC) = length(M.nzval)/length(M)
 dens(M::AbstractArray) = length(findn(M)[1])/length(M)
 
-function sparsify{T}(U::AbstractArray{T}, check=true)
+function sparsify!{Tv,Ti}(M::SparseMatrixCSC{Tv,Ti}, eps=eps(Tv))
+   n = nnz(M)
+   for i in eachindex(M.nzval)
+      if abs(M.nzval[i]) < eps
+         M.nzval[i] = zero(Tv)
+      end
+   end
+   dropzeros!(M)
+   m = nnz(M)
+
+   info("Sparsified density:", rpad(dens(U), 15), "→", rpad(dens(W),15))
+
+   return M
+end
+
+function sparsify!{T}(U::AbstractArray{T}, eps=eps(T); check=true)
+   if check
       W = deepcopy(U)
-    W[abs.(W) .< eps(T)] .= zero(T)
+   else
-    if check && rank(W) != rank(U)
+      W = U
+   end
+
+   W[abs.(W) .< eps] .= zero(T)
+
+   if check
       info("Sparsification would decrease the rank!")
       W = U
    else
-      info("Sparsified:", rpad(dens(U), 10), "\t", rpad(dens(W),10))
+      info("Sparsified density:", rpad(dens(U), 15), "→", rpad(dens(W),15))
    end
    W = sparse(W)
-    dropzeros!(W)
    return W
 end
 
+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))