ajout des news fichiers

src/algorithm/factorization/faust_HierarchicalFact.h

+#ifndef __FAUST_HIERARCHICAL_FACT_H__
+#define __FAUST_HIERARCHICAL_FACT_H__
+
+#include "faust_constant.h"
+#include <vector>
+#include "faust_Palm4MSA.h"
+#include "faust_Params.h"
+
+#ifdef __COMPILE_FPPIMERS__
+  #include "faust_Timer.h"
+#endif
+
+/*! \class HierarchicalFact
+   * \brief template class implementing hierarchical factorisation algorithm
+    : <br>  factorization of a data matrix into multiple factors using Faust::Palm4MSA algorithm mixed with a hierarchical approach.
+   *\tparam FPP scalar numeric type, e.g float or double
+   */
+
+template<Device DEVICE> class BlasHandle;
+template<Device DEVICE> class SpBlasHandle;
+
+
+namespace Faust
+{
+
+    template<typename FPP,Device DEVICE> class ConstraintGeneric;
+    template<typename FPP,Device DEVICE> class Palm4MSA;
+    template<typename FPP,Device DEVICE> class MatDense;
+    template<typename FPP,Device DEVICE> class MatSparse;
+    template<typename FPP,Device DEVICE> class Transform;
+
+    template<typename FPP> class StoppingCriterion;
+
+
+    template<typename FPP,Device DEVICE>
+    class HierarchicalFact
+    {
+
+       public:
+
+          HierarchicalFact(const Faust::Params<FPP,DEVICE>& params_, BlasHandle<DEVICE> cublasHandle, SpBlasHandle<DEVICE> cusparseHandle);
+          void get_facts(Faust::Transform<FPP,DEVICE> &)const;
+          void get_facts(std::vector<Faust::MatSparse<FPP,DEVICE> >&)const;
+          void get_facts(std::vector<Faust::MatDense<FPP,DEVICE> >& fact)const{fact = palm_global.get_facts();}
+          void compute_facts();
+          FPP get_lambda()const{return palm_global.get_lambda();}
+          const std::vector<std::vector< FPP> >& get_errors()const;
+
+
+        private:
+          void init();
+          void next_step();
+          void compute_errors();
+
+
+        private:
+          const std::vector< std::vector<const Faust::ConstraintGeneric<FPP,DEVICE>*> > cons;
+          bool isUpdateWayR2L;
+          bool isFactSideLeft;
+          bool isVerbose;
+          int ind_fact ; //indice de factorisation (!= Faust::Palm4MSA::ind_fact : indice de facteur)
+          int nb_fact; // nombre de factorisations (!= Faust::Palm4MSA::nb_fact : nombre de facteurs)
+          Faust::Palm4MSA<FPP,DEVICE> palm_2;
+          Faust::Palm4MSA<FPP,DEVICE> palm_global;
+          const FPP default_lambda; // initial value of lambda for factorization into two factors
+          //std::vector<Faust::MatDense<FPP,DEVICE> > S;
+          std::vector<const Faust::ConstraintGeneric<FPP,DEVICE>*> cons_tmp_global;
+          bool isFactorizationComputed;
+          std::vector<std::vector<FPP> > errors;
+          static const char * class_name;
+          BlasHandle<DEVICE> cublas_handle;
+          SpBlasHandle<DEVICE> cusparse_handle;
+
+
+    #ifdef __COMPILE_TIMERS__
+       public:
+          static Faust::Timer t_init;
+          static Faust::Timer t_next_step;
+
+        void print_timers()const;
+        //void print_prox_timers()const;
+    #endif
+
+
+
+    };
+
+}
+
+#include "faust_HierarchicalFact.hpp"
+
+#endif

src/algorithm/factorization/faust_HierarchicalFact.hpp

+#ifndef __HIERARCHICAL_FACT_CU_HPP__
+#define __HIERARCHICAL_FACT_CU_HPP__
+
+//#include "faust_HierarchicalFact.h"
+#ifdef __COMPILE_TIMERS__
+    #include "faust_Timer.h"
+#endif
+
+#ifdef __COMPILE_GPU__
+	#include "faust_MatSparse_gpu.h"
+	#include "faust_Transform_gpu.h"
+#else
+	#include "faust_MatSparse.h"
+	#include "faust_Transform.h"
+#endif
+
+#include "faust_exception.h"
+using namespace std;
+
+//Faust::HierarchicalFact::Faust::HierarchicalFact(){} // voir avec Luc les parametres par defaut
+
+template<typename FPP,Device DEVICE>
+const char * Faust::HierarchicalFact<FPP,DEVICE>::class_name="Faust::HierarchicalFact";
+
+template<typename FPP,Device DEVICE>
+Faust::HierarchicalFact<FPP,DEVICE>::HierarchicalFact(const Faust::Params<FPP,DEVICE>& params_, BlasHandle<DEVICE> cublasHandle, SpBlasHandle<DEVICE> cusparseHandle):
+   ind_fact(0),
+   cons(params_.cons),
+   isUpdateWayR2L(params_.isUpdateWayR2L),
+   isFactSideLeft(params_.isFactSideLeft),
+   isVerbose(params_.isVerbose),
+   nb_fact(params_.nb_fact-1),
+   palm_2(Palm4MSA<FPP,DEVICE>(params_, cublasHandle, false)),
+   palm_global(Palm4MSA<FPP,DEVICE>(params_, cublasHandle, true)),
+   cons_tmp_global(vector<const Faust::ConstraintGeneric<FPP,DEVICE>*>()),
+   default_lambda(params_.init_lambda),
+   isFactorizationComputed(false),
+   errors(std::vector<std::vector<FPP> >(2,std::vector<FPP >(params_.nb_fact-1,0.0))),
+   cublas_handle(cublasHandle),
+   cusparse_handle(cusparseHandle){}
+
+
+template<typename FPP,Device DEVICE>
+void Faust::HierarchicalFact<FPP,DEVICE>::init()
+{
+#ifdef __COMPILE_TIMERS__
+t_init.start();
+#endif
+
+   cons_tmp_global.clear();
+   if(isFactSideLeft)
+      cons_tmp_global.push_back(cons[0][ind_fact]);
+   else
+      cons_tmp_global.push_back(cons[1][ind_fact]);
+
+
+   palm_global.set_constraint(cons_tmp_global);
+   palm_global.init_fact(1);
+
+
+#ifdef __COMPILE_TIMERS__
+t_init.stop();
+#endif
+
+}
+
+template<typename FPP,Device DEVICE>
+void Faust::HierarchicalFact<FPP,DEVICE>::next_step()
+{
+#ifdef __COMPILE_TIMERS__
+t_next_step.start();
+#endif
+
+   if(isFactorizationComputed)
+   {
+      handleError(class_name,"next_step : factorization has already been computed");
+   }
+
+   vector<const Faust::ConstraintGeneric<FPP,DEVICE>*> cons_tmp_2(2);
+   cons_tmp_2[0]=cons[0][ind_fact];
+   cons_tmp_2[1]=cons[1][ind_fact];
+
+
+   palm_2.set_constraint(cons_tmp_2);
+
+   palm_2.init_fact(2);
+
+   palm_2.set_lambda(default_lambda);
+
+#ifdef __COMPILE_TIMERS__
+palm_2.init_local_timers();
+#endif
+   //while(palm_2.do_continue())
+    //  palm_2.next_step();
+	palm_2.compute_facts();
+
+
+#ifdef __COMPILE_TIMERS__
+palm_2.print_local_timers();
+#endif
+   palm_global.update_lambda_from_palm(palm_2);
+
+
+   if (isFactSideLeft)
+   {
+      cons_tmp_global[0]=cons[0][ind_fact];
+      typename vector<const Faust::ConstraintGeneric<FPP,DEVICE>*>::iterator it;
+      it = cons_tmp_global.begin();
+      cons_tmp_global.insert(it+1,cons[1][ind_fact]);
+   }
+   else
+   {
+      typename vector<const Faust::ConstraintGeneric<FPP,DEVICE>*>::iterator it;
+      it = cons_tmp_global.begin();
+      cons_tmp_global.insert(it+ind_fact,cons[0][ind_fact]);
+      cons_tmp_global[ind_fact+1]=cons[1][ind_fact];
+   }
+
+   palm_global.set_constraint(cons_tmp_global);
+
+
+   palm_global.init_fact_from_palm(palm_2, isFactSideLeft);
+
+#ifdef __COMPILE_TIMERS__
+palm_global.init_local_timers();
+#endif
+   //while(palm_global.do_continue())
+    //  palm_global.next_step();
+	palm_global.compute_facts();
+#ifdef __COMPILE_TIMERS__
+palm_global.print_local_timers();
+#endif
+
+   palm_2.set_data(palm_global.get_res(isFactSideLeft, ind_fact));
+
+
+   compute_errors();
+
+
+   ind_fact++;
+
+
+#ifdef __COMPILE_TIMERS__
+t_next_step.stop();
+palm_2.print_prox_timers();
+#endif
+}
+
+template<typename FPP,Device DEVICE>
+void Faust::HierarchicalFact<FPP,DEVICE>::get_facts(Faust::Transform<FPP,DEVICE> & fact)const
+{
+	std::vector<Faust::MatSparse<FPP,DEVICE> > spfacts;
+	get_facts(spfacts);
+	Faust::Transform<FPP,DEVICE> res(spfacts);
+	fact = res;
+}
+
+template<typename FPP,Device DEVICE>
+void Faust::HierarchicalFact<FPP,DEVICE>::get_facts(std::vector<Faust::MatSparse<FPP,DEVICE> >& sparse_facts)const
+{
+   /*if(!isFactorizationComputed)
+   {
+      cerr << "Error in Faust::HierarchicalFact<FPP,DEVICE>::get_facts : factorization has not been computed" << endl;
+      exit(EXIT_FAILURE);
+   }*/
+
+   const std::vector<Faust::MatDense<FPP,DEVICE> >& full_facts = palm_global.get_facts();
+   sparse_facts.resize(full_facts.size());
+   for (int i=0 ; i<sparse_facts.size() ; i++)
+{
+      sparse_facts[i].init(full_facts[i],cusparse_handle);
+}
+}
+
+
+
+template<typename FPP,Device DEVICE>
+void Faust::HierarchicalFact<FPP,DEVICE>::compute_facts()
+{
+   if(isFactorizationComputed)
+   {
+      handleError(class_name,"compute_facts : factorization has already been computed");
+   }
+
+  init();
+  for (int i=0 ; i<=nb_fact-1 ; i++)
+  {
+     cout << "Faust::HierarchicalFact<FPP,DEVICE>::compute_facts : factorisation "<<i+1<<"/"<<nb_fact <<endl;
+     next_step();
+  }
+
+  isFactorizationComputed = true;
+
+}
+
+
+template<typename FPP,Device DEVICE>
+const std::vector<std::vector< FPP> >& Faust::HierarchicalFact<FPP,DEVICE>::get_errors()const
+{
+    if(!isFactorizationComputed)
+    {
+        handleError(class_name,"get_errors() : Factorization has not been computed");
+    }
+    return errors;
+}
+
+template<typename FPP,Device DEVICE>
+void Faust::HierarchicalFact<FPP,DEVICE>::compute_errors()
+{
+   vector<Faust::MatSparse<FPP,DEVICE> > sp_facts;
+   get_facts(sp_facts);
+
+
+
+   Faust::Transform<FPP,DEVICE> faust_Transform_tmp(sp_facts, get_lambda());
+   const Faust::MatDense<FPP,DEVICE> estimate_mat = faust_Transform_tmp.get_product(cublas_handle, cusparse_handle);
+
+   Faust::MatDense<FPP,DEVICE> data(palm_global.get_data());
+
+   FPP data_norm = data.norm();
+
+
+   data -= estimate_mat;
+
+   errors[0][ind_fact] =  estimate_mat.norm()/data_norm;
+   errors[1][ind_fact] =  faust_Transform_tmp.get_total_nnz()/data.getNbRow()/data.getNbCol();
+
+}
+
+
+#ifdef __COMPILE_TIMERS__
+template<typename FPP,Device DEVICE> Faust::Timer Faust::HierarchicalFact<FPP,DEVICE>::t_init;
+template<typename FPP,Device DEVICE> Faust::Timer Faust::HierarchicalFact<FPP,DEVICE>::t_next_step;
+
+template<typename FPP,Device DEVICE>
+void Faust::HierarchicalFact<FPP,DEVICE>::print_timers()const
+{
+   palm_global.print_global_timers();
+   cout << "timers in Faust::HierarchicalFact :" << endl;
+   cout << "t_init      = " << t_init.get_time()      << " s for "<< t_init.get_nb_call()      << " calls" << endl;
+   cout << "t_next_step = " << t_next_step.get_time() << " s for "<< t_next_step.get_nb_call() << " calls" << endl<<endl;
+
+Faust::MatDense<FPP,DEVICE> M_tmp;
+Faust::MatSparse<FPP,DEVICE> S_tmp;
+M_tmp.print_timers();
+S_tmp.print_timers();
+
+}
+#endif
+
+#endif

src/algorithm/factorization/faust_Palm4MSA.h

+#ifndef __FAUST_PALM4MSA_H__
+#define __FAUST_PALM4MSA_H__
+
+#include <iostream>
+#include "faust_constant.h"
+#include <vector>
+#include "faust_StoppingCriterion.h"
+
+
+#ifdef __COMPILE_TIMERS__
+    #include "faust_Timer.h"
+#endif
+
+#include "faust_MatDense.h"
+
+template<Device DEVICE> class BlasHandle;
+
+namespace Faust
+{
+
+
+    template<typename FPP,Device DEVICE> class MatDense;
+    template<typename FPP,Device DEVICE> class Transform;
+
+    template<typename FPP,Device DEVICE> class ConstraintGeneric;
+    template<typename FPP,Device DEVICE> class Params;
+    template<typename FPP,Device DEVICE> class ParamsPalm;
+    template<typename FPP> class StoppingCriterion;
+
+
+    /*! \class Palm4MSA
+       * \brief template class implementing Palm4MSA (PALM for Multi-layer Sparse Approximation) factorization algorithm
+        : <br>
+        factorization of a data matrix (dense) into multiple factors (sparse) using PALM
+       *
+       *\tparam FPP scalar numeric type, e.g float or double
+       */
+
+
+    template<typename FPP,Device DEVICE>
+    class Palm4MSA
+    {
+
+       public:
+            /*!
+         *  \brief
+         * initialize Palm4MSA from Faust::Params (HierarchicalFact parameter)
+         *\tparam isGlobal_ : if true, the Palm4MSA StoppingCriterion stop_crit attribute is initialize from params_.stop_crit_global <br> and if false, it is initialize from stop_crit_2facts
+         */
+          Palm4MSA(const Faust::Params<FPP,DEVICE>& params_, const BlasHandle<DEVICE> blasHandle, const bool isGlobal_);
+          Palm4MSA(const Faust::ParamsPalm<FPP,DEVICE>& params_palm_, const BlasHandle<DEVICE> blasHandle, const bool isGlobal_=false);
+
+          void set_constraint(const std::vector<const Faust::ConstraintGeneric<FPP,DEVICE>*> const_vec_){const_vec=const_vec_;isConstraintSet=true;}
+          void set_data(const Faust::MatDense<FPP,DEVICE>& data_){data=data_;}
+          void set_lambda(const FPP lambda_){lambda = lambda_;}
+
+          /*!
+         *  \brief
+         * useful in HierarchicalFact, update lambda of palm_global from palm_2
+         */
+          void update_lambda_from_palm(const Palm4MSA& palm){lambda *= palm.lambda;}
+
+          /*!
+         *  \brief
+         * compute the factorisation
+         */
+          void compute_facts();
+
+          /*!
+         *  \brief
+         * return the multiplicative scalar lambda
+         */
+          FPP get_lambda()const{return lambda;}
+
+          FPP get_RMSE()const{return error.norm()/sqrt((double)(data.getNbRow()*data.getNbCol()));}
+          const Faust::MatDense<FPP,DEVICE>& get_res(bool isFactSideLeft_, int ind_)const{return isFactSideLeft_ ? S[0] : S[ind_+1];}
+          const Faust::MatDense<FPP,DEVICE>& get_data()const{return data;}
+          void get_facts(Faust::Transform<FPP,DEVICE> & faust_fact) const;
+
+
+              /*!
+         *  \brief
+         * initialize the factors to the default value,
+         * the first factor to be factorised is set to zero matrix
+         * whereas all the other are set to identity
+         */
+          void init_fact(int nb_facts_);
+          void next_step();
+          bool do_continue(){bool cont=stop_crit.do_continue(++ind_ite); if(!cont){ind_ite=-1;isConstraintSet=false;}return cont;} // CAUTION !!! pre-increment of ind_ite: the value in stop_crit.do_continue is ind_ite+1, not ind_ite
+          //bool do_continue()const{return stop_crit.do_continue(++ind_ite, error);};
+
+          /*!
+         *  \brief
+         * useful in HierarchicalFact, update the factors of palm_global from palm_2
+         */
+          void init_fact_from_palm(const Palm4MSA& palm, bool isFactSideLeft);
+          const std::vector<Faust::MatDense<FPP,DEVICE> >& get_facts()const {return S;}
+
+          ~Palm4MSA(){}
+
+       private:
+          void check_constraint_validity();
+          void compute_c();
+          void compute_grad_over_c();
+          void compute_projection();
+          void update_L();
+          void update_R();
+          void compute_lambda();
+          static const char * class_name;
+          static const FPP lipschitz_multiplicator;
+
+       public:
+          Faust::StoppingCriterion<FPP> stop_crit;
+
+
+       private:
+          // modif AL AL
+          Faust::MatDense<FPP,DEVICE> data;
+
+          FPP lambda;
+          int nb_fact; // number of factors
+          std::vector<Faust::MatDense<FPP,DEVICE> > S; // contains S_0^i, S_1^i, ...
+
+          // RorL_vec matches R if (!isUpdateWayR2L)
+          // RorL_vec matches L if (isUpdateWayR2L)
+          std::vector<Faust::MatDense<FPP,DEVICE> > RorL;
+          // LorR_mat matches L if (!isUpdateWayR2L)
+          // LorR_mat matches R if (isUpdateWayR2L)
+
+          // modif AL AL
+          Faust::MatDense<FPP,DEVICE> LorR;
+          //Faust::MatDense<FPP,DEVICE>& LorR;
+
+
+
+          std::vector<const Faust::ConstraintGeneric<FPP,DEVICE>*> const_vec; // vector of constraints of size nfact
+          int ind_fact; //indice de facteur (!= HierarchicalFact::ind_fact : indice de factorisation)
+          int ind_ite;
+          // FPP lipschitz_multiplicator;
+          const bool verbose;
+          const bool isUpdateWayR2L;
+          const bool isConstantStepSize;
+          bool isCComputed;
+          bool isGradComputed;
+          bool isProjectionComputed;
+          bool isLastFact;
+          bool isConstraintSet;
+          const bool isGlobal;
+          bool isInit; // only used for global factorization (if isGlobal)
+          Faust::MatDense<FPP,DEVICE> grad_over_c;
+          FPP c;
+          Faust::MatDense<FPP,DEVICE> error; // error = lambda*L*S*R - data
+          BlasHandle<DEVICE> blas_handle;
+
+
+
+
+
+
+    #ifdef __COMPILE_TIMERS__
+       public:
+          Faust::Timer_gpu t_local_compute_projection;
+          Faust::Timer_gpu t_local_compute_grad_over_c;
+          Faust::Timer_gpu t_local_compute_c;
+          Faust::Timer_gpu t_local_compute_lambda;
+          Faust::Timer_gpu t_local_update_R;
+          Faust::Timer_gpu t_local_update_L;
+          Faust::Timer_gpu t_local_check;
+          Faust::Timer_gpu t_local_init_fact;
+          Faust::Timer_gpu t_local_next_step;
+          Faust::Timer_gpu t_local_init_fact_from_palm;
+
+
+          static Faust::Timer_gpu t_global_compute_projection;
+          static Faust::Timer_gpu t_global_compute_grad_over_c;
+          static Faust::Timer_gpu t_global_compute_c;
+          static Faust::Timer_gpu t_global_compute_lambda;
+          static Faust::Timer_gpu t_global_update_R;
+          static Faust::Timer_gpu t_global_update_L;
+          static Faust::Timer_gpu t_global_check;
+          static Faust::Timer_gpu t_global_init_fact;
+          static Faust::Timer_gpu t_global_next_step;
+          static Faust::Timer_gpu t_global_init_fact_from_palm;
+
+          static Faust::Timer_gpu t_prox_const;
+          static Faust::Timer_gpu t_prox_sp;
+          static Faust::Timer_gpu t_prox_spcol;
+          static Faust::Timer_gpu t_prox_splin;
+          static Faust::Timer_gpu t_prox_normcol;
+
+          static int nb_call_prox_const;
+          static int nb_call_prox_sp;
+          static int nb_call_prox_spcol;
+          static int nb_call_prox_splin;
+          static int nb_call_prox_normcol;
+
+
+
+
+
+
+
+       void init_local_timers();
+
+       void print_global_timers()const;
+       void print_local_timers()const;
+       void print_prox_timers() const;
+    #endif
+
+    };
+
+}
+
+#include "faust_Palm4MSA.hpp"
+
+
+#endif
+
+

src/algorithm/factorization/faust_StoppingCriterion.h

+#ifndef __FAUST_STOPPING_CRITERION__
+#define __FAUST_STOPPING_CRITERION__
+
+#include "faust_constant.h"
+#include <iostream>
+
+
+namespace Faust
+{
+
+    template<typename T>
+    class StoppingCriterion
+    {
+       public:
+          StoppingCriterion():
+             isCriterionError(false),
+             nb_it(500){}
+
+          StoppingCriterion(int nb_it_):
+             isCriterionError(false),nb_it(nb_it_){check_validity();}
+
+          StoppingCriterion(T errorThreshold_, int maxIteration_=10000):
+             isCriterionError(true),errorThreshold(errorThreshold_),maxIteration(maxIteration_){check_validity();}
+
+          StoppingCriterion(bool isCriterionError_);
+
+          ~StoppingCriterion(){}
+
+          bool do_continue(int current_ite, T error=-2.0)const;
+          int get_crit() const{return nb_it;}
+       private:
+          void check_validity()const;
+
+       private:
+          // if isCriterionError then criterion is error else criterion is number of iteration
+          bool  isCriterionError;
+          int nb_it;   // number of iterations if !isCriterionError
+          T errorThreshold;
+          int maxIteration;
+          // only used as stopping criterion, if isCriterionError, when error is still greater than
+          static const char * class_name;
+    };
+
+}
+
+#include "faust_StoppingCriterion.hpp"
+
+#endif

src/algorithm/factorization/faust_StoppingCriterion.hpp

+#ifndef __STOPPING_CRITERION_HPP__
+#define __STOPPING_CRITERION_HPP__
+
+//#include "faust_StoppingCriterion.h"
+#include <iostream>
+#include <cstdlib>
+#include "faust_exception.h"
+
+template<typename T>
+const char * Faust::StoppingCriterion<T>::class_name="Faust::StoppingCriterion::";
+
+template<typename T>
+Faust::StoppingCriterion<T>::StoppingCriterion(bool isCriterionError_) : isCriterionError(isCriterionError_)
+{
+   if (isCriterionError_)
+   {
+      errorThreshold = 0.3;
+      maxIteration = 10000;
+   }
+   else
+      nb_it = 500;
+}
+
+template<typename T>
+void Faust::StoppingCriterion<T>::check_validity()const
+{
+   if (isCriterionError)
+   {
+      if (errorThreshold>1 || maxIteration < 0)
+      {
+        handleError(class_name,"check_validity : errorThreshold must be strictly greater than 1 and maxIteration must be strictly positive");
+      }
+   }
+   else if (nb_it < 0)
+   {
+     handleError(class_name,"::check_validity : nb_it must be positive");
+   }
+}
+
+// current_ite in zero-based indexing
+template<typename T>
+bool Faust::StoppingCriterion<T>::do_continue(int current_ite, T current_error /* = -2.0 */)const
+{
+
+   if (!isCriterionError) // if criterion is number of iteration, current_error does not matter
+      return current_ite<nb_it ? true : false;
+   else if (isCriterionError && current_error != -2.0)
+      if (current_error < errorThreshold)
+         return false;
+      else if (current_ite <  maxIteration) // and current_error >= errorThreshold
+         return true;
+      else // if current_error >= errorThreshold and current_ite >= maxIteration
+      {
+         std::cerr << "warning in Faust::StoppingCriterion<T>::do_continue : number of maximum iterations has been reached and current error is still greater than the threshold" << std::endl;
+         return true;
+      }
+   else // if criterion is error and current_error has not been initialized
+   {
+     handleError(class_name,"check_validity : when stopping criterion is error, the current error needs to be given as second parameter");
+   }
+}
+
+#endif