Update pyfaust eigtj (and fgft_givens) to respect the conventions established by numpy.linalg.eigh.

Returning a vector array for the eigenvalues.
Change the order of returned tuple elements.
Updating the doc and tests.

misc/test/src/Python/test_FaustPy.py

@@ -7,6 +7,7 @@ import sys
 import numpy as np
 from scipy.io import savemat,loadmat
 from numpy.linalg import norm
+from scipy.sparse import spdiags
 import math
 
 class TestFaustPy(unittest.TestCase):
@@ -906,7 +907,8 @@ class TestFaustFactory(unittest.TestCase):
         L = L.astype(np.float64)
         J = \
         int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
-        F, D = pyfaust.fact.fgft_givens(L, J, nGivens_per_fac=1, verbosity=0)
+        D, F = pyfaust.fact.fgft_givens(L, J, nGivens_per_fac=1, verbosity=0)
+        D = spdiags(D, [0], L.shape[0], L.shape[0])
         print("Lap norm:", norm(L, 'fro'))
         err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
         print("err: ", err)
@@ -922,14 +924,16 @@ class TestFaustFactory(unittest.TestCase):
         J = \
                 int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
         t = int(L.shape[0]/2)
-        F, D = fgft_givens(L, J, nGivens_per_fac=t, verbosity=0)
+        D, F = fgft_givens(L, J, nGivens_per_fac=t, verbosity=0)
+        D = spdiags(D, [0], L.shape[0], L.shape[0])
         print("Lap norm:", norm(L, 'fro'))
         err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
         print("err: ", err)
         # the error reference is from the C++ test,
         # misc/test/src/C++/GivensFGFTParallel.cpp.in
         self.assertAlmostEqual(err, 0.0410448, places=7)
-        F2, D2 = eigtj(L, J, nGivens_per_fac=t, verbosity=0)
+        D2, F2 = eigtj(L, J, nGivens_per_fac=t, verbosity=0)
+        D2 = spdiags(D, [0], L.shape[0], L.shape[0])
         print("Lap norm:", norm(L, 'fro'))
         err2 = norm((F2*D.todense())*F2.T.todense()-L,"fro")/norm(L,"fro")
         print("err2: ", err2)
@@ -946,7 +950,8 @@ class TestFaustFactory(unittest.TestCase):
         L = L.astype(np.float64)
         J = \
         int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
-        F, D = pyfaust.fact.fgft_givens(csr_matrix(L), J, nGivens_per_fac=0, verbosity=0)
+        D, F = pyfaust.fact.fgft_givens(csr_matrix(L), J, nGivens_per_fac=0, verbosity=0)
+        D = spdiags(D, [0], L.shape[0], L.shape[0])
         print("Lap norm:", norm(L, 'fro'))
         err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
         print("err: ", err)
@@ -963,14 +968,16 @@ class TestFaustFactory(unittest.TestCase):
         J = \
                 int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
         t = int(L.shape[0]/2)
-        F, D = fgft_givens(csr_matrix(L), J, nGivens_per_fac=t, verbosity=0)
+        D, F = fgft_givens(csr_matrix(L), J, nGivens_per_fac=t, verbosity=0)
+        D = spdiags(D, [0], L.shape[0], L.shape[0])
         print("Lap norm:", norm(L, 'fro'))
         err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
         print("err: ", err)
         # the error reference is from the C++ test,
         # misc/test/src/C++/GivensFGFTParallel.cpp.in
         self.assertAlmostEqual(err, 0.0410448, places=7)
-        F2, D2 = eigtj(L, J, nGivens_per_fac=t, verbosity=0)
+        D2, F2 = eigtj(L, J, nGivens_per_fac=t, verbosity=0)
+        D2 = spdiags(D, [0], L.shape[0], L.shape[0])
         print("Lap norm:", norm(L, 'fro'))
         err2 = norm((F2*D.todense())*F2.T.todense()-L,"fro")/norm(L,"fro")
         print("err2: ", err2)

wrapper/python/pyfaust/fact.py

@@ -80,9 +80,9 @@ def svdtj(M, J, t=1):
 
 def eigtj(M, maxiter, tol=0, relerr=True,  nGivens_per_fac=None, verbosity=0):
     """
-    Runs the truncated Jacobi algorithm to compute the eigenvalues of M (returned in W) and the corresponding eigenvectors (in Faust V) using the truncated Jacobi algorithm.
+    Runs the truncated Jacobi algorithm to compute the eigenvalues of M (returned in W) and the corresponding transform of eigenvectors (in Faust V columns).
 
-    The output is such that V*W.todense()*V.T approximates M.
+    The output is such that dot(dot(V,diag(W)), V.T) approximates M.
 
     The trade-off between accuracy and sparsity can be set through the
     parameters maxiter and nGivens_per_fac.
@@ -104,13 +104,14 @@ def eigtj(M, maxiter, tol=0, relerr=True,  nGivens_per_fac=None, verbosity=0):
 
 
     Returns:
-        The tuple (V,W):
-            - V the Faust object representing the approximate eigenvector
-            transform. The last factor of V is a permutation matrix.
+        The tuple (W,V):
+           - W (numpy.ndarray) the vector of the eigenvalues of M (in ascendant order).
+           - V the Faust object representing the approximate eigenvector
+            transform. The column V[:, i] is the eigenvector
+            corresponding to the eigenvalue W[i].<br/>
+            The last factor of V is a permutation matrix.
             The goal of this factor is to apply to the columns of V the same order as
             eigenvalues set in W.
-            - W (scipy.sparse.dia.dia_matrix) the diagonal matrix of the
-            approximate eigenvalues (in ascendant order along the rows/columns).
 
     References:
     [1]   Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast
@@ -129,7 +130,7 @@ def eigtj(M, maxiter, tol=0, relerr=True,  nGivens_per_fac=None, verbosity=0):
         demo_path = sep.join((get_data_dirpath(),'Laplacian_256_community.mat'))
         data_dict = loadmat(demo_path)
         Lap = data_dict['Lap'].astype(np.float)
-        Uhat, Dhat = eigtj(Lap, maxiter=Lap.shape[0]*100)
+        Dhat, Uhat = eigtj(Lap, maxiter=Lap.shape[0]*100)
         # Uhat is the Fourier matrix/eigenvectors approximation as a Faust
         # (200 factors + permutation mat.)
         # Dhat the eigenvalues diagonal matrix approx.
@@ -154,11 +155,11 @@ def fgft_givens(Lap, maxiter, tol=0.0, relerr=True, nGivens_per_fac=None, verbos
         verbosity: see eigtj
 
     Returns:
-        The tuple (FGFT,D):
+        The tuple (D, FGFT):
         - with FGFT being the Faust object representing
         the Fourier transform and,
-        - (scipy.sparse.dia.dia_matrix) D the eigenvalues of the Laplacian.
-        The eigenvalues are in ascendant order along the rows/columns.
+        - (numpy.ndarray) D a vector of the eigenvalues of the Laplacian (in
+        ascendant order).
 
 
     References:
@@ -168,8 +169,10 @@ def fgft_givens(Lap, maxiter, tol=0.0, relerr=True, nGivens_per_fac=None, verbos
     over Networks 2018, 4(2), pp 407-420.
     <https://hal.inria.fr/hal-01416110>
 
-    <b/> See also eigtj, fgft_palm
+    See also:
+        eigtj, fgft_palm
     """
+    if(nGivens_per_fac == None): nGivens_per_fac = int(Lap.shape[0]/2)
     if((isinstance(Lap, np.ndarray) and (Lap.T != Lap).any()) or Lap.shape[0] != Lap.shape[1]):
         raise ValueError(" the matrix/array must be square and should be symmetric.")
     if(not isinstance(maxiter, int)): raise TypeError("maxiter must be a int")
@@ -190,7 +193,7 @@ def fgft_givens(Lap, maxiter, tol=0.0, relerr=True, nGivens_per_fac=None, verbos
     else:
         raise TypeError("The matrix to diagonalize must be a"
                         " scipy.sparse.csr_matrix or a numpy array.")
-    return Faust(core_obj=core_obj), D
+    return D, Faust(core_obj=core_obj)
 
 def _check_fact_mat(funcname, M):
     if(not isinstance(M, np.ndarray)):

wrapper/python/src/_FaustCorePy.pyx

@@ -1381,8 +1381,9 @@ cdef class FaustFact:
                                                                    errIsRel)
 
         core._isReal = True
-        D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
-        return core, D_spdiag
+        #D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
+        #return core, D_spdiag
+        return core, D
 
     @staticmethod
     def fact_givens_fgft(Lap, J, t, verbosity=0, stoppingError = 0.0,
@@ -1415,6 +1416,7 @@ cdef class FaustFact:
                                                                            errIsRel)
 
         core._isReal = True
-        from scipy.sparse import spdiags
-        D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
-        return core, D_spdiag
+        #from scipy.sparse import spdiags
+        #D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
+        #return core, D_spdiag
+        return core, D