Implement/document matfaust.tools.omp for equivalency to pyfaust.tools.omp.

Also adjusting the pyfaust.demo and tools because of omp (documentation and arguments).

gen_doc/CMakeLists.txt

@@ -7,7 +7,7 @@ if(BUILD_DOCUMENTATION)
 			string(CONCAT DOXYGEN_FILE_PATTERNS "*.cpp *.hpp *.h *.cu *.hu")
 		endif()
 		if(BUILD_WRAPPER_MATLAB)
-			string(CONCAT DOXYGEN_FILE_PATTERNS ${DOXYGEN_FILE_PATTERNS} " Faust.m StoppingCriterion.m ConstraintGeneric.m ConstraintMat.m ConstraintReal.m ConstraintInt.m ConstraintName.m ParamsFact.m ParamsHierarchicalFact.m ParamsPalm4MSA.m FaustFactory.m hadamard.m quickstart.m fft.m bsl.m runtimecmp.m runall.m version.m faust_fact.m ParamsHierarchicalFactSquareMat.m ParamsHierarchicalFactRectMat.m license.m")
+			string(CONCAT DOXYGEN_FILE_PATTERNS ${DOXYGEN_FILE_PATTERNS} " Faust.m StoppingCriterion.m ConstraintGeneric.m ConstraintMat.m ConstraintReal.m ConstraintInt.m ConstraintName.m ParamsFact.m ParamsHierarchicalFact.m ParamsPalm4MSA.m FaustFactory.m hadamard.m quickstart.m fft.m bsl.m runtimecmp.m runall.m version.m faust_fact.m ParamsHierarchicalFactSquareMat.m ParamsHierarchicalFactRectMat.m license.m omp.m")
 		endif()
 		if(BUILD_WRAPPER_PYTHON)
 			string(CONCAT DOXYGEN_FILE_PATTERNS ${DOXYGEN_FILE_PATTERNS} " __init__.py factparams.py demo.py tools.py")

gen_doc/Doxyfile.in

@@ -519,7 +519,7 @@ EXCLUDE_SYMLINKS       = NO
 
 EXCLUDE_PATTERNS       = *faust_MatDense_gpu.h *faust_MatSparse_gpu.h 
 
-EXCLUDE_SYMBOLS = mtimes_trans subsasgn get_factor_nonopt reshape ABC ParamsPalm4MSAFGFT UpdateCholesky UpdateCholeskyFull UpdateCholeskySparse
+EXCLUDE_SYMBOLS = mtimes_trans subsasgn get_factor_nonopt reshape ABC ParamsPalm4MSAFGFT UpdateCholesky UpdateCholeskyFull UpdateCholeskySparse greed_omp_chol
 
 # The EXAMPLE_PATH tag can be used to specify one or more files or 
 # directories that contain example code fragments that are included (see 

wrapper/matlab/+matfaust/+tools/omp.m

+%> @package matfaust.tools @brief The matfaust tools namespace
+
+%===============================================================================
+%>  Runs the greedy OMP algorithm optimized by Cholesky decomposition.
+%===
+%> @b Usage
+%>
+%> &nbsp;&nbsp;&nbsp; @b x=omp(y, D) runs the algorithm with the default configuration which is equivalent to <b> x=omp(y, D, 'maxiter', length(y), 'verbose', false)</b>.<br/>
+%> &nbsp;&nbsp;&nbsp; <b> x=omp(y, D, 'maxiter', N)</b> stops the algorithm after N iterations.<br/>
+%> &nbsp;&nbsp;&nbsp; <b> x=omp(y, D, 'tol', 10**-16)</b> runs the algoritm until the relative error is lower or equal 10**-16. This is equivalent to <b>x=omp(y, D, 'tol', 10**-16, 'relerr', true)</b><br/>
+%> &nbsp;&nbsp;&nbsp; <b>x=omp(y, D, 'tol', 10**-16, 'relerr', false)</b> runs the algoritm until the absolute error is lower or equal 10**-16.<br/><br/>
+%>
+%>
+%> @param y The vector to approximate by D*x.
+%> @param D The dictionary as a matrix or a Faust.
+%> @param 'verbose', true (optional) To enable the verbosity.
+%> @param 'verbose', false (optional) To disable the verbosity (this is the default option).
+%> @param 'maxiter', number (optional) To define a maximum number of iterations (by default this length(y)).
+%> @param 'relerr', true (optional) For a stopping criterion based on the relative error (this is the default error).
+%> @param 'relerr', false (optional) For a stopping criterion based on the absolute error.
+%> @param 'tol', number (optional) To define the error value used for the relative or absolute error (by default this is 0, for not stopping on any error criterion).
+%>
+%>
+%>@return x the solution of y = D*x (according to the error).
+%===============================================================================
+function x = omp(y, D, varargin)
+	argc = length(varargin);
+	% set parameter default values
+	maxiter = length(y)
+	tol = 0
+	relerr = true
+	verbose = false
+	if(argc > 0)
+		for i=1:argc
+			switch(varargin{i})
+				case 'maxiter'
+					if(argc == i || ~ isscalar(varargin{i+1}))
+						error('maxiter keyword arg. is not followed by a number')
+					else
+						maxiter = real(floor((varargin{i+1}))) % real in case of cplx num
+					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 'verbose'
+					if(argc == i || ~ islogical(varargin{i+1}))
+						error('verbose keyword argument is not followed by a logical')
+					else
+						verbose = varargin{i+1}
+					end
+				otherwise
+					if(isstr(varargin{i}))
+						error('unrecognized argument')
+					end
+			end
+		end
+	end
+	if(tol > 0)
+		tol
+		if(relerr)
+			if(size(y,1) == 1)
+				y_sqr_norm = y*y'
+			else
+				y_sqr_norm = y'*y
+			end
+			x = greed_omp_chol(y, D, size(D,2), 'stopCrit', 'mse',  'stopTol', tol*y_sqr_norm/length(y), 'verbose', verbose);
+		else % absolute error
+			x = greed_omp_chol(y, D, size(D,2), 'stopCrit', 'mse',  'stopTol', tol^2/length(y), 'verbose', verbose);
+		end
+	else
+		% maxiter
+		x = greed_omp_chol(y, D, size(D,2), 'stopCrit', 'M', 'stopTol', maxiter, 'verbose', verbose)
+	end
+end

wrapper/python/pyfaust/demo.py

@@ -1134,8 +1134,9 @@ class bsl:
                     #print(matrix(data[:,i:i+1]).shape,  MEG.shape, M)
                     n = max(matrix(data[:,i:i+1]).shape)
                     solver_sol = greed_omp_chol(matrix(data[:,i:i+1]), MEG,
-                                                maxiter=min(sparsity,floor(n/4)),
+                                                maxiter=sparsity,
                                                 tol=10**-8*np.sqrt(n),
+                                                relerr=False,
                                                 verbose=False)
                     # about tol: it's the square root of 10**-16 (because
                     # the square of tol is used into greed_omp_chol) 

wrapper/python/pyfaust/tools.py

 # -*- coding: utf-8 -*-
 # @PYFAUST_LICENSE_HEADER@
 
-## @package pyfaust @brief <b> The FAµST tools module</b> 
+## @package pyfaust.tools @brief The pyfaust tools module
 
 import numpy as np
 from scipy.sparse.linalg import spsolve_triangular
@@ -27,7 +27,7 @@ def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False):
         verbose: to enable the verbosity (value to True).
 
     Returns:
-        x: the solution of y = D*x.
+        x: the solution of y = D*x (according to the error).
     """
     # check y is a numpy.matrix (or a matrix_csr ?)
     if(isinstance(y, np.ndarray)): y = matrix(y, copy=False)
@@ -56,10 +56,10 @@ def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False):
         if(not (type(maxiter) in [int, float])):
             raise ValueError("maxiter must be a number.")
 
-    tolerr = tol**2
     if(relerr):
-        tolerr *= (y.H*y)[0,0]
-
+        tolerr = tol * (y.H*y)[0,0]
+    else:
+        tolerr = tol**2
 
     # try if enough memory
     try: