Update matfaust.fact.eigtj as previous commit did for pyfaust, handle issue...

Update matfaust.fact.eigtj as previous commit did for pyfaust, handle issue with absolute and relative errors.

Errors are computed with squared norms in C++ but we want the stopping criterion based on error to be specified with simple frob. norms in API of wrappers.

misc/test/src/Python/test_FaustPy.py

@@ -1003,6 +1003,32 @@ class TestFaustFactory(unittest.TestCase):
         # misc/test/src/C++/GivensFGFTParallel.cpp.in
         self.assertEqual(err, err2)
 
+    def testeigtj_abserr(self):
+        from numpy.random import rand
+        from pyfaust.fact import eigtj
+        err = .01
+        M = rand(128,128)
+        M = M@M.T
+        D,U = eigtj(M, tol=err, relerr=False)
+        self.assertAlmostEqual(norm(M-U*np.diag(D)*U.T), err, places=3 )
+        M_cplx = M*np.complex(1,0) + rand(128,128)*np.complex(0,1)
+        M_cplx = M_cplx@np.matrix(M_cplx,copy=False).H
+        err=.1
+        D,U = eigtj(M_cplx, tol=err)
+        self.assertLessEqual(norm(M_cplx-U*np.diag(D)*U.H), err)
+
+    def testeigtj_relerr(self):
+        from numpy.random import rand
+        from pyfaust.fact import eigtj
+        err = .01
+        M = rand(128,128)
+        M = M@M.T
+        D,U = eigtj(M, tol=err)
+        self.assertLessEqual(norm(M-U*np.diag(D)*U.T)/norm(M), err)
+        M_cplx = M + rand(128,128)*np.complex(0,1)
+        M_cplx = M_cplx@np.matrix(M_cplx, copy=False).H
+        D,U = eigtj(M_cplx, tol=err)
+        self.assertLessEqual(norm(M_cplx-U*np.diag(D)*U.H)/norm(M_cplx), err)
 
     def testFactPalm4MSA_fgft(self):
         print("Test pyfaust.fact._palm4msa_fgft()")

src/algorithm/factorization/faust_GivensFGFT.h

@@ -6,6 +6,7 @@
 #include "faust_MatSparse.h"
 #include "faust_MatDense.h"
 #include "faust_Transform.h"
+#include <cfloat>
 #include <vector>
 
 namespace Faust {

src/algorithm/factorization/faust_GivensFGFT.hpp

@@ -388,7 +388,11 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_err()
 		err = Faust::fabs(err_d - err);
 		if(errIsRel) err /= err_d;
 		if(verbosity)
-			cout << "GivensFGFT factor : "<< ite <<  ", transform " << ((errIsRel)?"relative ":"absolute ") << "err.: " << err << endl;
+		{
+			cout << "factor : "<< ite <<  ", " << ((errIsRel)?"relative ":"absolute ") << "err.: " << err;
+			if(stoppingCritIsError) cout << " stoppingError: " << stoppingError << ")";
+			cout << endl;
+		}
 		errs.push_back(err);
 	}
 }
@@ -440,11 +444,14 @@ void GivensFGFT<FPP,DEVICE,FPP2>::compute_facts()
 {
 	is_D_ordered = false; // facts (re)computed then D must be reordered
 	ite = 0;
+	bool stopping = false;
 	while(J == 0 || ite < facts.size()) // when J == 0 the stopping criterion is the error against Lap
 	{
 		next_step();
 		ite++;
-		if(stoppingCritIsError && ((*(errs.end()-1) <= stoppingError || ite > 1 && errs[ite-1]-errs[ite-2] > 0)))
+		if(stopping = (ite > 1 && errs.size() > 2 && errs[ite-1]-errs[ite-2] > FLT_EPSILON))
+			if(verbosity>0) cerr << "warning: the algorithm stopped because the last error is greater than the previous one." << endl;
+		if(stopping || stoppingCritIsError && errs.size() > 0 && (*(errs.end()-1) - stoppingError ) < FLT_EPSILON)
 		{
 			facts.resize(ite);
 			break;

src/algorithm/factorization/faust_GivensFGFTComplex.h

@@ -6,6 +6,7 @@
 #include "faust_MatSparse.h"
 #include "faust_MatDense.h"
 #include "faust_Transform.h"
+#include <cfloat>
 #include <vector>
 
 namespace Faust {

src/algorithm/factorization/faust_GivensFGFTComplex.hpp

@@ -453,7 +453,11 @@ void GivensFGFTComplex<FPP,DEVICE,FPP2>::update_err()
 		err = Faust::fabs(err_d - err);
 		if(errIsRel) err /= err_d;
 		if(verbosity)
-			cout << "GivensFGFTComplex factor : "<< ite <<  ", transform " << ((errIsRel)?"relative ":"absolute ") << "err.: " << err << endl;
+		{
+			cout << "factor : "<< ite <<  ", " << ((errIsRel)?"relative ":"absolute ") << "err.: " << err;
+			if(stoppingCritIsError) cout << " stoppingError: " << stoppingError << ")";
+			cout << endl;
+		}
 		errs.push_back(err);
 	}
 }
@@ -506,11 +510,14 @@ void GivensFGFTComplex<FPP,DEVICE,FPP2>::compute_facts()
 {
 	is_D_ordered = false; // facts (re)computed then D must be reordered
 	ite = 0;
+	bool stopping = false;
 	while(J == 0 || ite < facts.size()) // when J == 0 the stopping criterion is the error against Lap
 	{
 		next_step();
 		ite++;
-		if(stoppingCritIsError && ((*(errs.end()-1) <= stoppingError || ite > 1 && errs[ite-1]-errs[ite-2] > 0)))
+		if(stopping = (ite > 1 && errs.size() > 2 && errs[ite-1]-errs[ite-2] > FLT_EPSILON))
+			if(verbosity>0) cerr << "warning: the algorithm stopped because the last error is greater than the previous one." << endl;
+		if(stopping || stoppingCritIsError && errs.size() > 0 && (*(errs.end()-1) - stoppingError ) < FLT_EPSILON)
 		{
 			facts.resize(ite);
 			break;

src/algorithm/factorization/faust_GivensFGFTParallel.hpp

@@ -166,6 +166,7 @@ void GivensFGFTParallel<FPP,DEVICE,FPP2>::update_fact()
 	this->fact_mod_row_ids.push_back(this->q);
 	this->fact_mod_col_ids.push_back(this->q);
 	this->fact_mod_values.push_back(cos(this->theta));
+	if(this->J == 0) this->facts.resize(this->ite+1);
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>

src/algorithm/factorization/faust_GivensFGFTParallelComplex.hpp

@@ -200,6 +200,7 @@ void GivensFGFTParallelComplex<FPP,DEVICE,FPP2>::update_fact()
 	this->fact_mod_row_ids.push_back(this->q);
 	this->fact_mod_col_ids.push_back(this->q);
 	this->fact_mod_values.push_back(c_qq);
+	if(this->J == 0) this->facts.resize(this->ite+1);
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>

wrapper/matlab/+matfaust/+fact/eigtj.m

@@ -52,7 +52,83 @@
 %> <p> @b See @b also fact.fgft_givens, fact.fgft_palm
 %>
 %==========================================================================================
-function [V,D] = eigtj(M, maxiter, varargin)
-	[V, D] = matfaust.fact.fgft_givens(M, maxiter, varargin{:});
-	%V = factors(V,1:numfactors(V)) % tmp workaround
+function [V,D] = eigtj(M, varargin)
+	import matfaust.Faust
+	nGivens_per_fac = 1; % default value
+	verbosity = 0; % default value
+%	if(~ ismatrix(M) || ~ isreal(M))
+%		error('M must be a real matrix.')
+%	end
+	if(size(M,1) ~= size(M,2))
+		error('M must be square')
+	end
+	bad_arg_err = 'bad number of arguments.';
+	maxiter = 0;
+	tol = 0;
+	relerr = true;
+	verbosity = 0;
+	argc = length(varargin);
+	order = 1; % ascending order
+	if(argc > 0)
+		for i=1:argc
+			switch(varargin{i})
+				case 'maxiter'
+					maxiter = varargin{i+1};
+					if(argc == i || ~ isnumeric(maxiter) || maxiter-floor(maxiter) > 0 || maxiter <= 0 )
+						error('maxiter keyword arg. is not followed by an integer')
+					end
+				case 'tol'
+					if(argc == i || ~ isscalar(varargin{i+1}))
+						error('tol keyword arg. is not followed by a number')
+					else
+						tol = real(varargin{i+1}); % real in case of cplx num
+					end
+				case 'relerr'
+					if(argc == i || ~ islogical(varargin{i+1}))
+						error('relerr keyword argument is not followed by a logical')
+					else
+						relerr = varargin{i+1};
+					end
+				case 'verbosity'
+					if(argc == i || ~ isscalar(varargin{i+1}))
+						error('verbose keyword argument is not followed by a number')
+					else
+						verbosity = floor(real(varargin{i+1}));
+					end
+				case 'nGivens_per_fac'
+					if(argc == i || ~ isscalar(varargin{i+1}) || ~ isnumeric(varargin{i+1}))
+						error('nGivens_per_fac must be followed by a positive integer.')
+					else
+						nGivens_per_fac = floor(abs(real(varargin{i+1})));
+						nGivens_per_fac = max(1, nGivens_per_fac);
+					end
+				case 'order'
+					if(argc == i || (~ strcmp(varargin{i+1}, 'ascend') && ~ strcmp(varargin{i+1}, 'descend') && ~ strcmp(varargin{i+1}, 'undef')))
+						error('order must be followed by a char array among ''ascend'', ''descend'' or ''undef''.')
+					else
+						order = varargin{i+1};
+						if(order(1) == 'a')
+							order = 1
+						elseif(order(1) == 'd')
+							order = -1
+						else
+							order = 0
+						end
+					end
+				otherwise
+					if(isstr(varargin{i}) && (~ strcmp(varargin{i}, 'ascend') && ~ strcmp(varargin{i}, 'descend') && ~ strcmp(varargin{i}, 'undef')) )
+						error([ varargin{i} ' unrecognized argument'])
+					end
+			end
+		end
+	end
+	if(maxiter == 0 && tol == 0)
+		error('Either maxiter or tol must be greater than zero.')
+	end
+	if(maxiter > 0)
+		nGivens_per_fac = min(nGivens_per_fac, maxiter);
+	end
+	[core_obj, D] = mexfgftgivensReal(M, maxiter, nGivens_per_fac, verbosity, tol, relerr, order);
+	D = sparse(diag(real(D)));
+	V = Faust(core_obj, isreal(M));
 end

wrapper/matlab/+matfaust/+fact/fgft_givens.m

+% experimental block start
 %==========================================================================================
 %> @brief Computes the FGFT of the Laplacian matrix Lap (using fact.eigtj).
 %>
@@ -29,74 +30,6 @@
 %>
 %==========================================================================================
 function [FGFT,D] = fgft_givens(Lap, maxiter, varargin)
-	import matfaust.Faust
-	nGivens_per_fac = 1; % default value
-	verbosity = 0; % default value
-%	if(~ ismatrix(Lap) || ~ isreal(Lap))
-%		error('Lap must be a real matrix.')
-%	end
-	if(size(Lap,1) ~= size(Lap,2))
-		error('Lap must be square')
-	end
-	if(~ isnumeric(maxiter) || maxiter-floor(maxiter) > 0 || maxiter <= 0)
-		error('maxiter must be a positive integer.')
-	end
-	bad_arg_err = 'bad number of arguments.';
-	tol = 0;
-	relerr = true;
-	verbosity = 0;
-	argc = length(varargin);
-	order = 1; % ascending order
-	if(argc > 0)
-		for i=1:argc
-			switch(varargin{i})
-				case 'tol'
-					if(argc == i || ~ isscalar(varargin{i+1}))
-						error('tol keyword arg. is not followed by a number')
-					else
-						tol = real(varargin{i+1}); % real in case of cplx num
-					end
-				case 'relerr'
-					if(argc == i || ~ islogical(varargin{i+1}))
-						error('relerr keyword argument is not followed by a logical')
-					else
-						relerr = varargin{i+1};
-					end
-				case 'verbosity'
-					if(argc == i || ~ isscalar(varargin{i+1}))
-						error('verbose keyword argument is not followed by a number')
-					else
-						verbosity = floor(real(varargin{i+1}));
-					end
-				case 'nGivens_per_fac'
-					if(argc == i || ~ isscalar(varargin{i+1}) || ~ isnumeric(varargin{i+1}))
-						error('nGivens_per_fac must be followed by a positive integer.')
-					else
-						nGivens_per_fac = floor(abs(real(varargin{i+1})));
-						nGivens_per_fac = min(nGivens_per_fac, maxiter);
-						nGivens_per_fac = max(1, nGivens_per_fac);
-					end
-				case 'order'
-					if(argc == i || (~ strcmp(varargin{i+1}, 'ascend') && ~ strcmp(varargin{i+1}, 'descend') && ~ strcmp(varargin{i+1}, 'undef')))
-						error('order must be followed by a char array among ''ascend'', ''descend'' or ''undef''.')
-					else
-						order = varargin{i+1};
-						if(order(1) == 'a')
-							order = 1
-						elseif(order(1) == 'd')
-							order = -1
-						else
-							order = 0
-						end
-					end
-				otherwise
-					if(isstr(varargin{i}) && (~ strcmp(varargin{i}, 'ascend') && ~ strcmp(varargin{i}, 'descend') && ~ strcmp(varargin{i}, 'undef')) )
-						error([ varargin{i} ' unrecognized argument'])
-					end
-			end
-		end
-	end
-	[core_obj, D] = mexfgftgivensReal(Lap, maxiter, nGivens_per_fac, verbosity, tol, relerr, order);
-	D = sparse(diag(real(D)));
-	FGFT = Faust(core_obj, isreal(Lap));
+	[FGFT, D] = matfaust.fact.eigtj(Lap, maxiter, varargin{:});
 end
+% experimental block end

wrapper/matlab/src/mexfgftgivens.cpp.in

@@ -84,6 +84,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 	if(nrhs >= 7)
 		order = (int) mxGetScalar(prhs[6]); //eigenvalues order
 
+	tol *= tol; // C++ backend works with squared norm error
+
 	const mxArray* matlab_matrix = prhs[0]; // Laplacian
 
 	if(mxIsComplex(matlab_matrix))

wrapper/python/pyfaust/fact.py

@@ -260,7 +260,7 @@ def fgft_givens(Lap, maxiter, tol=0.0, relerr=True, nGivens_per_fac=None,
     See also:
         eigtj, fgft_palm
     """
-    return eigtj(M, maxiter, tol, relerr, nGivens_per_fac, verbosity, order)
+    return eigtj(Lap, maxiter, tol, relerr, nGivens_per_fac, verbosity, order)
 # experimental block end
 
 def _check_fact_mat(funcname, M):

wrapper/python/src/_FaustCorePy.pyx

@@ -1584,12 +1584,13 @@ cdef class FaustFact:
                         " (to define a stopping  criterion)")
             maxiter = 0 
         if(nGivens_per_fac == None): nGivens_per_fac = int(M.shape[0]/2)
-        if((isinstance(M, np.ndarray) and (np.matrix(M, copy=False).H != M).any()) or M.shape[0] != M.shape[1]):
+        if(isinstance(M, np.ndarray) and \
+            not np.allclose(np.matrix(M, copy=False).H, M) or M.shape[0] != M.shape[1]):
             raise ValueError(" the matrix/array must be symmetric or hermitian.")
         if(not isinstance(maxiter, int)): raise TypeError("maxiter must be a int")
         if(not isinstance(nGivens_per_fac, int)): raise TypeError("nGivens_per_fac must be a int")
         nGivens_per_fac = max(nGivens_per_fac, 1)
-        nGivens_per_fac = min(nGivens_per_fac, maxiter)
+        if(maxiter > 0): nGivens_per_fac = min(nGivens_per_fac, maxiter)
         tol *= tol # the C++ impl. works on squared norms to measure errors
 
         M_is_real = M.dtype in [ 'float', 'float128',