From 7dbd11bebfbfcef5fb32c322d47dc082c1b6c580 Mon Sep 17 00:00:00 2001
From: kalmarek <kalmar@amu.edu.pl>
Date: Sun, 27 Aug 2017 18:58:01 +0200
Subject: [PATCH] eps-sparsify for both sparse and dense matrices

---
 src/Orbit-wise.jl | 44 +++++++++++++++++++++++++++++++++-----------
 1 file changed, 33 insertions(+), 11 deletions(-)

diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl
index 0fb6fd5..00b2acf 100644
--- a/src/Orbit-wise.jl
+++ b/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)
-    W = deepcopy(U)
-    W[abs.(W) .< eps(T)] .= zero(T)
-    if check && rank(W) != rank(U)
-        info("Sparsification would decrease the rank!")
-        W = U
-    else
-      info("Sparsified:", rpad(dens(U), 10), "\t", rpad(dens(W),10))
+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
-    W = sparse(W)
-    dropzeros!(W)
-    return W
+   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)
+   else
+      W = U
+   end
+
+   W[abs.(W) .< eps] .= zero(T)
+
+   if check
+      info("Sparsification would decrease the rank!")
+      W = U
+   else
+      info("Sparsified density:", rpad(dens(U), 15), "→", rpad(dens(W),15))
+   end
+   W = sparse(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))