From 33b45dd51832cdc5e96b58ad8e04cd176d57b3d2 Mon Sep 17 00:00:00 2001
From: kalmar <kalmar@amu.edu.pl>
Date: Wed, 26 Jul 2017 10:29:11 +0200
Subject: [PATCH 1/3] merge two sparsify functions

---
 src/Orbit-wise.jl         | 21 +++++++++++++--------
 src/OrbitDecomposition.jl | 15 +--------------
 2 files changed, 14 insertions(+), 22 deletions(-)

diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl
index a305db4..afc37b8 100644
--- a/src/Orbit-wise.jl
+++ b/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)
diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl
index b8e7eb7..d2963b9 100644
--- a/src/OrbitDecomposition.jl
+++ b/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

From 3ff649683a031f5cfab525acdf2f83b10e9e80cd Mon Sep 17 00:00:00 2001
From: kalmar <kalmar@amu.edu.pl>
Date: Wed, 26 Jul 2017 12:12:58 +0200
Subject: [PATCH 2/3] FIXXX: clamp to 0 only those eps-close to 0, not < eps

---
 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 afc37b8..3590b22 100644
--- a/src/Orbit-wise.jl
+++ b/src/Orbit-wise.jl
@@ -58,7 +58,7 @@ dens(M::AbstractArray) = sum(abs.(M) .!= 0)/length(M)
 
 function sparsify{T}(U::AbstractArray{T}, check=true)
     W = deepcopy(U)
-    W[W .< eps(T)] .= zero(T)
+    W[abs.(W) .< eps(T)] .= zero(T)
     if check && rank(W) != rank(U)
         info("Sparsification would decrease the rank!")
         W = U

From 66f860bac10bc5887cf0efa04cab7934d0bd0e53 Mon Sep 17 00:00:00 2001
From: kalmar <kalmar@amu.edu.pl>
Date: Wed, 26 Jul 2017 12:15:31 +0200
Subject: [PATCH 3/3] sparsify both on the forward and backward transform

---
 src/Orbit-wise.jl         | 2 +-
 src/OrbitDecomposition.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl
index 3590b22..bf85296 100644
--- a/src/Orbit-wise.jl
+++ b/src/Orbit-wise.jl
@@ -83,7 +83,7 @@ function init_orbit_data(logger, sett::Settings; radius=2)
    return 0
 end
 
-function transform(U::AbstractArray, V::AbstractArray; sparse=false)
+function transform(U::AbstractArray, V::AbstractArray; sparse=true)
    if sparse
       return sparsify(U'*V*U)
    else
diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl
index d2963b9..30b7ec3 100644
--- a/src/OrbitDecomposition.jl
+++ b/src/OrbitDecomposition.jl
@@ -123,7 +123,7 @@ function reconstruct_sol{T<:GroupElem, S<:AbstractArray}(mreps::Dict{T, S},
     for g in keys(mreps)
         A, B = mreps[g], mreps[inv(g)]
         for π in 1:length(Us)
-            recP .+= dims[π].* (A * Us[π]*Ps[π]*Ust[π] * B)
+            recP .+= sparsify(dims[π].* (A * Us[π]*Ps[π]*Ust[π] * B))
         end
     end
     recP .*= 1/length(keys(mreps))