diff --git a/misc/test/src/Python/test_FaustPy.py b/misc/test/src/Python/test_FaustPy.py index 038f9428422d5df2f1dad364d9c90d8b495b3bf1..570800fa0bea2548655abe5fdde0640c3b5924de 100644 --- a/misc/test/src/Python/test_FaustPy.py +++ b/misc/test/src/Python/test_FaustPy.py @@ -1191,6 +1191,113 @@ class TestFaustFactory(unittest.TestCase): self.assertLessEqual(count_nonzero(Mp[:,i]), k) self.assertAlmostEqual(norm(Mp), 1) + def test_splincol(self): + from pyfaust.factparams import splincol + from random import randint + from numpy.random import rand + from numpy import count_nonzero + from numpy.linalg import norm + min_n, min_m = 5, 5 + m = randint(min_m, 128) + n = randint(min_n, 128) + M = rand(m,n) + k = randint(1,m) + p = splincol((m,n),k) + Mp = p(M) + # TODO: define sparsity assertions to verify on Mp + self.assertAlmostEqual(norm(Mp), 1) + + def test_sp(self): + from pyfaust.factparams import sp + from random import randint + from numpy.random import rand + from numpy import count_nonzero + from numpy.linalg import norm + min_n, min_m = 5, 5 + m = randint(min_m, 128) + n = randint(min_n, 128) + M = rand(m,n) + k = randint(1,m*n) + p = sp((m,n),k) + Mp = p(M) + for i in range(0,n): + # np.savez('M.npz', M, Mp, k) + self.assertLessEqual(count_nonzero(Mp[:,i]), k) + self.assertAlmostEqual(norm(Mp), 1) + + def test_supp(self): + from pyfaust.factparams import supp + from numpy.random import rand + from numpy.random import randn, permutation as randperm + from numpy.linalg import norm + from random import randint + min_n, min_m = 5, 5 + m = randint(min_m, 128) + n = randint(min_n, 128) + M = rand(m,n) + k = randint(1, min(m,n)) + nnz_rinds = randperm(m)[:k] + nnz_cinds = randperm(n)[:k] + S = np.zeros((m,n)) + S[nnz_rinds, nnz_cinds] = 1 + p = supp(S) + pM = p(M) + # same nnz number + self.assertEqual(np.count_nonzero(pM), np.count_nonzero(S)) + # same support + self.assertTrue(np.allclose(pM != 0, S != 0)) + # pM normalized (according to fro-norm) + self.assertAlmostEqual(norm(pM), 1) + # same projection without normalization + p = supp(S, normalized=False) + pM = p(M) + self.assertTrue(np.allclose(pM[pM != 0], M[S != 0])) + + def test_const(self): + from pyfaust.factparams import const + from numpy.random import rand + from random import randint + min_n, min_m = 5, 5 + m = randint(min_m, 128) + n = randint(min_n, 128) + M = rand(m,n) + C = rand(m,n) + p = const(C, normalized=False) + pM = p(M) + self.assertTrue(np.allclose(C, pM)) + + def test_normcol(self): + from pyfaust.factparams import normcol + from numpy.random import rand + from numpy.random import randn, permutation as randperm + from numpy.linalg import norm + from random import randint, random + min_n, min_m = 5, 5 + m = randint(min_m, 128) + n = randint(min_n, 128) + M = rand(m,n)*random()*50 + k = random()*50 + p = normcol(M.shape,k, normalized=False) + pM = p(M) + for i in range(pM.shape[1]): + self.assertAlmostEqual(norm(pM[:,i]), k) + + def test_normlin(self): + from pyfaust.factparams import normlin + from numpy.random import rand + from numpy.random import randn, permutation as randperm + from numpy.linalg import norm + from random import randint, random + min_n, min_m = 5, 5 + m = randint(min_m, 128) + n = randint(min_n, 128) + M = rand(m,n)*random()*50 + k = random()*50 + p = normlin(M.shape,k, normalized=False) + pM = p(M) + for i in range(pM.shape[0]): + self.assertAlmostEqual(norm(pM[i,:]), k) + if __name__ == "__main__": if(len(sys.argv)> 1): # argv[1] is for adding a directory in PYTHONPATH diff --git a/src/algorithm/factorization/faust_GivensFGFT.h b/src/algorithm/factorization/faust_GivensFGFT.h index 51077efce28fd42c3268a56fa48fc164170f45b5..255100f3243241e557a1f8387ccf67d15af09826 100644 --- a/src/algorithm/factorization/faust_GivensFGFT.h +++ b/src/algorithm/factorization/faust_GivensFGFT.h @@ -21,6 +21,8 @@ namespace Faust { * \brief This class implements the Givens FGFT algorithm. * This algorithm is based on the classical Jacobi eigenvalues algorithm. * + * See parent class GivensFGFTGen for documentation about members. + * * References: * * [1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast @@ -34,73 +36,13 @@ namespace Faust { Faust::MatDense<FPP,DEVICE> C; /** \brief Column vector for the rowwise minimization of C (i.e. maximization of L). */ Faust::Vect<FPP,DEVICE> C_min_row; - /** \brief Pivot candidates q coordinates. */ -// int* q_candidates; /* default IndexType for underlying eigen matrix is int. */ protected: - /** \brief The number of targeted transform factors (when only one Givens rotation is stored per factor).*/ -// int J; - /** \brief Fourier matrix factorization matrices (Givens matrix). */ -// vector<Faust::MatSparse<FPP,DEVICE>> facts; - /** \brief Diagonalization approximate of Laplacian. */ -// Faust::Vect<FPP,DEVICE> D; - /** \brief Queue of errors (cf. calc_err()). */ -// vector<FPP2> errs; - /** \brief Pivot choices (p, q) coordinates. */ -// vector<pair<int,int>> coord_choices; -// FPP Lap_squared_fro_norm; - /** \brief L iteration factor: L_i = S^T L_{i-1} S, initialized from Lap (with S being facts[i]). */ -// Faust::MatGeneric<FPP, DEVICE> *L; /** \brief Rotation angle theta for the current iteration's Givens matrix. */ FPP2 theta; - /* Precomputed model identity matrix to init. facts[ite] before update. - * Identity matrix is completed later with cos/sin coefficients (in update_fact()). - */ - - /** \brief Defines the rows of facts[ite]. */ -// vector<int> fact_mod_row_ids; -// /** \brief Defines the columns of facts[ite]. */ -// vector<int> fact_mod_col_ids; -// /** \brief Defines the coefficients of facts[ite]. */ -// vector<FPP> fact_mod_values; - - - /** \brief Ordered indices of D to get increasing eigenvalues along the diagonal. */ -// vector<int> ord_indices; - - /** \brief inverse permutation of ord_indices (needed to retrieve start undefined order). */ -// vector<int> inv_ord_indices; - /** \brief True is the last fact (of facts) has been permuted */ -// bool last_fact_permuted; - /** \brief Cache for the ordered D. */ -// Faust::Vect<FPP,DEVICE> ordered_D; - /** \brief true if D has already been ordered (order_D() was called). */ -// bool is_D_ordered; -// int D_order_dir; - - /** \brief The level of verbosity (0 for nothing, 1 for iteration numbers,...) */ -// unsigned int verbosity; - - /** - * \brief Row index for the selected pivot in L. - */ -// int p, -// /** -// * \brief Column index for the selected pivot in L. -// */q; - - /** - * \brief Current iteration number (and index to the current factor to update). - */ -// unsigned int ite; - /** \brief In calc_theta() two values are calculated for theta, this boolean is set to true to always choose theta2 (useful for GivensFGFTParallel). */ bool always_theta2; -// bool stoppingCritIsError; -// double stoppingError; -// bool errIsRel; - /** * Function pointer to any step of the algorithm (internal purpose only). */ @@ -128,15 +70,6 @@ namespace Faust { */ virtual void next_step(); - /** - * \brief Algo. main loop (facts.size() iterations). - * - * It's the entry point to user code. - * It calls next_step() in order to execute one iteration. - * - */ -// void compute_facts(); - /** * \brief Algo. step 2.4 * @@ -145,29 +78,12 @@ namespace Faust { virtual void update_L(); protected: - /** \brief Algo. step 2.1. - */ virtual void choose_pivot(); - /** \brief Algo. step 2.1.1 - * - * Computes the max of L or sorts it. - */ virtual void max_L(); - /** \brief Algo. step 2.2. - * - * Computes theta angle for the current Givens factor. - * - */ void calc_theta(); - /** - * \brief Algo. step 2.3. - * - * Updates the current Givens factor according to the pivot chosen for the iteration. - * - */ virtual void update_fact(); virtual void update_L(MatDense<FPP,Cpu> &); @@ -188,134 +104,14 @@ namespace Faust { void update_L_first(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); void update_L_second(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); - /** - * \brief Algo. step 2.5. - * - * Updates the diagonal approximate D from L (the diagonal part is taken). - * - */ -// void update_D(); - /** - * \brief Algo. step 2.6. - * - * Computes the error of approximation for the current iteration. - * - */ void update_err(); - /** - * Order (sort) D ascendantly into ordered_D and keeps ordered indices in ord_indices. - */ -// void order_D(); - - /** - * Order D into ascendantly or descendantly into order_D and keeps ordered indices in ord_indices. - * - * \param order -1 to order descendantly, 1 to order ascendantly. - */ -// void order_D(int order); - - public: - /** - * Returns the ordered indices of D to get increasing eigenvalues along the diagonal. - * - * @see order_D() - */ -// const vector<int>& get_ord_indices(); - - - /** - * Returns the vector of errors computed by calc_err() during the algorithm iterations. - */ -// const vector<FPP2>& get_errs() const; - /** - * Returns the specific j-th iteration's error (computed in calc_err()). - */ -// FPP2 get_err(int j) const; - - /** - * Returns the diag. vector D in its current status (which is updated at each iteration). - */ -// const Faust::Vect<FPP,DEVICE>& get_D(const bool ord=false); - - /** - * \param ord 0 (no sort), 1 (ascendant sort), -1 (descendant sort). - */ -// const Faust::Vect<FPP,DEVICE>& get_D(const int ord=0); - - - /** \brief Returns the diagonal vector as a sparse matrix. - * - * \note This function makes copies, is not intented for repeated use (use get_D() for optimized calls). - * - **/ const Faust::MatSparse<FPP,DEVICE> get_Dspm(const bool ord=false); - /** - * Returns the diagonal by copying it in the buffer diag_data (should be allocated by the callee). - * - * \param ord true to get the data ordered by ascendant eigenvalues false otherwise. - */ -// void get_D(FPP* diag_data, const bool ord=false); - - /** - * - * \param ord: 0 for no sort, -1 for descending order, 1 for descending order. - */ -// void get_D(FPP* diag_data, const int ord=0); - - - /** - * Computes and returns the Fourier matrix approximate from the Givens factors computed up to this time. - * - */ -// const MatDense<FPP,DEVICE> compute_fourier(const bool ord=false); - - /** - * Returns the matrix L. - * - * @note this matrix is not the Laplacian matrix from which the algorithm has started from. This is a matrix initialized to the Laplacian before the first iteration and then is computed such that L = S^t L. S being the current Givens factor matrix updated (i.e. facts[ite]). - * - * @see get_Lap() - * - */ - const MatGeneric<FPP,DEVICE>& get_L() const ; - - /** - * Returns the vector of all pivot choices (p, q) coordinates made by the algorithm until the last iteration. - */ -// const vector<pair<int,int>>& get_coord_choices() const; - - /** - * Returns the j-th iteration's pivot choice (p, q). - */ -// void get_coord_choice(int j, int& p, int& q) const; - - /** - * Returns the Laplacian matrix unmodified (as it was when passed to the constructor). - */ -// const MatDense<FPP,DEVICE>& get_Lap() const; - - /** - * Returns the vector of Givens matrices at this stage of algorithm execution (terminated or not). - */ -// const vector<MatSparse<FPP,DEVICE>>& get_facts() const; - - - /** - * Returns a Faust::Transform object with copy of facts into it. - * - * \param ord true to get the Transform's facts ordering the last one columns according to ascendant eigenvalues, false to let facts as they went out from the algorithm (without reordering). - */ -// Faust::Transform<FPP,DEVICE> get_transform(bool ord); - /** - * \param ord -1 for descending order, 1 for ascending order. - */ -// Faust::Transform<FPP,DEVICE> get_transform(int ord); }; diff --git a/src/algorithm/factorization/faust_GivensFGFT.hpp b/src/algorithm/factorization/faust_GivensFGFT.hpp index b5cf23f6cd55117ba5080dfd519e67cc0ca58a6e..ec920ad6609413c3cb2553d3204f2e5f8f0568a0 100644 --- a/src/algorithm/factorization/faust_GivensFGFT.hpp +++ b/src/algorithm/factorization/faust_GivensFGFT.hpp @@ -1,12 +1,5 @@ using namespace Faust; -// handle Visual Studio's whim https://msdn.microsoft.com/en-us/library/4hwaceh6.aspx -#define _USE_MATH_DEFINES -#include <cmath> -#ifndef M_PI_4 -#define M_PI_4 (M_PI/4.0) -#endif - template<typename FPP, Device DEVICE, typename FPP2> void GivensFGFT<FPP,DEVICE,FPP2>::next_step() { @@ -348,18 +341,6 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_L_first(Eigen::SparseMatrix<FPP, RowMaj L.mat.innerVector(q) = tmp; } -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFT<FPP,DEVICE,FPP2>::update_D() -//{ -// // D = spdiag(diag(L)) -// for(int i=0;i<L->getNbRow();i++) -// D.getData()[i] = (*L)(i,i); -//#ifdef DEBUG_GIVENS -// D.Display(); -// cout << "D fro. norm: " << D.norm() << endl; -//#endif -//} - template<typename FPP, Device DEVICE, typename FPP2> void GivensFGFT<FPP,DEVICE,FPP2>::update_err() { @@ -397,299 +378,27 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_err() } } -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFT<FPP,DEVICE,FPP2>::order_D() -//{ -// order_D(1); -//} -// -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFT<FPP,DEVICE,FPP2>::order_D(const int order /* -1 for descending order, 1 for ascending order */) -//{ -// ordered_D = Faust::Vect<FPP,DEVICE>(D.size()); -// ord_indices.resize(0); -// for(int i=0;i<D.size();i++) -// ord_indices.push_back(i); -// sort(ord_indices.begin(), ord_indices.end(), [this, &order](int i, int j) { -// return order>0?D.getData()[i] < D.getData()[j]:(order <0?D.getData()[i] > D.getData()[j]:0); -// }); -// for(int i=0;i<ord_indices.size();i++) -// ordered_D.getData()[i] = D.getData()[ord_indices[i]]; -// // compute inverse permutation to keep easily possible retrieving undefined order -// inv_ord_indices.resize(ord_indices.size()); -// int j = 0, i = 0; -// while(j<ord_indices.size()) -// if(ord_indices[i] == j) -// { -// inv_ord_indices[j++] = i; -// i = 0; -// } -// else -// i++; -// is_D_ordered = true; -// D_order_dir = order; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<int>& GivensFGFT<FPP,DEVICE,FPP2>::get_ord_indices() -//{ -// if(! is_D_ordered) -// order_D(); -// return ord_indices; -//} - - -//template<typename FPP, Device DEVICE, typename FPP2> -//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(stopping = (ite > 1 && stoppingCritIsError && 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; -// } -// } -//} -// template<typename FPP, Device DEVICE, typename FPP2> -GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), /* Lap(Lap), facts(J>0?J:1), D(Lap.getNbRow()), */ C(Lap.getNbRow(), Lap.getNbCol()), /* errs(0), coord_choices(0), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false),*/ always_theta2(false) /*,verbosity(verbosity), stoppingCritIsError(stoppingError != 0.0), stoppingError(stoppingError), errIsRel(errIsRel), last_fact_permuted(false), Lap_squared_fro_norm(0), J(J)*/ - +GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol()), always_theta2(false) { - /* Matlab ref. code: - * facts = cell(1,J); - * n=size(Lap,1); - * L=Lap; - * C = 15*ones(n); - * err=zeros(1,J); - * coord_choices = zeros(2,J); - * - */ -// if(J == 0 && ! stoppingCritIsError) handleError("GivensFGFT", "Either J or stoppingError must be > 0"); + //see parent ctor C.setOnes(); C.scalarMultiply(15); // purely abitrary -// if(Lap.getNbCol() != Lap.getNbRow()) -// handleError("Faust::GivensFGFT", "Laplacian must be a square matrix."); -// -// // init the identity part of the factor buffer model -// // allocate the mem. space for the 4 additional rotation part coeffs -// for(int i=0;i<Lap.getNbRow();i++) -// { -// fact_mod_values.push_back(FPP(1)); -// fact_mod_col_ids.push_back(i); -// fact_mod_row_ids.push_back(i); -// } -// -// // init. D -// memset(D.getData(), 0, sizeof(FPP)*Lap.getNbRow()); -// -// L = new MatSparse<FPP,DEVICE>(Lap); } template<typename FPP, Device DEVICE, typename FPP2> -GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), /* Lap(Lap), facts(J>0?J:1), D(Lap.getNbRow()), */ C(Lap.getNbRow(), Lap.getNbCol()),/* errs(0), coord_choices(0), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false),*/ always_theta2(false) /*, verbosity(verbosity), stoppingCritIsError(stoppingError != 0.0), stoppingError(stoppingError), errIsRel(errIsRel), last_fact_permuted(false), Lap_squared_fro_norm(0), J(J) */ +GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol()), always_theta2(false) { - /* Matlab ref. code: - * facts = cell(1,J); - * n=size(Lap,1); - * L=Lap; - * C = 15*ones(n); - * err=zeros(1,J); - * coord_choices = zeros(2,J); - * - */ -// if(J == 0 && ! stoppingCritIsError) handleError("GivensFGFT", "Either J or stoppingError must be > 0"); + // see parent ctor C.setOnes(); C.scalarMultiply(15); // purely abitrary -// if(Lap.getNbCol() != Lap.getNbRow()) -// handleError("Faust::GivensFGFT", "Laplacian must be a square matrix."); -// -// // init the identity part of the factor buffer model -// // allocate the mem. space for the 4 additional rotation part coeffs -// for(int i=0;i<Lap.getNbRow();i++) -// { -// fact_mod_values.push_back(FPP(1)); -// fact_mod_col_ids.push_back(i); -// fact_mod_row_ids.push_back(i); -// } -// -// // init. D -// memset(D.getData(), 0, sizeof(FPP)*Lap.getNbRow()); -// -// L = new MatDense<FPP,DEVICE>(Lap); } -//template<typename FPP, Device DEVICE, typename FPP2> -//FPP2 GivensFGFT<FPP,DEVICE,FPP2>::get_err(int j) const -//{ -// if(j > 0 && j < errs.size()) -// return errs[j]; -// else -// throw out_of_range("GivensFGFT::get_err(j): j is out of range."); -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<FPP2>& GivensFGFT<FPP,DEVICE,FPP2>::get_errs() const -//{ -// return errs; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::Vect<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_D(const bool ord /* default to false */) -//{ -// if(ord) -// { -// if(!is_D_ordered) -// order_D(); -// return ordered_D; -// } -// return D; -//} -// -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::Vect<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_D(const int ord /* default to false */) -//{ -// if(ord != 0) -// { -// if(!is_D_ordered || ord != D_order_dir) -// order_D(ord); -// return ordered_D; -// } -// return D; -//} - template<typename FPP, Device DEVICE, typename FPP2> const Faust::MatSparse<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::get_Dspm(const bool ord /* default to false */) { MatSparse <FPP,DEVICE> spD; GivensFGFTGen<FPP,DEVICE,FPP2>::get_Dspm(spD, ord); return spD; -// const Faust::Vect<FPP,DEVICE>& D_ = this->get_D(ord); -// vector<int> nat_ord_indices; -// vector<FPP> diag_D; -// for(int i=0;i<D_.size();i++) -// { -// nat_ord_indices.push_back(i); -// diag_D.push_back(D_.getData()[i]); -// } -// MatSparse<FPP,DEVICE> spD(nat_ord_indices, nat_ord_indices, diag_D, nat_ord_indices.size(), nat_ord_indices.size()); -// return spD; } -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFT<FPP,DEVICE,FPP2>::get_D(FPP* diag_data, const bool ord /* default to false */) -//{ -// get_D(diag_data, ord?1:0); -//} -// -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFT<FPP,DEVICE,FPP2>::get_D(FPP* diag_data, const int ord /* default to false */) -//{ -// const Faust::Vect<FPP,DEVICE>& D_ = get_D(ord); -// const FPP* src_data_ptr = D_.getData(); -// memcpy(diag_data, src_data_ptr, sizeof(FPP)*D_.size()); -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::MatDense<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::compute_fourier(const bool ord /* default to false */) -//{ -// Faust::MatDense<FPP,Cpu> fourier(Lap.getNbRow(), Lap.getNbCol()); -// Faust::MatDense<FPP,Cpu>* ord_fourier; -// fourier.setEyes(); -// for(int i=facts.size()-1;i>=0;i--) -// facts[i].multiply(fourier, 'N'); -// if(ord) -// { -// if(!is_D_ordered) -// order_D(); -// ord_fourier = fourier.get_cols(ord_indices); -// fourier = *ord_fourier; -// delete ord_fourier; -// } -// return fourier; -//} - - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::MatGeneric<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_L() const -//{ -// return *L; -//} -// -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<pair<int,int>>& GivensFGFT<FPP,DEVICE,FPP2>::get_coord_choices() const -//{ -// return coord_choices; -//} -// -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFT<FPP,DEVICE,FPP2>::get_coord_choice(int j, int& p, int& q) const -//{ -// if(j > 0 && j < coord_choices.size()) -// { -// p = coord_choices[j].first; -// q = coord_choices[j].second; -// } -// else -// throw out_of_range("GivensFGFT::get_coord_choice(j,p,q): j is out of range."); -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::MatDense<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_Lap() const -//{ -// return Lap; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<Faust::MatSparse<FPP,DEVICE>>& GivensFGFT<FPP,DEVICE,FPP2>::get_facts() const -//{ -// return facts; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//Faust::Transform<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::get_transform(bool ord) -//{ -// return get_transform(ord?1:0); -//} -// -// -//template<typename FPP, Device DEVICE, typename FPP2> -//Faust::Transform<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::get_transform(int ord) -//{ -// //TODO: an optimization is possible by changing type of facts to vector<MatGeneric*> it would avoid copying facts into Transform and rather use them directly. It will need a destructor that deletes them eventually if they weren't transfered to a Transform object before. -// MatSparse<FPP,DEVICE> & last_fact = *(facts.end()-1); -// MatSparse<FPP,DEVICE> P(last_fact.getNbCol(), last_fact.getNbCol()); //last_fact permutation matrix -// // (to reorder eigenvector transform according to ordered D) -// -// if(last_fact_permuted) -// { -// // get back to undefined order -// for(int i=0;i<inv_ord_indices.size();i++) -// P.setCoeff(inv_ord_indices[i],i, FPP(1.0)); -// last_fact.multiplyRight(P); -// last_fact_permuted = false; -// } -// -// if(ord) -// { -// if(!is_D_ordered || ord != D_order_dir) -// order_D(ord); -// for(int i=0;i<ord_indices.size();i++) -// P.setCoeff(ord_indices[i],i, FPP(1.0)); -// // P.set_orthogonal(true); -// // facts.push_back(P); // we prefer to directly multiply the last factor by P -// last_fact_permuted = true; -// last_fact.multiplyRight(P); -// } -// Faust::Transform<FPP,DEVICE> t = Faust::Transform<FPP,DEVICE>(facts); -// // // remove the permutation factor if added temporarily for reordering -// // return ord?facts.erase(facts.end()-1),t:t; -// return t; -//} diff --git a/src/algorithm/factorization/faust_GivensFGFTComplex.h b/src/algorithm/factorization/faust_GivensFGFTComplex.h index 13698030962264b31315e9805aa8317ba4d5e14e..3c8078a74397f1341048662e18ffbec2d4e6f17c 100644 --- a/src/algorithm/factorization/faust_GivensFGFTComplex.h +++ b/src/algorithm/factorization/faust_GivensFGFTComplex.h @@ -20,6 +20,8 @@ namespace Faust { * \brief This class implements the Givens FGFT algorithm (Truncated Jacobi algorithm) for the complex matrix case (ideal case being the Hermitian matrix case). * This algorithm is based on the classical Jacobi eigenvalues algorithm. * + * See parent class GivensFGFTGen for documentation about members. + * * References: * * [1] Le Magoarou L., Gribonval R. and Tremblay N., "Approximate fast @@ -36,71 +38,11 @@ namespace Faust { Faust::MatDense<FPP,DEVICE> C; /** \brief Column vector for the rowwise minimization of C (i.e. maximization of L). */ Faust::Vect<FPP,DEVICE> C_min_row; - /** \brief Pivot candidates q coordinates. */ - //int* q_candidates; /* default IndexType for underlying eigen matrix is int. */ public: using DType = typename FPP::value_type; protected: - /** \brief The number of targeted transform factors (when only one Givens rotation is stored per factor).*/ -// int J; - /** \brief Fourier matrix/eigenvectors factorization matrices (Givens matrix). */ -// vector<Faust::MatSparse<FPP,DEVICE>> facts; - /** \brief Diagonalization approximate of Laplacian. */ - //Faust::Vect<DType,DEVICE> D; - /** \brief Queue of errors (cf. calc_err()). */ -// vector<FPP2> errs; - /** \brief Pivot choices (p, q) coordinates. */ -// vector<pair<int,int>> coord_choices; - /** \brief Graph Laplacian to diagonalize/approximate. */ -// Faust::MatGeneric<FPP, DEVICE>& Lap; -// FPP2 Lap_squared_fro_norm; - /** \brief L iteration factor: L_i = S^T L_{i-1} S, initialized from Lap (with S being facts[i]). */ -// Faust::MatGeneric<FPP, DEVICE> *L; - /** \brief Rotation angle theta for the current iteration's Givens matrix. */ + /** \brief Rotation angle for the current iteration's "Givens" matrix. */ FPP theta1, theta2; -// bool last_fact_permuted; - /* Precomputed model identity matrix to init. facts[ite] before update. - * Identity matrix is completed later with cos/sin coefficients (in update_fact()). - */ - - /** \brief Defines the rows of facts[ite]. */ -// vector<int> fact_mod_row_ids; -// /** \brief Defines the columns of facts[ite]. */ -// vector<int> fact_mod_col_ids; -// /** \brief Defines the coefficients of facts[ite]. */ -// vector<FPP> fact_mod_values; - - - /** \brief Ordered indices of D to get increasing eigenvalues along the diagonal. */ -// vector<int> ord_indices; - /** \brief inverse permutation of ord_indices (needed to retrieve start undefined order). */ -// vector<int> inv_ord_indices; - - /** \brief Cache for the ordered D. */ -// Faust::Vect<DType,DEVICE> ordered_D; - /** \brief true if D has already been ordered (order_D() was called). */ - //bool is_D_ordered; - //int D_order_dir; - - /** \brief The level of verbosity (0 for nothing, 1 for iteration numbers,...) */ - //unsigned int verbosity; - - /** - * \brief Row index for the selected pivot in L. - */ -// int p, -// /** -// * \brief Column index for the selected pivot in L. -// */q; -// - /** - * \brief Current iteration number (and index to the current factor to update). - */ - // unsigned int ite; - -// bool stoppingCritIsError; -// double stoppingError; -// bool errIsRel; /** * Function pointer to any step of the algorithm (internal purpose only). @@ -111,56 +53,18 @@ namespace Faust { public: const static unsigned int ERROR_CALC_PERIOD = 100; - /** Algorithm class constructor. - * \param Lap The Laplacian matrix to approximate/diagonalize. - * \param J The number of iterations, Givens rotations factors. - * */ GivensFGFTComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true); GivensFGFTComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true); - /** Destructor */ -// virtual ~GivensFGFTComplex() {delete L;}; - /** - * \brief Algo. main step. - * - * This function calls all step functions of the algorithm in the proper order to execute only one iteration. - * - * External code may call this function instead of compute_facts() in order to debug iterations taken one by one. - */ virtual void next_step(); - /** - * \brief Algo. main loop (facts.size() iterations). - * - * It's the entry point to user code. - * It calls next_step() in order to execute one iteration. - * - */ -// void compute_facts(); - - /** - * \brief Algo. step 2.4 - * - * Updates L after Givens factor update for the next iteration. - */ virtual void update_L(); protected: - /** \brief Algo. step 2.1. - */ virtual void choose_pivot(); - /** \brief Algo. step 2.1.1 - * - * Computes the max of L or sorts it. - */ virtual void max_L(); - /** \brief Algo. step 2.2. - * - * Computes theta angle for the current Givens factor. - * - */ void calc_theta(); /** @@ -168,12 +72,6 @@ namespace Faust { */ void check_pivot_image(FPP& c_pp, FPP& c_pq, FPP& c_qp, FPP& c_qq); - /** - * \brief Algo. step 2.3. - * - * Updates the current Givens factor according to the pivot chosen for the iteration. - * - */ virtual void update_fact(); virtual void update_L(MatDense<FPP,Cpu> &); @@ -194,136 +92,13 @@ namespace Faust { void update_L_first(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP& c_pp, const FPP& c_pq, const FPP& c_qp, const FPP& c_qq, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); void update_L_second(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP& c_pp, const FPP& c_pq, const FPP& c_qp, const FPP& c_qq, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); - /** - * \brief Algo. step 2.5. - * - * Updates the diagonal approximate D from L (the diagonal part is taken). - * - */ -// void update_D(); - /** - * \brief Algo. step 2.6. - * - * Computes the error of approximation for the current iteration. - * - */ void update_err(); - /** - * Order (sort) D ascendantly into ordered_D and keeps ordered indices in ord_indices. - */ -// void order_D(); - - /** - * Order D into ascendantly or descendantly into order_D and keeps ordered indices in ord_indices. - * - * \param order -1 to order descendantly, 1 to order ascendantly. - */ -// void order_D(int order); - - public: - /** - * Returns the ordered indices of D to get increasing eigenvalues along the diagonal. - * - * @see order_D() - */ -// const vector<int>& get_ord_indices(); - - - /** - * Returns the vector of errors computed by calc_err() during the algorithm iterations. - */ -// const vector<FPP2>& get_errs() const; - - /** - * Returns the specific j-th iteration's error (computed in calc_err()). - */ -// FPP2 get_err(int j) const; - - /** - * Returns the diag. vector D in its current status (which is updated at each iteration). - */ -// const Faust::Vect<DType /*FPP*/,DEVICE>& get_D(const bool ord=false); - - /** - * \param ord 0 (no sort), 1 (ascendant sort), -1 (descendant sort). - */ -// const Faust::Vect<DType /*FPP*/,DEVICE>& get_D(const int ord=0); - - - /** \brief Returns the diagonal vector as a sparse matrix. - * - * \note This function makes copies, is not intented for repeated use (use get_D() for optimized calls). - * - **/ const Faust::MatSparse<FPP,DEVICE> get_Dspm(const bool ord=false); - /** - * Returns the diagonal by copying it in the buffer diag_data (should be allocated by the callee). - * - * \param ord true to get the data ordered by ascendant eigenvalues false otherwise. - */ -// void get_D(DType/*FPP*/* diag_data, const bool ord=false); - - /** - * - * \param ord: 0 for no sort, -1 for descending order, 1 for descending order. - */ -// void get_D(DType/*FPP*/* diag_data, const int ord=0); - - - /** - * Computes and returns the Fourier matrix approximate from the Givens factors computed up to this time. - * - */ -// const MatDense<FPP,DEVICE> compute_fourier(const bool ord=false); - - /** - * Returns the matrix L. - * - * @note this matrix is not the Laplacian matrix from which the algorithm has started from. This is a matrix initialized to the Laplacian before the first iteration and then is computed such that L = S^t L. S being the current Givens factor matrix updated (i.e. facts[ite]). - * - * @see get_Lap() - * - */ -// const MatGeneric<FPP,DEVICE>& get_L() const ; - - /** - * Returns the vector of all pivot choices (p, q) coordinates made by the algorithm until the last iteration. - */ -// const vector<pair<int,int>>& get_coord_choices() const; - - /** - * Returns the j-th iteration's pivot choice (p, q). - */ -// void get_coord_choice(int j, int& p, int& q) const; - - /** - * Returns the Laplacian matrix unmodified (as it was when passed to the constructor). - */ -// const MatDense<FPP,DEVICE>& get_Lap() const; - - /** - * Returns the vector of Givens matrices at this stage of algorithm execution (terminated or not). - */ -// const vector<MatSparse<FPP,DEVICE>>& get_facts() const; - - - /** - * Returns a Faust::Transform object with copy of facts into it. - * - * \param ord true to get the Transform's facts ordering the last one columns according to ascendant eigenvalues, false to let facts as they went out from the algorithm (without reordering). - */ -// Faust::Transform<FPP,DEVICE> get_transform(bool ord); - /** - * \param ord -1 for descending order, 1 for ascending order. - */ -// Faust::Transform<FPP,DEVICE> get_transform(int ord); - - }; } diff --git a/src/algorithm/factorization/faust_GivensFGFTComplex.hpp b/src/algorithm/factorization/faust_GivensFGFTComplex.hpp index 796ea66b208b9985adcb6fc36d36e54b67c5a062..b44cb3719df12a36b4f384a3fb101e3bb22656f1 100644 --- a/src/algorithm/factorization/faust_GivensFGFTComplex.hpp +++ b/src/algorithm/factorization/faust_GivensFGFTComplex.hpp @@ -1,12 +1,5 @@ using namespace Faust; -// handle Visual Studio's whim https://msdn.microsoft.com/en-us/library/4hwaceh6.aspx -#define _USE_MATH_DEFINES -#include <cmath> -#ifndef M_PI_4 -#define M_PI_4 (M_PI/4.0) -#endif - template<typename FPP, Device DEVICE, typename FPP2> void GivensFGFTComplex<FPP,DEVICE,FPP2>::next_step() { @@ -414,12 +407,6 @@ void GivensFGFTComplex<FPP,DEVICE,FPP2>::update_L_first(Eigen::SparseMatrix<FPP, L.mat.innerVector(q) = tmp; } -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFTComplex<FPP,DEVICE,FPP2>::update_D() -//{ -// GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>::update_D(this->L); -//} - template<typename FPP, Device DEVICE, typename FPP2> void GivensFGFTComplex<FPP,DEVICE,FPP2>::update_err() { @@ -457,303 +444,26 @@ void GivensFGFTComplex<FPP,DEVICE,FPP2>::update_err() } } -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFTComplex<FPP,DEVICE,FPP2>::order_D() -//{ -// order_D(1); -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFTComplex<FPP,DEVICE,FPP2>::order_D(const int order /* -1 for descending order, 1 for ascending order */) -//{ -// ordered_D = Faust::Vect<DType /*FPP*/,DEVICE>(this->D.size()); -// ord_indices.resize(0); -// for(int i=0;i<this->D.size();i++) -// ord_indices.push_back(i); -// sort(ord_indices.begin(), ord_indices.end(), [this, &order](int i, int j) { -// return order>0?this->D.getData()[i]/*.real()*/ < this->D.getData()[j]/*.real()*/:(order <0?this->D.getData()[i]/*.real()*/ > this->D.getData()[j]/*.real()*/:0); -// }); -// for(int i=0;i<ord_indices.size();i++) -// { -// ordered_D.getData()[i] = this->D.getData()[ord_indices[i]]; -// } -// // compute inverse permutation to keep easily possible retrieving undefined order -// inv_ord_indices.resize(ord_indices.size()); -// int j = 0, i = 0; -// while(j<ord_indices.size()) -// if(ord_indices[i] == j) -// { -// inv_ord_indices[j++] = i; -// i = 0; -// } -// else -// i++; -// is_D_ordered = true; -// D_order_dir = order; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<int>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_ord_indices() -//{ -// if(! is_D_ordered) -// order_D(); -// return ord_indices; -//} - - -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFTComplex<FPP,DEVICE,FPP2>::compute_facts() -//{ -// this->is_D_ordered = false; // facts (re)computed then D must be reordered -// this->ite = 0; -// bool stopping = false; -// while(this->J == 0 || this->ite < this->facts.size()) // when J == 0 the stopping criterion is the error against Lap -// { -// next_step(); -// this->ite++; -// if(stopping = (this->ite > 1 && this->stoppingCritIsError && this->errs.size() > 2 && this->errs[this->ite-1]-this->errs[this->ite-2] > FLT_EPSILON)) -// if(this->verbosity>0) cerr << "warning: the algorithm stopped because the last error is greater than the previous one." << endl; -// if(stopping || this->stoppingCritIsError && this->errs.size() > 0 && (*(this->errs.end()-1) - this->stoppingError ) < FLT_EPSILON) -// { -// this->facts.resize(this->ite); -// break; -// } -// } -//} - template<typename FPP, Device DEVICE, typename FPP2> -GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel)/*, Lap(Lap), facts(J>0?J:1), D(Lap.getNbRow())*/, C(Lap.getNbRow(), Lap.getNbCol())/*,errs(0), coord_choices(0), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false), verbosity(verbosity), stoppingCritIsError(stoppingError != 0.0), stoppingError(stoppingError), errIsRel(errIsRel),last_fact_permuted(false), Lap_squared_fro_norm(FPP2(0)) J(J)*/ +GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol()) { - /* Matlab ref. code: - * facts = cell(1,J); - * n=size(Lap,1); - * L=Lap; - * C = 15*ones(n); - * err=zeros(1,J); - * coord_choices = zeros(2,J); - * - */ -// if(this->J == 0 && ! this->stoppingCritIsError) handleError("GivensFGFT", "Either J or stoppingError must be > 0"); + //see parent ctor this->C.setZeros(); - // C.setOnes(); - // C.scalarMultiply(-15); // purely abitrary -// if(Lap.getNbCol() != Lap.getNbRow()) -// handleError("Faust::GivensFGFTComplex", "Laplacian must be a square matrix."); - -// // init the identity part of the factor buffer model -// // allocate the mem. space for the 4 additional rotation part coeffs -// for(int i=0;i<Lap.getNbRow();i++) -// { -// fact_mod_values.push_back(FPP(1)); -// fact_mod_col_ids.push_back(i); -// fact_mod_row_ids.push_back(i); -// } - - // init. D -// memset(this->D.getData(), 0, sizeof(DType/*FPP*/)*Lap.getNbRow()); - -// L = new MatSparse<FPP,DEVICE>(Lap); } template<typename FPP, Device DEVICE, typename FPP2> -GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel) /* Lap(Lap), facts(J>0?J:1), D(Lap.getNbRow())*/, C(Lap.getNbRow(), Lap.getNbCol())/*, errs(0), coord_choices(0), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false), verbosity(verbosity), stoppingCritIsError(stoppingError != 0.0), stoppingError(stoppingError), errIsRel(errIsRel), last_fact_permuted(false), Lap_squared_fro_norm(FPP2(0)) J(J)*/ +GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol()) { - /* Matlab ref. code: - * facts = cell(1,J); - * n=size(Lap,1); - * L=Lap; - * C = 15*ones(n); - * err=zeros(1,J); - * coord_choices = zeros(2,J); - * - */ -// if(this->J == 0 && ! this->stoppingCritIsError) handleError("GivensFGFT", "Either J or stoppingError must be > 0"); + // see parent ctor this->C.setZeros(); -// C.setOnes(); -// C.scalarMultiply(FPP(-15)); // purely abitrary -// if(Lap.getNbCol() != Lap.getNbRow()) -// handleError("Faust::GivensFGFTComplex", "Laplacian must be a square matrix."); - -// // init the identity part of the factor buffer model -// // allocate the mem. space for the 4 additional rotation part coeffs -// for(int i=0;i<Lap.getNbRow();i++) -// { -// fact_mod_values.push_back(FPP(1)); -// fact_mod_col_ids.push_back(i); -// fact_mod_row_ids.push_back(i); -// } - - // init. D - //memset(this->D.getData(), 0, sizeof(DType/*FPP*/)*Lap.getNbRow()); - -// L = new MatDense<FPP,DEVICE>(Lap); } -//template<typename FPP, Device DEVICE, typename FPP2> -//FPP2 GivensFGFTComplex<FPP,DEVICE,FPP2>::get_err(int j) const -//{ -// if(j > 0 && j < errs.size()) -// return errs[j]; -// else -// throw out_of_range("GivensFGFTComplex::get_err(j): j is out of range."); -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<FPP2>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_errs() const -//{ -// return errs; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::Vect<typename FPP::value_type/*FPP*/,DEVICE>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_D(const bool ord /* default to false */) -//{ -// if(ord) -// { -// if(!is_D_ordered) -// order_D(); -// return ordered_D; -// } -// return D; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::Vect<typename FPP::value_type/*FPP*/,DEVICE>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_D(const int ord /* default to false */) -//{ -// if(ord != 0) -// { -// if(!is_D_ordered || ord != D_order_dir) -// order_D(ord); -// return ordered_D; -// } -// return D; -//} - template<typename FPP, Device DEVICE, typename FPP2> -const Faust::MatSparse<FPP,DEVICE> GivensFGFTComplex<FPP,DEVICE,FPP2>::get_Dspm(const bool ord /* default to false*/) +const Faust::MatSparse<FPP,DEVICE> GivensFGFTComplex<FPP,DEVICE,FPP2>::get_Dspm(const bool ord /* default to false*/) { MatSparse <FPP,DEVICE> spD; GivensFGFTGen<typename FPP::value_type,DEVICE,FPP2,FPP>::get_Dspm(spD, ord); return spD; -// const Faust::Vect<FPP,DEVICE>& D_ = this->get_D(ord); -// vector<int> nat_ord_indices; -// vector<FPP> diag_D; -// for(int i=0;i<D_.size();i++) -// { -// nat_ord_indices.push_back(i); -// diag_D.push_back(FPP(D_.getData()[i])); -// } -// MatSparse<FPP,DEVICE> spD(nat_ord_indices, nat_ord_indices, diag_D, nat_ord_indices.size(), nat_ord_indices.size()); -// return spD; } -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFTComplex<FPP,DEVICE,FPP2>::get_D(DType* /*FPP**/ diag_data, const bool ord /* default to false */) -//{ -// get_D(diag_data, ord?1:0); -//} -// -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFTComplex<FPP,DEVICE,FPP2>::get_D(DType* /*FPP**/ diag_data, const int ord /* default to false */) -//{ -// const Faust::Vect<FPP,DEVICE>& D_ = get_D(ord); -// const DType* /*FPP**/ src_data_ptr = D_.getData(); -// memcpy(diag_data, src_data_ptr, sizeof(FPP)*D_.size()); -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::MatDense<FPP,DEVICE> GivensFGFTComplex<FPP,DEVICE,FPP2>::compute_fourier(const bool ord /* default to false */) -//{ -// Faust::MatDense<FPP,Cpu> fourier(L->getNbRow(), L->getNbCol()); -// Faust::MatDense<FPP,Cpu>* ord_fourier; -// fourier.setEyes(); -// for(int i=facts.size()-1;i>=0;i--) -// facts[i].multiply(fourier, 'N'); -// if(ord) -// { -// if(!this->is_D_ordered) -// this->order_D(); -// ord_fourier = fourier.get_cols(this->ord_indices); -// fourier = *ord_fourier; -// delete ord_fourier; -// } -// return fourier; -//} - - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::MatGeneric<FPP,DEVICE>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_L() const -//{ -// return *L; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<pair<int,int>>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_coord_choices() const -//{ -// return coord_choices; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//void GivensFGFTComplex<FPP,DEVICE,FPP2>::get_coord_choice(int j, int& p, int& q) const -//{ -// if(j > 0 && j < coord_choices.size()) -// { -// p = coord_choices[j].first; -// q = coord_choices[j].second; -// } -// else -// throw out_of_range("GivensFGFTComplex::get_coord_choice(j,p,q): j is out of range."); -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const Faust::MatDense<FPP,DEVICE>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_Lap() const -//{ -// return Lap; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//const vector<Faust::MatSparse<FPP,DEVICE>>& GivensFGFTComplex<FPP,DEVICE,FPP2>::get_facts() const -//{ -// return this->facts; -//} - -//template<typename FPP, Device DEVICE, typename FPP2> -//Faust::Transform<FPP,DEVICE> GivensFGFTComplex<FPP,DEVICE,FPP2>::get_transform(bool ord) -//{ -// return get_transform(ord?1:0); -//} -// -// -//template<typename FPP, Device DEVICE, typename FPP2> -//Faust::Transform<FPP,DEVICE> GivensFGFTComplex<FPP,DEVICE,FPP2>::get_transform(int ord) -//{ -// //TODO: an optimization is possible by changing type of facts to vector<MatGeneric*> it would avoid copying facts into Transform and rather use them directly. It will need a destructor that deletes them eventually if they weren't transfered to a Transform object before. -// MatSparse<FPP,DEVICE> & last_fact = *(this->facts.end()-1); -// MatSparse<FPP,DEVICE> P(last_fact.getNbCol(), last_fact.getNbCol()); //last_fact permutation matrix -// // (to reorder eigenvector transform according to ordered D) -// -// if(last_fact_permuted) -// { -// // non-ordered eigenvector transform (ord == 0) or opposite order asked (ord != D_order_dir) -// // get back to undefined order -// for(int i=0;i<this->inv_ord_indices.size();i++) -// P.setCoeff(this->inv_ord_indices[i],i, FPP(1.0)); -// last_fact.multiplyRight(P); -// last_fact_permuted = false; -// } -// -// if(ord) -// { -// if(!this->is_D_ordered || ord != this->D_order_dir) -// this->order_D(ord); -// for(int i=0;i<this->ord_indices.size();i++) -// P.setCoeff(this->ord_indices[i],i, FPP(1.0)); -// // P.set_orthogonal(true); -// // facts.push_back(P); // we prefer to directly multiply the last factor by P -// last_fact_permuted = true; -// last_fact.multiplyRight(P); -// } -// Faust::Transform<FPP,DEVICE> t = Faust::Transform<FPP,DEVICE>(this->facts); -// // // remove the permutation factor if added temporarily for reordering -// // return ord?facts.erase(facts.end()-1),t:t; -// return t; -//} + diff --git a/src/algorithm/factorization/faust_GivensFGFTGen.h b/src/algorithm/factorization/faust_GivensFGFTGen.h index b6ca12ae0adcd9efc9c4360b83506ba118a96689..7aa756da7fcf92cb0fc7b6b09cd88e3d81052a8c 100644 --- a/src/algorithm/factorization/faust_GivensFGFTGen.h +++ b/src/algorithm/factorization/faust_GivensFGFTGen.h @@ -30,11 +30,11 @@ namespace Faust { * */ /** \brief Temporary storage matrix for maximization of L. */ -// Faust::MatDense<FPP4,DEVICE> C; + // Faust::MatDense<FPP4,DEVICE> C; /** \brief Column vector for the rowwise minimization of C (i.e. maximization of L). */ -// Faust::Vect<FPP4,DEVICE> C_min_row; + // Faust::Vect<FPP4,DEVICE> C_min_row; protected: - /** \brief Fourier matrix/eigenvectors factorization matrices (Givens matrix). */ + /** \brief Fourier matrix/eigenvectors factorization matrices (Givens matrices). */ vector<Faust::MatSparse<FPP4,DEVICE>> facts; /** \brief L iteration factor: L_i = S^T L_{i-1} S, initialized from Lap (with S being facts[i]). */ @@ -44,22 +44,20 @@ namespace Faust { /** \brief The number of targeted transform factors (when only one Givens rotation is stored per factor).*/ int J; - /** \brief Fourier matrix factorization matrices (Givens matrix). */ -// vector<Faust::MatSparse<FPP,DEVICE>> facts; - /** \brief Diagonalization approximate of Laplacian. */ + /** \brief approximate eigenvalues vector. */ Faust::Vect<FPP,DEVICE> D; - /** \brief Queue of errors (cf. calc_err()). */ + /** \brief Queue of errors (cf. update_err()). */ vector<FPP2> errs; /** \brief Pivot choices (p, q) coordinates. */ vector<pair<int,int>> coord_choices; /** \brief Graph Laplacian to diagonalize/approximate. */ Faust::MatGeneric<FPP4, DEVICE>& Lap; + /** \brief Laplacian dimension */ unsigned int dim_size; + /** \brief Laplacian Frobenius norm */ FPP2 Lap_squared_fro_norm; - /** \brief L iteration factor: L_i = S^T L_{i-1} S, initialized from Lap (with S being facts[i]). */ -// Faust::MatGeneric<FPP, DEVICE> *L; /** \brief Rotation angle theta for the current iteration's Givens matrix. */ -// FPP2 theta; + // FPP2 theta; /* Precomputed model identity matrix to init. facts[ite] before update. * Identity matrix is completed later with cos/sin coefficients (in update_fact()). @@ -84,6 +82,7 @@ namespace Faust { Faust::Vect<FPP,DEVICE> ordered_D; /** \brief true if D has already been ordered (order_D() was called). */ bool is_D_ordered; + /** \brief 1 if eigenvalues has been ordered in ascending order -1 ortherwise (other values is for undefined order). */ int D_order_dir; /** \brief The level of verbosity (0 for nothing, 1 for iteration numbers,...) */ @@ -102,33 +101,27 @@ namespace Faust { */ unsigned int ite; - /** \brief In calc_theta() two values are calculated for theta, this boolean is set to true to always choose theta2 (useful for GivensFGFTGenParallel). */ -// bool always_theta2; - + /** \brief true if the stopping criterion is error (otherwise it's the number of iterations/number of Givens matrices) */ bool stoppingCritIsError; + /** \brief error value according to algorithm stops if stoppingCritIsError == true */ double stoppingError; + /** \brief true if the stopping error is taken as relative error (absolute otherwise). */ bool errIsRel; - /** - * Function pointer to any step of the algorithm (internal purpose only). - */ -// typedef void (GivensFGFTGen<FPP,DEVICE,FPP2>::*substep_fun)(); - public: -// const static unsigned int ERROR_CALC_PERIOD = 100; + // const static unsigned int ERROR_CALC_PERIOD = 100; /** Algorithm class constructor. * \param Lap The Laplacian matrix to approximate/diagonalize. * \param J The number of iterations, Givens rotations factors. + * TODO: complete argument list * */ GivensFGFTGen(Faust::MatGeneric<FPP4,DEVICE>* Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel); GivensFGFTGen(Faust::MatSparse<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel); GivensFGFTGen(Faust::MatDense<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel); -// GivensFGFTGen(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true); -// GivensFGFTGen(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true); /** Destructor */ virtual ~GivensFGFTGen() {delete[] q_candidates; delete L;}; @@ -183,24 +176,24 @@ namespace Faust { */ virtual void update_fact()=0; -// virtual void update_L(MatDense<FPP,Cpu> &); -// -// virtual void update_L(MatSparse<FPP,Cpu> &); + // virtual void update_L(MatDense<FPP,Cpu> &); + // + // virtual void update_L(MatSparse<FPP,Cpu> &); /** * Computes the first S'*L (only used by update_L() in optimization enabled code). * */ -// void update_L_first(Faust::Vect<FPP,DEVICE>& L_vec_p, Faust::Vect<FPP,DEVICE>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, MatDense<FPP,DEVICE> & L); + // void update_L_first(Faust::Vect<FPP,DEVICE>& L_vec_p, Faust::Vect<FPP,DEVICE>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, MatDense<FPP,DEVICE> & L); /** * Computes L*S (only used by update_L() in optimization enabled code). * Must be called after update_L_first() to finish the update of L = S'*L*S. */ -// void update_L_second(Faust::Vect<FPP,DEVICE>& L_vec_p, Faust::Vect<FPP,DEVICE>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, MatDense<FPP,DEVICE> & L); + // void update_L_second(Faust::Vect<FPP,DEVICE>& L_vec_p, Faust::Vect<FPP,DEVICE>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, MatDense<FPP,DEVICE> & L); -// void update_L_first(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); -// void update_L_second(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); + // void update_L_first(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); + // void update_L_second(Eigen::SparseMatrix<FPP, RowMajor> & L_vec_p, Eigen::SparseMatrix<FPP, RowMajor>& L_vec_q, const FPP2& c, const FPP2& s, int p, int q, Faust::MatSparse<FPP,DEVICE> & L); /** * \brief Algo. step 2.5. * @@ -218,14 +211,14 @@ namespace Faust { virtual void update_err()=0; /** - * Order (sort) D ascendantly into ordered_D and keeps ordered indices in ord_indices. + * Sort D in descending order into ordered_D and keeps ordered indices in ord_indices. */ void order_D(); /** - * Order D into ascendantly or descendantly into order_D and keeps ordered indices in ord_indices. + * Sort D into order_D and keeps ordered indices in ord_indices. * - * \param order -1 to order descendantly, 1 to order ascendantly. + * \param order -1 for descending order, 1 to ascending order. */ void order_D(int order); @@ -260,15 +253,13 @@ namespace Faust { */ const Faust::Vect<FPP,DEVICE>& get_D(const int ord=0); - template<typename FPP3> - void get_Dspm(Faust::MatSparse<FPP3,DEVICE>&, const bool ord=false); - /** \brief Returns the diagonal vector as a sparse matrix. * * \note This function makes copies, is not intented for repeated use (use get_D() for optimized calls). * **/ -// const Faust::MatSparse<FPP,DEVICE> get_Dspm(const bool ord=false); + template<typename FPP3> + void get_Dspm(Faust::MatSparse<FPP3,DEVICE>&, const bool ord=false); /** * Returns the diagonal by copying it in the buffer diag_data (should be allocated by the callee). diff --git a/wrapper/python/pyfaust/factparams.py b/wrapper/python/pyfaust/factparams.py index c7fcf61eeb83a67132b506221e0b29d31e36a5a4..2528ae80f9c5a8335f0cb39707336da9baf9cdf2 100644 --- a/wrapper/python/pyfaust/factparams.py +++ b/wrapper/python/pyfaust/factparams.py @@ -285,7 +285,8 @@ class ConstraintMat(ConstraintGeneric): np.ndarray)): raise TypeError('ConstraintMat must receive a numpy matrix as cons_value ' 'argument.') - self.cons_value = self._cons_value + self.cons_value = np.asfortranarray(self._cons_value) + self._cons_value = self.cons_value if(normalized == None): if(self.name.name == ContraintName.CONST): # for const proj the default is to not normalize @@ -335,7 +336,7 @@ class ConstraintReal(ConstraintGeneric): >>> from numpy.linalg import norm >>> cons = ConstraintReal('normcol', 10, 10, 2.) # a short for ConstraintReal(ConstraintName(ConstraintName.NORMCOL), 10, 10, 2.) >>> M = rand(10,10)*10 - >>> norm(M[:,2],2) + >>> norm(M[:,2]) 17.041462424512272 >>> norm(cons.project(M)[:,2]) 2.0 diff --git a/wrapper/python/src/_FaustCorePy.pyx b/wrapper/python/src/_FaustCorePy.pyx index c7e36fe8afa4a72f8fe673962eb94a3f00c6e75e..4e4e3667d78999cf1905234925b0ff9e5567f040 100644 --- a/wrapper/python/src/_FaustCorePy.pyx +++ b/wrapper/python/src/_FaustCorePy.pyx @@ -858,10 +858,9 @@ cdef class ConstraintMatCore: 'float16', 'float32', 'float64', 'double'] - order = 'C' - if(np.isfortran(M)): - order = 'FORTRAN' - M_out = np.empty(M.shape, dtype=M.dtype, order=order) + M = np.asfortranarray(M) + + M_out = np.empty(M.shape, dtype=M.dtype, order='F') check_matrix(isReal, M) @@ -902,12 +901,11 @@ cdef class ConstraintRealCore: 'float16', 'float32', 'float64', 'double'] - - - order = 'C' - if(np.isfortran(M)): - order = 'FORTRAN' - M_out = np.empty(M.shape, dtype=M.dtype, order=order) +# order = 'C' +# if(np.isfortran(M)): +# order = 'FORTRAN' + M = np.asfortranarray(M) + M_out = np.empty(M.shape, dtype=M.dtype, order='F') check_matrix(isReal, M) if(isReal):