From d1afe518c11c7c49a9da002f9da403fd843f2e2b Mon Sep 17 00:00:00 2001
From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com>
Date: Thu, 22 May 2025 08:58:46 -0400
Subject: [PATCH 1/3] rename omp to ompursuit internally to prevent naming
 clash in c++20

---
 spams/spams.py                   | 4 ++--
 spams_wrap/spams.h               | 4 ++--
 spams_wrap/spams/decomp/decomp.h | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/spams/spams.py b/spams/spams.py
index 4669f2dd..62e50b9f 100644
--- a/spams/spams.py
+++ b/spams/spams.py
@@ -185,9 +185,9 @@ def omp(X,D,L=None,eps= None,lambda1 = None,return_reg_path = False, numThreads
         if str(type(lambda1)) != "<type 'numpy.ndarray'>":
             lambda1 = np.array([lambda1],dtype=X.dtype)
     if return_reg_path:
-        ((indptr,indices,data,shape),path) = spams_wrap.omp(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
+        ((indptr,indices,data,shape),path) = spams_wrap.ompursuit(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
     else:
-        (indptr,indices,data,shape) = spams_wrap.omp(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
+        (indptr,indices,data,shape) = spams_wrap.ompursuit(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
     alpha = ssp.csc_matrix((data,indices,indptr),shape)
     if return_reg_path:
         return (alpha,path)
diff --git a/spams_wrap/spams.h b/spams_wrap/spams.h
index f3e980c1..ce2a4cf4 100644
--- a/spams_wrap/spams.h
+++ b/spams_wrap/spams.h
@@ -279,10 +279,10 @@ SpMatrix<T> *_omp(Matrix<T> *X,Matrix<T> *D,Matrix<T> **path,bool return_reg_pat
     if(return_reg_path) {
       *path = new Matrix<T>(K,scalar_L);
       (*path)->setZeros();
-      omp((Matrix<T> &)(*X),(Matrix<T> &)(*D),(SpMatrix<T> &)(*alpha),pL,pE,pLambda,vecL,vecEps,vecLambda,numThreads,*path);
+      ompursuit((Matrix<T> &)(*X),(Matrix<T> &)(*D),(SpMatrix<T> &)(*alpha),pL,pE,pLambda,vecL,vecEps,vecLambda,numThreads,*path);
     } else {
       *path = NULL;
-      omp((Matrix<T> &)(*X),(Matrix<T> &)(*D),(SpMatrix<T> &)(*alpha),pL,pE,pLambda,vecL,vecEps,vecLambda,numThreads);
+      ompursuit((Matrix<T> &)(*X),(Matrix<T> &)(*D),(SpMatrix<T> &)(*alpha),pL,pE,pLambda,vecL,vecEps,vecLambda,numThreads);
     }
     return alpha;
 }
diff --git a/spams_wrap/spams/decomp/decomp.h b/spams_wrap/spams/decomp/decomp.h
index 2c8ccec5..a4e5f4f0 100644
--- a/spams_wrap/spams/decomp/decomp.h
+++ b/spams_wrap/spams/decomp/decomp.h
@@ -51,7 +51,7 @@ static char nonUnit='n';
 ///   * optimized for a large number of signals (precompute the Gramm matrix
 
 template <typename T>
-void omp(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
+void ompursuit(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
       const int *L, const T* eps, const T* lambda, const bool vecL = false,
       const bool vecEps = false, const bool Lambda=false, const int numThreads=-1,
       Matrix<T>* path = NULL);
@@ -320,7 +320,7 @@ void coreSOMP(const Matrix<T>& X, const Matrix<T>& D, const Matrix<T>& G,
 ///   * optimized for a big number of signals (precompute the Gramm matrix
 
 template <typename T>
-void omp(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
+void ompursuit(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
       const int* pL, const T* peps, const T* pLambda,
       const bool vecL, const bool vecEps,
       const bool vecLambda, const int numThreads, Matrix<T>* path) {

From 6d1d7c61b018e74c362b231b4584a4f1d22fe03f Mon Sep 17 00:00:00 2001
From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com>
Date: Thu, 22 May 2025 10:44:17 -0400
Subject: [PATCH 2/3] fix templating warning for c++20

---
 spams/spams.py                     |   4 +-
 spams_wrap/spams/dictLearn/dicts.h |  89 +++---
 spams_wrap/spams/linalg/linalg.h   | 360 ++++++++++++------------
 spams_wrap/spams/prox/fista.h      | 436 ++++++++++++++---------------
 4 files changed, 444 insertions(+), 445 deletions(-)

diff --git a/spams/spams.py b/spams/spams.py
index 62e50b9f..4669f2dd 100644
--- a/spams/spams.py
+++ b/spams/spams.py
@@ -185,9 +185,9 @@ def omp(X,D,L=None,eps= None,lambda1 = None,return_reg_path = False, numThreads
         if str(type(lambda1)) != "<type 'numpy.ndarray'>":
             lambda1 = np.array([lambda1],dtype=X.dtype)
     if return_reg_path:
-        ((indptr,indices,data,shape),path) = spams_wrap.ompursuit(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
+        ((indptr,indices,data,shape),path) = spams_wrap.omp(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
     else:
-        (indptr,indices,data,shape) = spams_wrap.ompursuit(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
+        (indptr,indices,data,shape) = spams_wrap.omp(X,D,return_reg_path,given_L,L,given_eps,eps,given_lambda1,lambda1,numThreads)
     alpha = ssp.csc_matrix((data,indices,indptr),shape)
     if return_reg_path:
         return (alpha,path)
diff --git a/spams_wrap/spams/dictLearn/dicts.h b/spams_wrap/spams/dictLearn/dicts.h
index 5a9f08c6..a2cb184e 100644
--- a/spams_wrap/spams/dictLearn/dicts.h
+++ b/spams_wrap/spams/dictLearn/dicts.h
@@ -1,5 +1,5 @@
 
-/* Software SPAMS v2.1 - Copyright 2009-2011 Julien Mairal 
+/* Software SPAMS v2.1 - Copyright 2009-2011 Julien Mairal
  *
  * This file is part of SPAMS.
  *
@@ -70,7 +70,7 @@ void regul_error(char *buffer, int bufsize,const char *message) {
     // calculate size
     for(unsigned int i = 0;i < NBREGUL;i++)
       size += strlen(regul_table[i].name) + 1;
-  } 
+  }
   if (size >= bufsize) {
     n1 = bufsize - 1;
     strncpy(buffer,"Invalid regularization\n",n1);
@@ -90,7 +90,7 @@ void regul_error(char *buffer, int bufsize,const char *message) {
 
 template <typename T> struct ParamDictLearn {
    public:
-      ParamDictLearn() : 
+      ParamDictLearn() :
          mode(PENALTY),
 	   regul(FISTA::RIDGE),
 	   tol(1.e-4),
@@ -131,7 +131,7 @@ template <typename T> struct ParamDictLearn {
          updateD(true),
          updateW(true),
          updateTheta(true),
-         logName(NULL), 
+         logName(NULL),
          iter_updateD(1) { };
       ~ParamDictLearn() { delete[](logName); };
       int iter;
@@ -141,7 +141,7 @@ template <typename T> struct ParamDictLearn {
       T tol;
       bool ista;
       bool fixed_step;
-      bool posAlpha; 
+      bool posAlpha;
       constraint_type_D modeD;
       bool posD;
       mode_compute modeParam;
@@ -192,17 +192,17 @@ template <typename T> class Trainer {
             const int NUM_THREADS=-1);
       /// Constructor with existing structure
       Trainer(const Matrix<T>& A, const Matrix<T>& B, const Matrix<T>& D,
-            const int itercount, const int batchsize, 
+            const int itercount, const int batchsize,
             const int NUM_THREADS);
 
       /// train or retrain using the matrix X
       /* train with graph or tree */
       void train(const Data<T>& X, const ParamDictLearn<T>& param);
       void train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
-		 const GraphStruct<T>* graph_st = NULL, 
+		 const GraphStruct<T>* graph_st = NULL,
 		 const TreeStruct<T>* tree_st = NULL);
       void trainOffline(const Data<T>& X, const ParamDictLearn<T>& param);
-      
+
       /// Accessors
       void getA(Matrix<T>& A) const { A.copy(_A);};
       void getB(Matrix<T>& B) const { B.copy(_B);};
@@ -211,14 +211,14 @@ template <typename T> class Trainer {
 
    private:
       /// Forbid lazy copies
-      explicit Trainer<T>(const Trainer<T>& trainer);
+      explicit Trainer(const Trainer<T>& trainer);
       /// Forbid lazy copies
       Trainer<T>& operator=(const Trainer<T>& trainer);
 
       /// clean the dictionary
       void cleanDict(const Data<T>& X, Matrix<T>& G,
             const bool posD = false,
-            const constraint_type_D modeD = L2, const T gamma1 = 0, 
+            const constraint_type_D modeD = L2, const T gamma1 = 0,
             const T gamma2 = 0,
             const T maxCorrel =
             0.999999);
@@ -238,7 +238,7 @@ template <typename T> class Trainer {
 
 /// Empty constructor
 template <typename T> Trainer<T>::Trainer() : _k(0), _initialDict(false),
-   _itercount(0), _batchsize(256) { 
+   _itercount(0), _batchsize(256) {
       _NUM_THREADS=1;
 #ifdef _OPENMP
       _NUM_THREADS =  MIN(MAX_THREADS,omp_get_num_procs());
@@ -248,9 +248,9 @@ template <typename T> Trainer<T>::Trainer() : _k(0), _initialDict(false),
 
 /// Constructor with data
 template <typename T> Trainer<T>::Trainer(const int k, const
-      int batchsize, const int NUM_THREADS) : _k(k), 
-   _initialDict(false), _itercount(0),_batchsize(batchsize), 
-   _NUM_THREADS(NUM_THREADS) { 
+      int batchsize, const int NUM_THREADS) : _k(k),
+   _initialDict(false), _itercount(0),_batchsize(batchsize),
+   _NUM_THREADS(NUM_THREADS) {
       if (_NUM_THREADS == -1) {
          _NUM_THREADS=1;
 #ifdef _OPENMP
@@ -260,7 +260,7 @@ template <typename T> Trainer<T>::Trainer(const int k, const
    };
 
 /// Constructor with initial dictionary
-template <typename T> Trainer<T>::Trainer(const Matrix<T>& D, 
+template <typename T> Trainer<T>::Trainer(const Matrix<T>& D,
       const int batchsize, const int NUM_THREADS) : _k(D.n()),
      _initialDict(true),_itercount(0),_batchsize(batchsize),
    _NUM_THREADS(NUM_THREADS) {
@@ -277,7 +277,7 @@ template <typename T> Trainer<T>::Trainer(const Matrix<T>& D,
 
 /// Constructor with existing structure
 template <typename T> Trainer<T>::Trainer(const Matrix<T>& A, const Matrix<T>&
-      B, const Matrix<T>& D, const int itercount, const int batchsize, 
+      B, const Matrix<T>& D, const int itercount, const int batchsize,
       const int NUM_THREADS) : _k(D.n()),_initialDict(true),_itercount(itercount),
    _batchsize(batchsize),
     _NUM_THREADS(NUM_THREADS) {
@@ -293,7 +293,7 @@ template <typename T> Trainer<T>::Trainer(const Matrix<T>& A, const Matrix<T>&
    };
 
 template <typename T>
-void Trainer<T>::cleanDict(const Data<T>& X, Matrix<T>& G, 
+void Trainer<T>::cleanDict(const Data<T>& X, Matrix<T>& G,
       const bool posD,
       const constraint_type_D modeD, const T gamma1,
       const T gamma2,
@@ -364,9 +364,9 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
             cout << "Online Dictionary Learning with exponential decay t0: " << t0 << " rho: " << rho << endl;
          }
       }
-      if (param.posD) 
+      if (param.posD)
          cout << "Positivity constraints on D activated" << endl;
-      if (param.posAlpha) 
+      if (param.posAlpha)
          cout << "Positivity constraints on alpha activated" << endl;
       if (param.modeD != L2) cout << "Sparse dictionaries, mode: " << param.modeD << ", gamma1: " << param.gamma1 << ", gamma2: " << param.gamma2 << endl;
       cout << "mode Alpha " << param.mode << endl;
@@ -385,7 +385,7 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
    const int M = X.n();
    const int K = _k;
    const int n = X.m();
-   const int L = param.mode == SPARSITY ? static_cast<int>(param.lambda) : 
+   const int L = param.mode == SPARSITY ? static_cast<int>(param.lambda) :
       param.mode == PENALTY && param.lambda == 0 && param.lambda2 > 0 && !param.posAlpha ? K : MIN(n,K);
    const int batchsize= param.batch ? M : MIN(_batchsize,M);
 
@@ -397,7 +397,7 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
       flush(cout);
    }
 
-   if (_D.m() != n || _D.n() != K) 
+   if (_D.m() != n || _D.n() != K)
       _initialDict=false;
 
    srandom(0);
@@ -509,7 +509,7 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
          flush(cout);
       }
       time.stop();
-      if (param.iter < 0 && 
+      if (param.iter < 0 &&
             time.getElapsed() > T(-param.iter)) break;
       if (param.log) {
          int seconds=static_cast<int>(floor(log(time.getElapsed())*5));
@@ -522,10 +522,10 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
          }
       }
       time.start();
-      
+
       Matrix<T> G;
       _D.XtX(G);
-      if (param.clean) 
+      if (param.clean)
          this->cleanDict(X,G,param.posD,
                param.modeD,param.gamma1,param.gamma2);
       G.addDiag(MAX(param.lambda2,1e-10));
@@ -591,7 +591,7 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
             }
          }
          int count2=0;
-         for (int k = 0; k<L; ++k) 
+         for (int k = 0; k<L; ++k)
             if (ind[k] == -1) {
                break;
             } else {
@@ -696,8 +696,8 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
          Beven.scal(scal);
          Aodd.scal(scal);
          Bodd.scal(scal);
-         if ((_itercount > 0 && i*batchsize < M) 
-               || (_itercount == 0 && t0 != 0 && 
+         if ((_itercount > 0 && i*batchsize < M)
+               || (_itercount == 0 && t0 != 0 &&
                   i*batchsize < 10000)) {
             Aorig.scal(scal);
             Borig.scal(scal);
@@ -738,7 +738,7 @@ void Trainer<T>::train(const Data<T>& X, const ParamDictLearn<T>& param) {
                      newd.normalize2();
                      di.copy(newd);
                   }
-               } else if (param.clean && 
+               } else if (param.clean &&
                      ((_itercount+i)*batchsize) > 10000) {
                   _D.refCol(k,di);
                   di.setZeros();
@@ -786,9 +786,9 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
             cout << "Online Dictionary Learning with exponential decay t0: " << t0 << " rho: " << rho << endl;
          }
       }
-      if (param.posD) 
+      if (param.posD)
          cout << "Positivity constraints on D activated" << endl;
-      if (param.posAlpha) 
+      if (param.posAlpha)
          cout << "Positivity constraints on alpha activated" << endl;
       if (param.modeD != L2) cout << "Sparse dictionaries, mode: " << param.modeD << ", gamma1: " << param.gamma1 << ", gamma2: " << param.gamma2 << endl;
       cout << "mode Alpha " << param.mode << endl;
@@ -807,7 +807,7 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
    const int M = X.n();
    const int K = _k;
    const int n = X.m();
-   const int L = param.mode == SPARSITY ? static_cast<int>(param.lambda) : 
+   const int L = param.mode == SPARSITY ? static_cast<int>(param.lambda) :
       param.mode == PENALTY && param.lambda == 0 && param.lambda2 > 0 && !param.posAlpha ? K : MIN(n,K);
    const int batchsize= param.batch ? M : MIN(_batchsize,M);
 
@@ -819,7 +819,7 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
       flush(cout);
    }
 
-   if (_D.m() != n || _D.n() != K) 
+   if (_D.m() != n || _D.n() != K)
       _initialDict=false;
 
    srandom(0);
@@ -958,7 +958,7 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
          flush(cout);
       }
       time.stop();
-      if (param.iter < 0 && 
+      if (param.iter < 0 &&
             time.getElapsed() > T(-param.iter)) break;
       if (param.log) {
          int seconds=static_cast<int>(floor(log(time.getElapsed())*5));
@@ -971,10 +971,10 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
          }
       }
       time.start();
-      
+
       // !!      Matrix<T> G;
       _D.XtX(G);
-      if (param.clean) 
+      if (param.clean)
          this->cleanDict(X,G,param.posD,
                param.modeD,param.gamma1,param.gamma2);
       G.addDiag(MAX(param.lambda2,1e-10));
@@ -1038,7 +1038,7 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
 	   losses[numT]->init(Xj);
 	   regularizers[numT]->reset();
 	   //	   alpha.refCol(j,alphai);
-	   FISTA::FISTA_Generic(*(losses[numT]),*(regularizers[numT]),alpha0i,alphai,optim_infoi,param_fista); 
+	   FISTA::FISTA_Generic(*(losses[numT]),*(regularizers[numT]),alpha0i,alphai,optim_infoi,param_fista);
 	   alphai.toSparse(spcoeffj);
 	 } else {
             if (param.mode == SPARSITY) {
@@ -1051,7 +1051,7 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
          }
 	 if(param.mode != FISTAMODE) {
 	   INTM count2=0;
-	   for (int k = 0; k<L; ++k) 
+	   for (int k = 0; k<L; ++k)
 	     if (ind[k] == -1) {
 	       break;
 	     } else {
@@ -1157,8 +1157,8 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
          Beven.scal(scal);
          Aodd.scal(scal);
          Bodd.scal(scal);
-         if ((_itercount > 0 && i*batchsize < M) 
-               || (_itercount == 0 && t0 != 0 && 
+         if ((_itercount > 0 && i*batchsize < M)
+               || (_itercount == 0 && t0 != 0 &&
                   i*batchsize < 10000)) {
             Aorig.scal(scal);
             Borig.scal(scal);
@@ -1199,7 +1199,7 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
                      newd.normalize2();
                      di.copy(newd);
                   }
-               } else if (param.clean && 
+               } else if (param.clean &&
                      ((_itercount+i)*batchsize) > 10000) {
                   _D.refCol(k,di);
                   di.setZeros();
@@ -1234,7 +1234,7 @@ void Trainer<T>::train_fista(const Data<T>& X, const ParamDictLearn<T>& param,
 
 
 template <typename T>
-void writeLog(const Matrix<T>& D, const T time, int iter, 
+void writeLog(const Matrix<T>& D, const T time, int iter,
       char* name) {
    std::ofstream f;
    f.precision(12);
@@ -1253,7 +1253,7 @@ void writeLog(const Matrix<T>& D, const T time, int iter,
 
 
 template <typename T>
-void Trainer<T>::trainOffline(const Data<T>& X, 
+void Trainer<T>::trainOffline(const Data<T>& X,
       const ParamDictLearn<T>& param) {
 
    int sparseD = param.modeD == L1L2 ? 2 : 6;
@@ -1335,7 +1335,7 @@ void Trainer<T>::trainOffline(const Data<T>& X,
    for (int i = 0; i<JJ; ++i) {
       if (J < 0 && time.getElapsed() > T(-J)) break;
       _D.XtX(G);
-      if (param.clean) 
+      if (param.clean)
          this->cleanDict(X,G,param.posD,
                param.modeD,param.gamma1,param.gamma2);
       int j;
@@ -1433,7 +1433,7 @@ void Trainer<T>::trainOffline(const Data<T>& X,
       }
    }
    _D.XtX(G);
-   if (param.clean) 
+   if (param.clean)
       this->cleanDict(X,G,param.posD,param.modeD,
             param.gamma1,param.gamma2);
    time.printElapsed();
@@ -1447,4 +1447,3 @@ void Trainer<T>::trainOffline(const Data<T>& X,
 
 
 #endif
-
diff --git a/spams_wrap/spams/linalg/linalg.h b/spams_wrap/spams/linalg/linalg.h
index 169735e1..8f842559 100755
--- a/spams_wrap/spams/linalg/linalg.h
+++ b/spams_wrap/spams/linalg/linalg.h
@@ -1,4 +1,4 @@
-/* Software SPAMS v2.1 - Copyright 2009-2011 Julien Mairal 
+/* Software SPAMS v2.1 - Copyright 2009-2011 Julien Mairal
  *
  * This file is part of SPAMS.
  *
@@ -17,7 +17,7 @@
  */
 
 /* \file
- *                toolbox Linalg 
+ *                toolbox Linalg
  *
  *                by Julien Mairal
  *                julien.mairal@inria.fr
@@ -66,12 +66,12 @@ typedef std::list< INTM > group;
 typedef std::list< group > list_groups;
 typedef std::vector< group > vector_groups;
 
-template <typename T> 
+template <typename T>
 static inline bool isZero(const T lambda) {
    return static_cast<double>(abs<T>(lambda)) < 1e-99;
 }
 
-template <typename T> 
+template <typename T>
 static inline bool isEqual(const T lambda1, const T lambda2) {
    return static_cast<double>(abs<T>(lambda1-lambda2)) < 1e-99;
 }
@@ -145,23 +145,23 @@ template <typename T> class AbstractMatrixB {
             const T alpha = 1.0, const T beta = 0.0) const = 0;
 
       /// perform b = alpha*A*x + beta*b, when x is sparse
-      virtual void mult(const SpVector<T>& x, Vector<T>& b, 
+      virtual void mult(const SpVector<T>& x, Vector<T>& b,
             const T alpha = 1.0, const T beta = 0.0) const = 0;
 
-      virtual void mult(const Vector<T>& x, Vector<T>& b, 
+      virtual void mult(const Vector<T>& x, Vector<T>& b,
             const T alpha = 1.0, const T beta = 0.0) const = 0;
 
       /// perform C = a*A*B + b*C, possibly transposing A or B.
-      virtual void mult(const Matrix<T>& B, Matrix<T>& C, 
+      virtual void mult(const Matrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const = 0;
 
-      virtual void mult(const SpMatrix<T>& B, Matrix<T>& C, 
+      virtual void mult(const SpMatrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const = 0;
 
       /// perform C = a*B*A + b*C, possibly transposing A or B.
-      virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C, 
+      virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const = 0;
 
@@ -177,8 +177,8 @@ template <typename T> class AbstractMatrixB {
 
       void ridgeCG(const Vector<T>& b, Vector<T>& x, const T lambda, const T tol, const int itermax) const;
       void ridgeCG(const Vector<T>& b, const Vector<T>& delta, Vector<T>& x, const T lambda, const T tol, const int itermax) const;
-      void ridgeCG(const Matrix<T>& b, Matrix<T>& x, const T lambda, const T tol, const int itermax, const int numThreads = -1) const; 
-      void ridgeCG(const Matrix<T>& b, const Matrix<T>& delta, Matrix<T>& x, const T lambda, const T tol, const int itermax, const int numThreads = -1) const; 
+      void ridgeCG(const Matrix<T>& b, Matrix<T>& x, const T lambda, const T tol, const int itermax, const int numThreads = -1) const;
+      void ridgeCG(const Matrix<T>& b, const Matrix<T>& delta, Matrix<T>& x, const T lambda, const T tol, const int itermax, const int numThreads = -1) const;
 
       virtual ~AbstractMatrixB() { };
 };
@@ -247,7 +247,7 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
    /// Reference the column i to i+n into the Matrix mat
    inline void refSubMat(INTM i, INTM n, Matrix<T>& mat) const;
    /// extract a sub-matrix of a symmetric matrix
-   inline void subMatrixSym(const Vector<INTM>& indices, 
+   inline void subMatrixSym(const Vector<INTM>& indices,
          Matrix<T>& subMatrix) const;
    /// reference a modifiable reference to the data, DANGEROUS
    inline T* rawX() const { return _X; };
@@ -282,7 +282,7 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
    inline void set(const T a);
    /// Clear the matrix
    inline void clear();
-   /// Put white Gaussian noise in the matrix 
+   /// Put white Gaussian noise in the matrix
    inline void setAleat();
    /// set the matrix to the identity;
    inline void eye();
@@ -349,12 +349,12 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
    /// using two iterations of the power method
    inline void eigLargestSymApprox(const Vector<T>& u0,
          Vector<T>& u) const;
-   /// find the eigenvector corresponding to the eivenvalue with the 
+   /// find the eigenvector corresponding to the eivenvalue with the
    /// largest magnitude when the current matrix is symmetric,
-   /// using the power method. It 
-   /// returns the eigenvalue. u0 is an initial guess for the 
+   /// using the power method. It
+   /// returns the eigenvalue. u0 is an initial guess for the
    /// eigenvector.
-   inline T eigLargestMagnSym(const Vector<T>& u0, 
+   inline T eigLargestMagnSym(const Vector<T>& u0,
          Vector<T>& u) const;
    /// returns the value of the eigenvalue with the largest magnitude
    /// using the power iteration.
@@ -370,17 +370,17 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
    /// perform b = A'x, when x is sparse
    inline void multTrans(const SpVector<T>& x, Vector<T>& b, const T alpha =1.0, const T beta = 0.0) const;
    /// perform b = alpha*A*x+beta*b
-   inline void mult(const Vector<T>& x, Vector<T>& b, 
+   inline void mult(const Vector<T>& x, Vector<T>& b,
          const T alpha = 1.0, const T beta = 0.0) const;
    /// perform b = alpha*A*x + beta*b, when x is sparse
-   inline void mult(const SpVector<T>& x, Vector<T>& b, 
+   inline void mult(const SpVector<T>& x, Vector<T>& b,
          const T alpha = 1.0, const T beta = 0.0) const;
    /// perform C = a*A*B + b*C, possibly transposing A or B.
-   inline void mult(const Matrix<T>& B, Matrix<T>& C, 
+   inline void mult(const Matrix<T>& B, Matrix<T>& C,
          const bool transA = false, const bool transB = false,
          const T a = 1.0, const T b = 0.0) const;
    /// perform C = a*B*A + b*C, possibly transposing A or B.
-   inline void multSwitch(const Matrix<T>& B, Matrix<T>& C, 
+   inline void multSwitch(const Matrix<T>& B, Matrix<T>& C,
          const bool transA = false, const bool transB = false,
          const T a = 1.0, const T b = 0.0) const;
    /// perform C = A*B, when B is sparse
@@ -399,7 +399,7 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
    /// XXt = A*A'
    inline void XXt(Matrix<T>& XXt) const;
    /// XXt = A*A' where A is an upper triangular matrix
-   inline void upperTriXXt(Matrix<T>& XXt, 
+   inline void upperTriXXt(Matrix<T>& XXt,
          const INTM L) const;
    /// extract the diagonal
    inline void diag(Vector<T>& d) const;
@@ -508,7 +508,7 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
    inline void NadarayaWatson(const Vector<INTM>& ind, const T sigma);
    /// performs soft-thresholding of the vector
    inline void blockThrshold(const T nu, const INTM sizeGroup);
-   /// performs sparse projections of the columns 
+   /// performs sparse projections of the columns
    inline void sparseProject(Matrix<T>& out, const T thrs,   const int mode = 1, const T lambda1 = 0,
          const T lambda2 = 0, const T lambda3 = 0, const bool pos = false, const int numThreads=-1);
    inline void transformFilter();
@@ -518,7 +518,7 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
    inline void toSparse(SpMatrix<T>& matrix) const;
    /// make a sparse copy of the current matrix
    inline void toSparseTrans(SpMatrix<T>& matrixTrans);
-   /// make a reference of the matrix to a vector vec 
+   /// make a reference of the matrix to a vector vec
    inline void toVect(Vector<T>& vec) const;
    /// Accessor
    inline INTM V() const { return 1;};
@@ -532,7 +532,7 @@ template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, pu
 
    protected:
    /// Forbid lazy copies
-   explicit Matrix<T>(const Matrix<T>& matrix);
+   explicit Matrix(const Matrix<T>& matrix);
    /// Forbid lazy copies
    Matrix<T>& operator=(const Matrix<T>& matrix);
 
@@ -560,7 +560,7 @@ template<typename T> class Vector {
    /// Constructor with existing data
    Vector(T* X, INTM n);
    /// Copy constructor
-   explicit Vector<T>(const Vector<T>& vec);
+   explicit Vector(const Vector<T>& vec);
 
    /// Destructor
    virtual ~Vector();
@@ -611,7 +611,7 @@ template<typename T> class Vector {
    inline void refSubVec(INTM i, INTM n, Vector<T>& mat) const { mat.setData(_X+i,n); };
 
    /// put a random permutation of size n (for integral vectors)
-   inline void randperm(int n);  
+   inline void randperm(int n);
    /// put random values in the vector (White Gaussian Noise)
    inline void setAleat();
    /// clear the vector
@@ -661,15 +661,15 @@ template<typename T> class Vector {
    inline void div(const Vector<T>& x, const Vector<T>& y);
    /// A <- x .^ 2
    inline void sqr(const Vector<T>& x);
-   /// A <- 1 ./ sqrt(x) 
+   /// A <- 1 ./ sqrt(x)
    inline void sqr();
-   /// A <- 1 ./ sqrt(A) 
+   /// A <- 1 ./ sqrt(A)
    inline void Sqrt(const Vector<T>& x);
-   /// A <- 1 ./ sqrt(x) 
+   /// A <- 1 ./ sqrt(x)
    inline void Sqrt();
-   /// A <- 1 ./ sqrt(x) 
+   /// A <- 1 ./ sqrt(x)
    inline void Invsqrt(const Vector<T>& x);
-   /// A <- 1 ./ sqrt(A) 
+   /// A <- 1 ./ sqrt(A)
    inline void Invsqrt();
    /// A <- 1./x
    inline void inv(const Vector<T>& x);
@@ -751,7 +751,7 @@ template<typename T> class Vector {
 
 
    /// Conversion
-   /// make a sparse copy 
+   /// make a sparse copy
    inline void toSparse(SpVector<T>& vec) const;
    /// extract the rows of a matrix corresponding to a binary mask
    inline void copyMask(Vector<T>& out, Vector<bool>& mask) const;
@@ -759,7 +759,7 @@ template<typename T> class Vector {
 
 
    private:
-   /// = operator, 
+   /// = operator,
    Vector<T>& operator=(const Vector<T>& vec);
 
    /// if the data has been externally allocated
@@ -848,10 +848,10 @@ template<typename T> class SpMatrix : public Data<T>, public AbstractMatrixB<T>
    /// XAt <- X*A'
    inline void XAt(const Matrix<T>& X, Matrix<T>& XAt) const;
    /// XAt <- X(:,indices)*A(:,indices)'
-   inline void XAt(const Matrix<T>& X, Matrix<T>& XAt, 
+   inline void XAt(const Matrix<T>& X, Matrix<T>& XAt,
          const Vector<INTM>& indices) const;
    /// XAt <- sum_i w_i X(:,i)*A(:,i)'
-   inline void wXAt( const Vector<T>& w, const Matrix<T>& X, 
+   inline void wXAt( const Vector<T>& w, const Matrix<T>& X,
          Matrix<T>& XAt, const int numthreads=-1) const;
    inline void XtX(Matrix<T>& XtX) const;
 
@@ -861,17 +861,17 @@ template<typename T> class SpMatrix : public Data<T>, public AbstractMatrixB<T>
    inline void multTrans(const SpVector<T>& x, Vector<T>& y,
          const T alpha = 1.0, const T beta = 0.0) const;
    /// perform b = alpha*A*x + beta*b, when x is sparse
-   inline void mult(const SpVector<T>& x, Vector<T>& b, 
+   inline void mult(const SpVector<T>& x, Vector<T>& b,
          const T alpha = 1.0, const T beta = 0.0) const;
    /// perform b = alpha*A*x + beta*b, when x is sparse
-   inline void mult(const Vector<T>& x, Vector<T>& b, 
+   inline void mult(const Vector<T>& x, Vector<T>& b,
          const T alpha = 1.0, const T beta = 0.0) const;
    /// perform C = a*A*B + b*C, possibly transposing A or B.
-   inline void mult(const Matrix<T>& B, Matrix<T>& C, 
+   inline void mult(const Matrix<T>& B, Matrix<T>& C,
          const bool transA = false, const bool transB = false,
          const T a = 1.0, const T b = 0.0) const;
    /// perform C = a*B*A + b*C, possibly transposing A or B.
-   inline void multSwitch(const Matrix<T>& B, Matrix<T>& C, 
+   inline void multSwitch(const Matrix<T>& B, Matrix<T>& C,
          const bool transA = false, const bool transB = false,
          const T a = 1.0, const T b = 0.0) const;
    /// perform C = a*B*A + b*C, possibly transposing A or B.
@@ -919,7 +919,7 @@ template<typename T> class SpMatrix : public Data<T>, public AbstractMatrixB<T>
    bool _externAlloc;
    /// data
    T* _v;
-   /// row indices 
+   /// row indices
    INTM* _r;
    /// indices of the beginning of columns
    INTM* _pB;
@@ -975,9 +975,9 @@ template <typename T> class SpVector {
    inline T* rawX() const { return _v; };
    inline INTM* rawR() const { return _r; };
 
-   /// 
+   ///
    inline INTM L() const { return _L; };
-   /// 
+   ///
    inline void setL(const INTM L) { _L=L; };
    /// a <- a.^2
    inline void sqr();
@@ -1006,7 +1006,7 @@ template <typename T> class SpVector {
    explicit SpVector(const SpVector<T>& vector);
    SpVector<T>& operator=(const SpVector<T>& vector);
 
-   /// external allocation 
+   /// external allocation
    bool _externAlloc;
    /// data
    T* _v;
@@ -1036,7 +1036,7 @@ template<typename T> class ProdMatrix : public AbstractMatrix<T> {
 
       /// set_matrices
       inline void setMatrices(const Matrix<T>& D, const bool high_memory=true);
-      inline void setMatrices(const Matrix<T>& D, 
+      inline void setMatrices(const Matrix<T>& D,
             const Matrix<T>& X, const bool high_memory=true);
       /// compute DtX(:,i)
       inline void copyCol(const INTM i, Vector<T>& DtXi) const;
@@ -1070,7 +1070,7 @@ template<typename T> class ProdMatrix : public AbstractMatrix<T> {
 
 
 /* ************************************
- * Implementation of the class Matrix 
+ * Implementation of the class Matrix
  * ************************************/
 
 /// Constructor with existing data X of an m x n matrix
@@ -1153,7 +1153,7 @@ template <typename T> inline void Matrix<T>::getData(Vector<T>& x, const INTM i)
    this->copyCol(i,x);
 };
 
-template <typename T> inline void Matrix<T>::getGroup(Matrix<T>& data, 
+template <typename T> inline void Matrix<T>::getGroup(Matrix<T>& data,
       const vector_groups& groups, const INTM i) const {
    const group& gr = groups[i];
    const INTM N = gr.size();
@@ -1172,7 +1172,7 @@ template <typename T> inline void Matrix<T>::refCol(INTM i, Vector<T>& x) const
    x.clear();
    x._X=_X+i*_m;
    x._n=_m;
-   x._externAlloc=true; 
+   x._externAlloc=true;
 };
 
 /// Reference the column i to i+n INTMo the Matrix mat
@@ -1282,7 +1282,7 @@ template <typename T> inline void Matrix<T>::clear() {
    _externAlloc=true;
 };
 
-/// Put white Gaussian noise in the matrix 
+/// Put white Gaussian noise in the matrix
 template <typename T> inline void Matrix<T>::setAleat() {
    for (INTM i = 0; i<_n*_m; ++i) _X[i]=normalDistrib<T>();
 };
@@ -1307,7 +1307,7 @@ template <typename T> inline void Matrix<T>::normalize() {
          this->refCol(i,d);
          d.setAleat();
          d.normalize();
-      } 
+      }
    }
 };
 
@@ -1318,7 +1318,7 @@ template <typename T> inline void Matrix<T>::normalize2() {
       if (norm > 1.0) {
          T invNorm=1.0/norm;
          cblas_scal<T>(_m,invNorm,_X+_m*i,1);
-      } 
+      }
    }
 };
 
@@ -1336,12 +1336,12 @@ template <typename T> inline void Matrix<T>::center() {
 template <typename T> inline void Matrix<T>::center_rows() {
    Vector<T> mean_rows(_m);
    mean_rows.setZeros();
-   for (INTM i = 0; i<_n; ++i) 
-      for (INTM j = 0; j<_m; ++j) 
+   for (INTM i = 0; i<_n; ++i)
+      for (INTM j = 0; j<_m; ++j)
          mean_rows[j] += _X[i*_m+j];
    mean_rows.scal(T(1.0)/_n);
-   for (INTM i = 0; i<_n; ++i) 
-      for (INTM j = 0; j<_m; ++j) 
+   for (INTM i = 0; i<_n; ++i)
+      for (INTM j = 0; j<_m; ++j)
          _X[i*_m+j] -= mean_rows[j];
 };
 
@@ -1546,7 +1546,7 @@ template <typename T> inline void Matrix<T>::addToCols(
       const Vector<T>& cent) {
    Vector<T> col;
    for (INTM i = 0; i<_n; ++i) {
-      this->refCol(i,col);      
+      this->refCol(i,col);
       col.add(cent[i]);
    }
 };
@@ -1555,7 +1555,7 @@ template <typename T> inline void Matrix<T>::addVecToCols(
       const Vector<T>& vec, const T a) {
    Vector<T> col;
    for (INTM i = 0; i<_n; ++i) {
-      this->refCol(i,col);      
+      this->refCol(i,col);
       col.add(vec,a);
    }
 };
@@ -1593,7 +1593,7 @@ template <typename T> inline void Matrix<T>::singularValues(Vector<T>& u) const
       syev<T>(no,lower,_n,XtX.rawX(),_n,u.rawX());
       u.thrsPos();
       u.Sqrt();
-   } else if (_n > 10*_m) { 
+   } else if (_n > 10*_m) {
       Matrix<T> XXt;
       this->XXt(XXt);
       syev<T>(no,lower,_m,XXt.rawX(),_m,u.rawX());
@@ -1623,7 +1623,7 @@ template <typename T> inline void Matrix<T>::svd(Matrix<T>& U, Vector<T>& S, Mat
       Vt.transpose(V);
       Vector<T> inveigs;
       inveigs.copy(S);
-      for (INTM i = 0; i<num_eig; ++i) 
+      for (INTM i = 0; i<num_eig; ++i)
          if (S[i] > 1e-10) {
             inveigs[i]=T(1.0)/S[i];
          } else {
@@ -1638,7 +1638,7 @@ template <typename T> inline void Matrix<T>::svd(Matrix<T>& U, Vector<T>& S, Mat
       U.mult(*this,V,true,false);
       Vector<T> inveigs;
       inveigs.copy(S);
-      for (INTM i = 0; i<num_eig; ++i) 
+      for (INTM i = 0; i<num_eig; ++i)
          if (S[i] > 1e-10) {
             inveigs[i]=T(1.0)/S[i];
          } else {
@@ -1693,10 +1693,10 @@ template <typename T> inline void Matrix<T>::eigLargestSymApprox(
    }
 };
 
-/// find the eigenvector corresponding to the eivenvalue with the 
+/// find the eigenvector corresponding to the eivenvalue with the
 /// largest magnitude when the current matrix is symmetric,
-/// using the power method. It 
-/// returns the eigenvalue. u0 is an initial guess for the 
+/// using the power method. It
+/// returns the eigenvalue. u0 is an initial guess for the
 /// eigenvector.
 template <typename T> inline T Matrix<T>::eigLargestMagnSym(
       const Vector<T>& u0, Vector<T>& u) const {
@@ -1748,7 +1748,7 @@ template <typename T> inline void Matrix<T>::invSym() {
 };
 
 /// perform b = alpha*A'x + beta*b
-template <typename T> inline void Matrix<T>::multTrans(const Vector<T>& x, 
+template <typename T> inline void Matrix<T>::multTrans(const Vector<T>& x,
       Vector<T>& b, const T a, const T c) const {
    b.resize(_n);
    //   assert(x._n == _m && b._n == _n);
@@ -1756,7 +1756,7 @@ template <typename T> inline void Matrix<T>::multTrans(const Vector<T>& x,
 };
 
 /// perform b = A'x, when x is sparse
-template <typename T> inline void Matrix<T>::multTrans(const SpVector<T>& x, 
+template <typename T> inline void Matrix<T>::multTrans(const SpVector<T>& x,
       Vector<T>& b, const T alpha, const T beta) const {
    b.resize(_n);
    Vector<T> col;
@@ -1788,7 +1788,7 @@ template <typename T> inline void Matrix<T>::multTrans(
 };
 
 /// perform b = alpha*A*x+beta*b
-template <typename T> inline void Matrix<T>::mult(const Vector<T>& x, 
+template <typename T> inline void Matrix<T>::mult(const Vector<T>& x,
       Vector<T>& b, const T a, const T c) const {
    //  assert(x._n == _n && b._n == _m);
    b.resize(_m);
@@ -1796,7 +1796,7 @@ template <typename T> inline void Matrix<T>::mult(const Vector<T>& x,
 };
 
 /// perform b = alpha*A*x + beta*b, when x is sparse
-template <typename T> inline void Matrix<T>::mult(const SpVector<T>& x, 
+template <typename T> inline void Matrix<T>::mult(const SpVector<T>& x,
       Vector<T>& b, const T a, const T a2) const {
    if (!a2) {
       b.setZeros();
@@ -1815,7 +1815,7 @@ template <typename T> inline void Matrix<T>::mult(const SpVector<T>& x,
 };
 
 /// perform C = a*A*B + b*C, possibly transposing A or B.
-template <typename T> inline void Matrix<T>::mult(const Matrix<T>& B, 
+template <typename T> inline void Matrix<T>::mult(const Matrix<T>& B,
       Matrix<T>& C, const bool transA, const bool transB,
       const T a, const T b) const {
    CBLAS_TRANSPOSE trA,trB;
@@ -1831,11 +1831,11 @@ template <typename T> inline void Matrix<T>::mult(const Matrix<T>& B,
    }
    if (transB) {
       trB = CblasTrans;
-      n = B._m; 
+      n = B._m;
       //  assert(B._n == k);
    } else {
       trB = CblasNoTrans;
-      n = B._n; 
+      n = B._n;
       // assert(B._m == k);
    }
    C.resize(m,n);
@@ -1845,7 +1845,7 @@ template <typename T> inline void Matrix<T>::mult(const Matrix<T>& B,
 
 /// perform C = a*B*A + b*C, possibly transposing A or B.
 template <typename T>
-inline void Matrix<T>::multSwitch(const Matrix<T>& B, Matrix<T>& C, 
+inline void Matrix<T>::multSwitch(const Matrix<T>& B, Matrix<T>& C,
       const bool transA, const bool transB,
       const T a, const T b) const {
    B.mult(*this,C,transB,transA,a,b);
@@ -2098,7 +2098,7 @@ template <typename T> inline T Matrix<T>::asum() const {
 template <typename T> inline T Matrix<T>::trace() const {
    T sum=T();
    INTM m = MIN(_n,_m);
-   for (INTM i = 0; i<m; ++i) 
+   for (INTM i = 0; i<m; ++i)
       sum += _X[i*_m+i];
    return sum;
 };
@@ -2126,7 +2126,7 @@ template <typename T> inline T Matrix<T>::norm_inf_2_col() const {
    for (INTM i = 0; i<_n; ++i) {
       refCol(i,col);
       T norm_col = col.nrm2();
-      if (norm_col > max) 
+      if (norm_col > max)
          max = norm_col;
    }
    return max;
@@ -2148,10 +2148,10 @@ template <typename T> inline void Matrix<T>::norm_2_rows(
       Vector<T>& norms) const {
    norms.resize(_m);
    norms.setZeros();
-   for (INTM i = 0; i<_n; ++i) 
-      for (INTM j = 0; j<_m; ++j) 
+   for (INTM i = 0; i<_n; ++i)
+      for (INTM j = 0; j<_m; ++j)
          norms[j] += _X[i*_m+j]*_X[i*_m+j];
-   for (INTM j = 0; j<_m; ++j) 
+   for (INTM j = 0; j<_m; ++j)
       norms[j]=sqrt(norms[j]);
 };
 
@@ -2160,8 +2160,8 @@ template <typename T> inline void Matrix<T>::norm_2sq_rows(
       Vector<T>& norms) const {
    norms.resize(_m);
    norms.setZeros();
-   for (INTM i = 0; i<_n; ++i) 
-      for (INTM j = 0; j<_m; ++j) 
+   for (INTM i = 0; i<_n; ++i)
+      for (INTM j = 0; j<_m; ++j)
          norms[j] += _X[i*_m+j]*_X[i*_m+j];
 };
 
@@ -2192,8 +2192,8 @@ template <typename T> inline void Matrix<T>::norm_inf_cols(Vector<T>& norms) con
 template <typename T> inline void Matrix<T>::norm_inf_rows(Vector<T>& norms) const {
    norms.resize(_m);
    norms.setZeros();
-   for (INTM i = 0; i<_n; ++i) 
-      for (INTM j = 0; j<_m; ++j) 
+   for (INTM i = 0; i<_n; ++i)
+      for (INTM j = 0; j<_m; ++j)
          norms[j] = MAX(abs<T>(_X[i*_m+j]),norms[j]);
 };
 
@@ -2201,8 +2201,8 @@ template <typename T> inline void Matrix<T>::norm_inf_rows(Vector<T>& norms) con
 template <typename T> inline void Matrix<T>::norm_l1_rows(Vector<T>& norms) const {
    norms.resize(_m);
    norms.setZeros();
-   for (INTM i = 0; i<_n; ++i) 
-      for (INTM j = 0; j<_m; ++j) 
+   for (INTM i = 0; i<_n; ++i)
+      for (INTM j = 0; j<_m; ++j)
          norms[j] += abs<T>(_X[i*_m+j]);
 };
 
@@ -2219,7 +2219,7 @@ template <typename T> inline void Matrix<T>::norm_2sq_cols(
    }
 };
 
-template <typename T> 
+template <typename T>
 inline void Matrix<T>::sum_cols(Vector<T>& sum) const {
    sum.resize(_m);
    sum.setZeros();
@@ -2349,7 +2349,7 @@ template <typename T> inline void Matrix<T>::blockThrshold(const T nu,
    }
 }
 
-template <typename T> inline void Matrix<T>::sparseProject(Matrix<T>& Y, 
+template <typename T> inline void Matrix<T>::sparseProject(Matrix<T>& Y,
       const T thrs,   const int mode, const T lambda1,
       const T lambda2, const T lambda3, const bool pos,
       const int numThreads) {
@@ -2361,7 +2361,7 @@ template <typename T> inline void Matrix<T>::sparseProject(Matrix<T>& Y,
    }
 
    int i;
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
    for (i = 0; i< _n; ++i) {
 #ifdef _OPENMP
       int numT=omp_get_thread_num();
@@ -2417,7 +2417,7 @@ template <typename T> inline void Matrix<T>::rank1Update(
 };
 
 template <typename T>
-inline void Matrix<T>::rank1Update_mult(const Vector<T>& vec1, 
+inline void Matrix<T>::rank1Update_mult(const Vector<T>& vec1,
       const Vector<T>& vec1b,
       const SpVector<T>& vec2,
       const T alpha) {
@@ -2486,7 +2486,7 @@ template <typename T> inline void Matrix<T>::rank1Update(
 };
 
 
-/// compute x, such that b = Ax, 
+/// compute x, such that b = Ax,
 template <typename T> inline void Matrix<T>::conjugateGradient(
       const Vector<T>& b, Vector<T>& x, const T tol, const int itermax) const {
    Vector<T> R,P,AP;
@@ -2516,7 +2516,7 @@ template <typename T> inline void Matrix<T>::drop(char* fileName) const {
    f.open(fileName, ofstream::trunc);
    std::cout << "Matrix written in " << fileName << std::endl;
    for (INTM i = 0; i<_n; ++i) {
-      for (INTM j = 0; j<_m; ++j) 
+      for (INTM j = 0; j<_m; ++j)
          f << _X[i*_m+j] << " ";
       f << std::endl;
    }
@@ -2580,7 +2580,7 @@ template <typename T> inline void Matrix<T>::toSparse(SpMatrix<T>& out) const {
       pB=new INTM[_n+1];
    }
    INTM* pE=pB+1;
-   for (INTM i = 0; i<_n*_m; ++i) 
+   for (INTM i = 0; i<_n*_m; ++i)
       if (_X[i] != 0) ++count;
    INTM* r;
    T* v;
@@ -2621,7 +2621,7 @@ template <typename T> inline void Matrix<T>::toSparseTrans(
       pB=new INTM[_m+1];
    }
    INTM* pE=pB+1;
-   for (INTM i = 0; i<_n*_m; ++i) 
+   for (INTM i = 0; i<_n*_m; ++i)
       if (_X[i] != 0) ++count;
    INTM* r;
    T* v;
@@ -2651,7 +2651,7 @@ template <typename T> inline void Matrix<T>::toSparseTrans(
    out._externAlloc=false;
 };
 
-/// make a reference of the matrix to a vector vec 
+/// make a reference of the matrix to a vector vec
 template <typename T> inline void Matrix<T>::toVect(
       Vector<T>& vec) const {
    vec.clear();
@@ -2663,7 +2663,7 @@ template <typename T> inline void Matrix<T>::toVect(
 /// merge two dictionaries
 template <typename T> inline void Matrix<T>::merge(const Matrix<T>& B,
       Matrix<T>& C) const {
-   const INTM K =_n; 
+   const INTM K =_n;
    Matrix<T> G;
    this->mult(B,G,true,false);
    std::list<INTM> list;
@@ -2829,7 +2829,7 @@ inline void Vector<T>::logspace(const INTM n, const T a, const T b) {
 template <typename T>
 inline INTM Vector<T>::nnz() const {
    INTM sum=0;
-   for (INTM i = 0; i<_n; ++i) 
+   for (INTM i = 0; i<_n; ++i)
       if (_X[i] != T()) ++sum;
    return sum;
 };
@@ -2969,26 +2969,26 @@ template <typename T> inline void Vector<T>::hardThrshold(const T nu) {
 
 /// performs thresholding of the vector
 template <typename T> inline void Vector<T>::thrsmax(const T nu) {
-   for (INTM i = 0; i<_n; ++i) 
+   for (INTM i = 0; i<_n; ++i)
       _X[i]=MAX(_X[i],nu);
 }
 
 /// performs thresholding of the vector
 template <typename T> inline void Vector<T>::thrsmin(const T nu) {
-   for (INTM i = 0; i<_n; ++i) 
+   for (INTM i = 0; i<_n; ++i)
       _X[i]=MIN(_X[i],nu);
 }
 
 /// performs thresholding of the vector
 template <typename T> inline void Vector<T>::thrsabsmin(const T nu) {
-   for (INTM i = 0; i<_n; ++i) 
+   for (INTM i = 0; i<_n; ++i)
       _X[i]=MAX(MIN(_X[i],nu),-nu);
 }
 
 /// performs thresholding of the vector
 template <typename T> inline void Vector<T>::thrshold(const T nu) {
-   for (INTM i = 0; i<_n; ++i) 
-      if (abs<T>(_X[i]) < nu) 
+   for (INTM i = 0; i<_n; ++i)
+      if (abs<T>(_X[i]) < nu)
          _X[i]=0;
 }
 /// performs soft-thresholding of the vector
@@ -3368,7 +3368,7 @@ template <typename T> inline T Vector<T>::asum() const {
 
 template <typename T> inline T Vector<T>::lzero() const {
    INTM count=0;
-   for (INTM i = 0; i<_n; ++i) 
+   for (INTM i = 0; i<_n; ++i)
       if (_X[i] != 0) ++count;
    return count;
 };
@@ -3384,7 +3384,7 @@ template <typename T> inline T Vector<T>::afused() const {
 /// returns the sum of the vector
 template <typename T> inline T Vector<T>::sum() const {
    T sum=T();
-   for (INTM i = 0; i<_n; ++i) sum +=_X[i]; 
+   for (INTM i = 0; i<_n; ++i) sum +=_X[i];
    return sum;
 };
 
@@ -3393,7 +3393,7 @@ template <typename T> inline void Vector<T>::sign(Vector<T>& signs) const {
    T* prSign=signs.rawX();
    for (INTM i = 0; i<_n; ++i) {
       if (_X[i] == 0) {
-         prSign[i]=0.0; 
+         prSign[i]=0.0;
       } else {
          prSign[i] = _X[i] > 0 ? 1.0 : -1.0;
       }
@@ -3569,7 +3569,7 @@ inline void Vector<T>::sparseProject(Vector<T>& out, const T thrs, const int mod
    } else if (mode == 4) {
       /// min_u 0.5||b-u||_2^2  + lambda1||u||_1 / ||u||_2^2 <= thrs
       out.copy(*this);
-      if (pos) 
+      if (pos)
          out.thrsPos();
       out.softThrshold(lambda1);
       T nrm=out.nrm2sq();
@@ -3582,7 +3582,7 @@ inline void Vector<T>::sparseProject(Vector<T>& out, const T thrs, const int mod
       //      if (nrm > thrs)
       //         out.scal(sqr_alt<T>(thrs/nrm));
       //  } else if (mode == 6) {
-      /// min_u 0.5||b-u||_2^2  + lambda1||u||_1 +lambda2 Fused(u) +0.5lambda_3 ||u||_2^2 
+      /// min_u 0.5||b-u||_2^2  + lambda1||u||_1 +lambda2 Fused(u) +0.5lambda_3 ||u||_2^2
       this->fusedProjectHomotopy(out,lambda1,lambda2,lambda3,true);
 } else if (mode==6) {
    /// min_u ||b-u||_2^2  /  lambda1||u||_1 +lambda2 Fused(u) + 0.5lambda3||u||_2^2 <= thrs
@@ -3591,7 +3591,7 @@ inline void Vector<T>::sparseProject(Vector<T>& out, const T thrs, const int mod
    /// min_u ||b-u||_2^2 / (1-lambda1)*||u||_2^2 + lambda1||u||_1 <= thrs
    if (lambda1 < 1e-10) {
       out.copy(*this);
-      if (pos) 
+      if (pos)
          out.thrsPos();
       out.normalize2();
       out.scal(sqrt(thrs));
@@ -3619,7 +3619,7 @@ inline void Vector<T>::l1l2projectb(Vector<T>& out, const T thrs, const T gamma,
    } else if (mode == 3) {
       /// min_u 0.5||b-u||_2^2  + gamma||u||_1 / ||u||_2^2 <= thrs
       out.copy(*this);
-      if (pos) 
+      if (pos)
          out.thrsPos();
       out.softThrshold(gamma);
       T nrm=out.nrm2();
@@ -3632,7 +3632,7 @@ inline void Vector<T>::l1l2projectb(Vector<T>& out, const T thrs, const T gamma,
 /// min_u ||b-u||_2^2 / ||u||_1 + (gamma/2) ||u||_2^2 <= thrs
 template <typename T>
    inline void Vector<T>::l1l2project(Vector<T>& out, const T thrs, const T gamma, const bool pos) const {
-      if (gamma == 0) 
+      if (gamma == 0)
          return this->l1project(out,thrs,pos);
       out.copy(*this);
       if (pos) {
@@ -3713,7 +3713,7 @@ static inline T fusedHomotopyAux(const bool& sign1,
 };
 
 template <typename T>
-inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha, 
+inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
       const T lambda1,const T lambda2,const T lambda3,
       const bool penalty) {
    T* pr_DtR=_X;
@@ -3747,7 +3747,7 @@ inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
    alpha.set(pr_gamma[0]);
    /// update DtR
    this->sub(alpha);
-   for (INTM j = K-2; j>=0; --j) 
+   for (INTM j = K-2; j>=0; --j)
       pr_DtR[j] += pr_DtR[j+1];
 
    pr_DtR[0]=0;
@@ -3761,7 +3761,7 @@ inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
    /// Solve the Lasso using simplified LARS
    for (INTM i = 1; i<K; ++i) {
       /// exit if constraINTMs are satisfied
-      /// min_u ||b-u||_2^2  +  lambda1||u||_1 +lambda2 Fused(u) + 0.5lambda3||u||_2^2 
+      /// min_u ||b-u||_2^2  +  lambda1||u||_1 +lambda2 Fused(u) + 0.5lambda3||u||_2^2
       if (penalty && currentLambda <= lambda2) break;
       if (!penalty) {
          /// min_u ||b-u||_2^2  /  lambda1||u||_1 +lambda2 Fused(u) + 0.5lambda3||u||_2^2 <= 1.0
@@ -3775,7 +3775,7 @@ inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
       /// Update pr_ind and pr_c
       if (newAtom) {
          INTM j;
-         for (j = 1; j<i; ++j) 
+         for (j = 1; j<i; ++j)
             if (pr_ind[j] > currentInd) break;
          for (INTM k = i; k>j; --k) {
             pr_ind[k]=pr_ind[k-1];
@@ -3801,9 +3801,9 @@ inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
          }
          pr_u[i] = pr_signs[i-1] ? -pr_c[i-1] : pr_c[i-1];
          pr_u[i] += pr_signs[i] ? pr_c[i-1]+pr_c[i] : -pr_c[i-1]-pr_c[i];
-      } 
+      }
 
-      // Compute Du 
+      // Compute Du
       pr_Du[0]=pr_u[0];
       for (INTM k = 1; k<pr_ind[1]; ++k)
          pr_Du[k]=pr_Du[0];
@@ -3813,16 +3813,16 @@ inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
             pr_Du[k]=pr_Du[pr_ind[j]];
       }
 
-      /// Compute DDu 
+      /// Compute DDu
       DDu.copy(Du);
-      for (INTM j = K-2; j>=0; --j) 
+      for (INTM j = K-2; j>=0; --j)
          pr_DDu[j] += pr_DDu[j+1];
 
       /// Check constraINTMs
       T max_step1 = INFINITY;
       if (penalty) {
          max_step1 = currentLambda-lambda2;
-      } 
+      }
 
       /// Check changes of sign
       T max_step2 = INFINITY;
@@ -3849,7 +3849,7 @@ inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
       currentInd = scores.fmin();
       max_step3 = pr_scores[currentInd];
       T step = MIN(max_step1,MIN(max_step3,max_step2));
-      if (step == 0 || step == INFINITY) break; 
+      if (step == 0 || step == INFINITY) break;
 
       /// Update gamma, alpha, DtR, currentLambda
       for (INTM j = 0; j<=i; ++j) {
@@ -3860,10 +3860,10 @@ inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
       currentLambda -= step;
       if (step == max_step2) {
          /// Update signs,pr_ind, pr_c
-         for (INTM k = step_out; k<=i; ++k) 
+         for (INTM k = step_out; k<=i; ++k)
             pr_ind[k]=pr_ind[k+1];
          pr_ind[i]=K;
-         for (INTM k = step_out; k<=i; ++k) 
+         for (INTM k = step_out; k<=i; ++k)
             pr_signs[k]=pr_signs[k+1];
          pr_c[step_out-1]=T(1.0)/(pr_ind[step_out]-pr_ind[step_out-1]);
          pr_c[step_out]=T(1.0)/(pr_ind[step_out+1]-pr_ind[step_out]);
@@ -3892,7 +3892,7 @@ inline void Vector<T>::fusedProject(Vector<T>& alpha, const T lambda1, const T l
 
    T total_alpha =alpha.sum();
    /// Modification of beta
-   for (INTM i = K-2; i>=0; --i) 
+   for (INTM i = K-2; i>=0; --i)
       pr_beta[i]+=pr_beta[i+1];
 
    for (INTM i = 0; i<itermax; ++i) {
@@ -4062,7 +4062,7 @@ inline void Vector<T>::applyBayerPattern(const int offset) {
 };
 
 
-/// make a sparse copy 
+/// make a sparse copy
 template <typename T> inline void Vector<T>::toSparse(
       SpVector<T>& vec) const {
    INTM L=0;
@@ -4111,7 +4111,7 @@ inline void Matrix<T>::copyMask(Matrix<T>& out, Vector<bool>& mask) const {
 
 
 /* ****************************
- * Implementation of SpMatrix 
+ * Implementation of SpMatrix
  * ****************************/
 
 
@@ -4154,7 +4154,7 @@ template <typename T> SpMatrix<T>::~SpMatrix() {
 };
 
 /// reference the column i INTMo vec
-template <typename T> inline void SpMatrix<T>::refCol(INTM i, 
+template <typename T> inline void SpMatrix<T>::refCol(INTM i,
       SpVector<T>& vec) const {
    if (vec._nzmax > 0) vec.clear();
    vec._v=_v+_pB[i];
@@ -4192,7 +4192,7 @@ template<typename T>
 void SpMatrix<T>::getData(Vector<T>& data, const INTM index) const {
    data.resize(_m);
    data.setZeros();
-   for (INTM i = _pB[index]; i< _pB[index+1]; ++i) 
+   for (INTM i = _pB[index]; i< _pB[index+1]; ++i)
       data[_r[i]]=_v[i];
 };
 
@@ -4272,7 +4272,7 @@ template <typename T> inline void SpMatrix<T>::clear() {
 };
 
 /// resize the matrix
-template <typename T> inline void SpMatrix<T>::resize(const INTM m, 
+template <typename T> inline void SpMatrix<T>::resize(const INTM m,
       const INTM n, const INTM nzmax) {
    if (n == _n && m == _m && nzmax == _nzmax) return;
    this->clear();
@@ -4317,7 +4317,7 @@ inline void SpMatrix<T>::multTrans(const Vector<T>& x, Vector<T>& y,
 
 /// perform b = alpha*A*x + beta*b, when x is sparse
 template <typename T>
-inline void SpMatrix<T>::multTrans(const SpVector<T>& x, Vector<T>& y, 
+inline void SpMatrix<T>::multTrans(const SpVector<T>& x, Vector<T>& y,
       const T alpha, const T beta) const {
    y.resize(_n);
    if (beta) {
@@ -4356,7 +4356,7 @@ inline void SpMatrix<T>::mult(const Vector<T>& x, Vector<T>& y,
 
 /// perform b = alpha*A*x + beta*b, when x is sparse
 template <typename T>
-inline void SpMatrix<T>::mult(const SpVector<T>& x, Vector<T>& y, 
+inline void SpMatrix<T>::mult(const SpVector<T>& x, Vector<T>& y,
       const T alpha, const T beta) const {
    y.resize(_m);
    if (beta) {
@@ -4376,7 +4376,7 @@ inline void SpMatrix<T>::mult(const SpVector<T>& x, Vector<T>& y,
 
 /// perform C = a*A*B + b*C, possibly transposing A or B.
 template <typename T>
-inline void SpMatrix<T>::mult(const Matrix<T>& B, Matrix<T>& C, 
+inline void SpMatrix<T>::mult(const Matrix<T>& B, Matrix<T>& C,
       const bool transA, const bool transB,
       const T a, const T b) const {
    if (transA) {
@@ -4444,7 +4444,7 @@ inline void SpMatrix<T>::mult(const Matrix<T>& B, Matrix<T>& C,
 
 /// perform C = a*A*B + b*C, possibly transposing A or B.
 template <typename T>
-inline void SpMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C, 
+inline void SpMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
       const bool transA, const bool transB,
       const T a, const T b) const {
    if (transA) {
@@ -4512,7 +4512,7 @@ inline void SpMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
 
 /// perform C = a*B*A + b*C, possibly transposing A or B.
 template <typename T>
-inline void SpMatrix<T>::multSwitch(const Matrix<T>& B, Matrix<T>& C, 
+inline void SpMatrix<T>::multSwitch(const Matrix<T>& B, Matrix<T>& C,
       const bool transA, const bool transB,
       const T a, const T b) const {
    B.mult(*this,C,transB,transA,a,b);
@@ -4544,37 +4544,37 @@ inline void SpMatrix<T>::copyRow(const INTM ind, Vector<T>& x) const {
    }
 };
 
-template <typename T> 
+template <typename T>
 inline void SpMatrix<T>::addVecToCols(
       const Vector<T>& vec, const T a) {
    const T* pr_vec = vec.rawX();
    if (isEqual(a,T(1.0))) {
-      for (INTM i = 0; i<_n; ++i) 
-         for (INTM j = _pB[i]; j<_pE[i]; ++j) 
+      for (INTM i = 0; i<_n; ++i)
+         for (INTM j = _pB[i]; j<_pE[i]; ++j)
             _v[j] += pr_vec[_r[j]];
    } else {
-      for (INTM i = 0; i<_n; ++i) 
-         for (INTM j = _pB[i]; j<_pE[i]; ++j) 
+      for (INTM i = 0; i<_n; ++i)
+         for (INTM j = _pB[i]; j<_pE[i]; ++j)
             _v[j] += a*pr_vec[_r[j]];
    }
 };
 
-template <typename T> 
+template <typename T>
 inline void SpMatrix<T>::addVecToColsWeighted(
       const Vector<T>& vec, const T* weights, const T a) {
    const T* pr_vec = vec.rawX();
    if (isEqual(a,T(1.0))) {
-      for (INTM i = 0; i<_n; ++i) 
-         for (INTM j = _pB[i]; j<_pE[i]; ++j) 
+      for (INTM i = 0; i<_n; ++i)
+         for (INTM j = _pB[i]; j<_pE[i]; ++j)
             _v[j] += pr_vec[_r[j]]*weights[j-_pB[i]];
    } else {
-      for (INTM i = 0; i<_n; ++i) 
-         for (INTM j = _pB[i]; j<_pE[i]; ++j) 
+      for (INTM i = 0; i<_n; ++i)
+         for (INTM j = _pB[i]; j<_pE[i]; ++j)
             _v[j] += a*pr_vec[_r[j]]*weights[j-_pB[i]];
    }
 };
 
-template <typename T> 
+template <typename T>
 inline void SpMatrix<T>::sum_cols(Vector<T>& sum) const {
    sum.resize(_m);
    sum.setZeros();
@@ -4613,7 +4613,7 @@ template <typename T> inline void SpMatrix<T>::AAt(Matrix<T>& aat) const {
    }
 
    cblas_copy<T>(K*K,aatT,1,aat._X,1);
-   for (i = 1; i<NUM_THREADS; ++i) 
+   for (i = 1; i<NUM_THREADS; ++i)
       cblas_axpy<T>(K*K,1.0,aatT+K*K*i,1,aat._X,1);
    aat.fillSymmetric();
    delete[](aatT);
@@ -4663,7 +4663,7 @@ template <typename T> inline void SpMatrix<T>::AAt(Matrix<T>& aat,
    }
 
    cblas_copy<T>(K*K,aatT,1,aat._X,1);
-   for (i = 1; i<NUM_THREADS; ++i) 
+   for (i = 1; i<NUM_THREADS; ++i)
       cblas_axpy<T>(K*K,1.0,aatT+K*K*i,1,aat._X,1);
    aat.fillSymmetric();
    delete[](aatT);
@@ -4698,7 +4698,7 @@ template <typename T> inline void SpMatrix<T>::wAAt(const Vector<T>& w,
    }
 
    cblas_copy<T>(K*K,aatT,1,aat._X,1);
-   for (i = 1; i<NUM_THREADS; ++i) 
+   for (i = 1; i<NUM_THREADS; ++i)
       cblas_axpy<T>(K*K,1.0,aatT+K*K*i,1,aat._X,1);
    aat.fillSymmetric();
    delete[](aatT);
@@ -4732,7 +4732,7 @@ template <typename T> inline void SpMatrix<T>::XAt(const Matrix<T>& X,
    }
 
    cblas_copy<T>(n*K,XatT,1,XAt._X,1);
-   for (i = 1; i<NUM_THREADS; ++i) 
+   for (i = 1; i<NUM_THREADS; ++i)
       cblas_axpy<T>(n*K,1.0,XatT+n*K*i,1,XAt._X,1);
    delete[](XatT);
 };
@@ -4766,7 +4766,7 @@ template <typename T> inline void SpMatrix<T>::XAt(const Matrix<T>& X,
    }
 
    cblas_copy<T>(n*K,XatT,1,XAt._X,1);
-   for (i = 1; i<NUM_THREADS; ++i) 
+   for (i = 1; i<NUM_THREADS; ++i)
       cblas_axpy<T>(n*K,1.0,XatT+n*K*i,1,XAt._X,1);
    delete[](XatT);
 };
@@ -4805,7 +4805,7 @@ template <typename T> inline void SpMatrix<T>::wXAt(const Vector<T>& w,
    }
 
    cblas_copy<T>(n*K,XatT,1,XAt._X,1);
-   for (i = 1; i<NUM_THREADS; ++i) 
+   for (i = 1; i<NUM_THREADS; ++i)
       cblas_axpy<T>(n*K,1.0,XatT+n*K*i,1,XAt._X,1);
    delete[](XatT);
 };
@@ -4837,7 +4837,7 @@ template <typename T> inline void SpMatrix<T>::toFullTrans(
 
 
 /// use the data from v, r for _v, _r
-template <typename T> inline void SpMatrix<T>::convert(const Matrix<T>&vM, 
+template <typename T> inline void SpMatrix<T>::convert(const Matrix<T>&vM,
       const Matrix<INTM>& rM, const INTM K) {
    const INTM M = rM.n();
    const INTM L = rM.m();
@@ -4914,7 +4914,7 @@ inline void SpMatrix<T>::norm_1_cols(Vector<T>& norms) const {
 
 
 /* ***************************
- * Implementation of SpVector 
+ * Implementation of SpVector
  * ***************************/
 
 
@@ -4973,13 +4973,13 @@ template <typename T> inline void SpVector<T>::print(const string& name) const {
 /// create a reference on the vector r
 template <typename T> inline void SpVector<T>::refIndices(
       Vector<INTM>& indices) const {
-   indices.setPointer(_r,_L);   
+   indices.setPointer(_r,_L);
 };
 
 /// creates a reference on the vector val
 template <typename T> inline void SpVector<T>::refVal(
       Vector<T>& val) const {
-   val.setPointer(_v,_L);   
+   val.setPointer(_v,_L);
 };
 
 /// a <- a.^2
@@ -5071,7 +5071,7 @@ template <typename T> void inline SpVector<T>::toSpMatrix(
          out_r[i]=_r[i]-col*m;
       }
    }
-   for (current_col++ ; current_col < n+1; ++current_col) 
+   for (current_col++ ; current_col < n+1; ++current_col)
       out_pB[current_col]=_L;
 };
 
@@ -5088,9 +5088,9 @@ template <typename T> void inline SpVector<T>::toFull(Vector<T>& out)
  * ****************************/
 
 template <typename T> ProdMatrix<T>::ProdMatrix()  {
-   _DtX= NULL; 
-   _X=NULL; 
-   _D=NULL; 
+   _DtX= NULL;
+   _X=NULL;
+   _D=NULL;
    _high_memory=true;
    _n=0;
    _m=0;
@@ -5114,7 +5114,7 @@ template <typename T> ProdMatrix<T>::ProdMatrix(const Matrix<T>& D, const Matrix
 template <typename T> inline void ProdMatrix<T>::setMatrices(const Matrix<T>& D, const Matrix<T>& X,
       const bool high_memory)  {
    _high_memory=high_memory;
-   _m = D.n(); 
+   _m = D.n();
    _n = X.n();
    if (high_memory) {
       D.mult(X,*_DtX,true,false);
@@ -5129,7 +5129,7 @@ template <typename T> inline void ProdMatrix<T>::setMatrices(const Matrix<T>& D,
 template <typename T> inline void ProdMatrix<T>::setMatrices(
       const Matrix<T>& D, const bool high_memory) {
    _high_memory=high_memory;
-   _m = D.n(); 
+   _m = D.n();
    _n = D.n();
    if (high_memory) {
       D.XtX(*_DtX);
@@ -5137,7 +5137,7 @@ template <typename T> inline void ProdMatrix<T>::setMatrices(
       _X=&D;
       _D=&D;
       _DtX=NULL;
-   } 
+   }
    _addDiag=0;
 };
 
@@ -5150,7 +5150,7 @@ template <typename T> inline void ProdMatrix<T>::copyCol(const INTM i, Vector<T>
       _X->refCol(i,Xi);
       _D->multTrans(Xi,DtXi);
       if (_addDiag && _m == _n) DtXi[i] += _addDiag;
-   } 
+   }
 };
 
 /// compute DtX(:,i)
@@ -5163,7 +5163,7 @@ template <typename T> inline void ProdMatrix<T>::extract_rawCol(const INTM i,T*
       _X->refCol(i,Xi);
       _D->multTrans(Xi,vDtXi);
       if (_addDiag && _m == _n) DtXi[i] += _addDiag;
-   } 
+   }
 };
 
 template <typename T> inline void ProdMatrix<T>::add_rawCol(const INTM i,T* DtXi,
@@ -5176,7 +5176,7 @@ template <typename T> inline void ProdMatrix<T>::add_rawCol(const INTM i,T* DtXi
       _X->refCol(i,Xi);
       _D->multTrans(Xi,vDtXi,a,T(1.0));
       if (_addDiag && _m == _n) DtXi[i] += a*_addDiag;
-   } 
+   }
 };
 
 template <typename T> void inline ProdMatrix<T>::addDiag(const T diag) {
@@ -5254,7 +5254,7 @@ template <typename T> class SubMatrix : public AbstractMatrix<T> {
       AbstractMatrix<T>* _matrix;
 };
 
-template <typename T> 
+template <typename T>
 SubMatrix<T>::SubMatrix(AbstractMatrix<T>& G, Vector<INTM>& indI, Vector<INTM>& indJ) {
    _matrix = &G;
    _indicesI.copy(indI);
@@ -5287,7 +5287,7 @@ template <typename T> void inline SubMatrix<T>::extract_rawCol(const INTM i, T*
    }
 };
 
-template <typename T> inline void SubMatrix<T>::copyCol(const INTM i, 
+template <typename T> inline void SubMatrix<T>::copyCol(const INTM i,
       Vector<T>& DtXi) const {
    this->extract_rawCol(i,DtXi.rawX());
 };
@@ -5309,7 +5309,7 @@ template <typename T> void inline SubMatrix<T>::diag(Vector<T>& diag) const {
    }
 };
 
-template <typename T> inline T SubMatrix<T>::operator()(const INTM index1, 
+template <typename T> inline T SubMatrix<T>::operator()(const INTM index1,
       const INTM index2) const {
    return (*_matrix)(_indicesI[index1],_indicesJ[index2]);
 }
@@ -5330,23 +5330,23 @@ template <typename T> class ShiftMatrix : public AbstractMatrixB<T> {
             const T alpha = 1.0, const T beta = 0.0) const;
 
       /// perform b = alpha*A*x + beta*b, when x is sparse
-      virtual void mult(const SpVector<T>& x, Vector<T>& b, 
+      virtual void mult(const SpVector<T>& x, Vector<T>& b,
             const T alpha = 1.0, const T beta = 0.0) const;
 
-      virtual void mult(const Vector<T>& x, Vector<T>& b, 
+      virtual void mult(const Vector<T>& x, Vector<T>& b,
             const T alpha = 1.0, const T beta = 0.0) const;
 
       /// perform C = a*A*B + b*C, possibly transposing A or B.
-      virtual void mult(const Matrix<T>& B, Matrix<T>& C, 
+      virtual void mult(const Matrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const;
 
-      virtual void mult(const SpMatrix<T>& B, Matrix<T>& C, 
+      virtual void mult(const SpMatrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const;
 
       /// perform C = a*B*A + b*C, possibly transposing A or B.
-      virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C, 
+      virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const;
 
@@ -5363,7 +5363,7 @@ template <typename T> class ShiftMatrix : public AbstractMatrixB<T> {
       virtual ~ShiftMatrix() {  };
 
    private:
-      void center() { 
+      void center() {
          Vector<T> ones(_m);
          ones.set(T(1.0)/_m);
          this->multTrans(ones,_means);
@@ -5457,7 +5457,7 @@ template <typename T> void ShiftMatrix<T>::mult(const Matrix<T>&
    cerr << "Shift Matrix is used in inadequate setting" << endl;
 }
 
-template <typename T> void ShiftMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C, 
+template <typename T> void ShiftMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
       const bool transA, const bool transB, const T a, const T b) const {
    cerr << "Shift Matrix is used in inadequate setting" << endl;
 }
@@ -5502,7 +5502,7 @@ template <typename T> void ShiftMatrix<T>::print(const string& name) const {
 /// Matrix with shifts
 template <typename T> class DoubleRowMatrix : public AbstractMatrixB<T> {
    public:
-      DoubleRowMatrix(const AbstractMatrixB<T>& inputmatrix) : _inputmatrix(&inputmatrix) { 
+      DoubleRowMatrix(const AbstractMatrixB<T>& inputmatrix) : _inputmatrix(&inputmatrix) {
          _n=inputmatrix.n();
          _m=2*inputmatrix.m();
       };
@@ -5514,23 +5514,23 @@ template <typename T> class DoubleRowMatrix : public AbstractMatrixB<T> {
             const T alpha = 1.0, const T beta = 0.0) const;
 
       /// perform b = alpha*A*x + beta*b, when x is sparse
-      virtual void mult(const SpVector<T>& x, Vector<T>& b, 
+      virtual void mult(const SpVector<T>& x, Vector<T>& b,
             const T alpha = 1.0, const T beta = 0.0) const;
 
-      virtual void mult(const Vector<T>& x, Vector<T>& b, 
+      virtual void mult(const Vector<T>& x, Vector<T>& b,
             const T alpha = 1.0, const T beta = 0.0) const;
 
       /// perform C = a*A*B + b*C, possibly transposing A or B.
-      virtual void mult(const Matrix<T>& B, Matrix<T>& C, 
+      virtual void mult(const Matrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const;
 
-      virtual void mult(const SpMatrix<T>& B, Matrix<T>& C, 
+      virtual void mult(const SpMatrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const;
 
       /// perform C = a*B*A + b*C, possibly transposing A or B.
-      virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C, 
+      virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C,
             const bool transA = false, const bool transB = false,
             const T a = 1.0, const T b = 0.0) const;
 
@@ -5557,7 +5557,7 @@ template <typename T> void DoubleRowMatrix<T>::multTrans(const
       Vector<T>& x, Vector<T>& b, const T alpha, const T beta) const {
    const INTM mm = _inputmatrix->m();
    Vector<T> tmp(mm);
-   for (INTM i = 0; i<mm; ++i) 
+   for (INTM i = 0; i<mm; ++i)
       tmp[i]=x[2*i]+x[2*i+1];
    _inputmatrix->multTrans(tmp,b,alpha,beta);
 };
@@ -5607,7 +5607,7 @@ template <typename T> void DoubleRowMatrix<T>::mult(const Matrix<T>&
    cerr << "Double Matrix is used in inadequate setting" << endl;
 }
 
-template <typename T> void DoubleRowMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C, 
+template <typename T> void DoubleRowMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
       const bool transA, const bool transB, const T a, const T b) const {
    FLAG(4)
    cerr << "Double Matrix is used in inadequate setting" << endl;
@@ -5778,7 +5778,7 @@ template <typename T> void AbstractMatrixB<T>::ridgeCG(const Matrix<T>& mb,
    init_omp(numThreads);
    const int num_probs=mb.n();
    int i = 0;
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
    for (i =0; i<num_probs; ++i) {
       Vector<T> b;
       mb.refCol(i,b);
@@ -5794,7 +5794,7 @@ template <typename T> void AbstractMatrixB<T>::ridgeCG(const Matrix<T>& mb, cons
    init_omp(numThreads);
    const int num_probs=mb.n();
    int i = 0;
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
    for (i =0; i<num_probs; ++i) {
       Vector<T> b;
       mb.refCol(i,b);
diff --git a/spams_wrap/spams/prox/fista.h b/spams_wrap/spams/prox/fista.h
index 22bd08c4..38e8b4c9 100644
--- a/spams_wrap/spams/prox/fista.h
+++ b/spams_wrap/spams/prox/fista.h
@@ -1,5 +1,5 @@
 
-/* Software SPAMS v2.1 - Copyright 2009-2011 Julien Mairal 
+/* Software SPAMS v2.1 - Copyright 2009-2011 Julien Mairal
  *
  * This file is part of SPAMS.
  *
@@ -135,18 +135,18 @@ namespace FISTA {
    };
 
    bool regul_for_matrices(const regul_t& regul) {
-      return regul==L1L2 || regul==L1LINF || regul==L1L2_L1 || regul==L1LINF_L1 
+      return regul==L1L2 || regul==L1LINF || regul==L1L2_L1 || regul==L1LINF_L1
          || regul==TREEMULT || regul==GRAPHMULT || regul==L1LINFCR ||
          regul==TRACE_NORM || regul==RANK;
    }
 
-   template <typename T> struct ParamFISTA { 
-      ParamFISTA() { num_threads=1; max_it=100; L0=T(0.1); gamma=T(1.5); tol=T(1e-10); 
+   template <typename T> struct ParamFISTA {
+      ParamFISTA() { num_threads=1; max_it=100; L0=T(0.1); gamma=T(1.5); tol=T(1e-10);
          it0=10; max_iter_backtracking=1000; loss=SQUARE; compute_gram=false; admm=false; lin_admm=false;
-         intercept=false; regul=RIDGE; resetflow=false; delta=0; lambda2=0; lambda3=0; verbose=false; 
+         intercept=false; regul=RIDGE; resetflow=false; delta=0; lambda2=0; lambda3=0; verbose=false;
          pos=false; clever=true; a=T(1.0); b=T(0.0); c=T(1.0);
          log=false; logName=NULL; ista=false; subgrad=false;
-         length_names=30; 
+         length_names=30;
          name_regul=new char[length_names];
          name_loss=new char[length_names];
          is_inner_weights=false;
@@ -162,10 +162,10 @@ namespace FISTA {
          ngroups=0;
          linesearch_mode=0;
       }
-      ~ParamFISTA() { 
+      ~ParamFISTA() {
          if (!copied) {
-            delete[](name_regul); 
-            delete[](name_loss); 
+            delete[](name_regul);
+            delete[](name_loss);
          }
       };
       int num_threads;
@@ -213,7 +213,7 @@ namespace FISTA {
       int linesearch_mode;
    };
 
-   template <typename T> struct ParamReg { 
+   template <typename T> struct ParamReg {
       ParamReg() { size_group=1; lambda2d1 = 0; lambda=0; lambda3d1 = 0; pos=false; intercept=false; num_cols=1; graph_st=NULL; tree_st=NULL;
       graph_path_st=NULL; resetflow=false; clever=false; linf=true; transpose=false; ngroups=0;
       groups=NULL; };
@@ -237,12 +237,12 @@ namespace FISTA {
 
    template <typename T>
       bool param_for_admm(const ParamFISTA<T>& param) {
-         return (param.admm) && (param.loss==SQUARE || param.loss == HINGE) 
+         return (param.admm) && (param.loss==SQUARE || param.loss == HINGE)
             && (param.regul==GRAPH_L2 || param.regul==GRAPH || param.regul == NONE);
       };
 
 
-   template <typename T,  typename F = Matrix<T>, typename D = Vector<T> , 
+   template <typename T,  typename F = Matrix<T>, typename D = Vector<T> ,
             typename E = Vector<T> >
                class SplittingFunction {
                   public:
@@ -263,7 +263,7 @@ namespace FISTA {
                      virtual void add_mult_design_matrix(const E& prim, E& out, const T fact) const { };
 
                   private:
-                     explicit SplittingFunction<T,F,D,E>(const SplittingFunction<T,F,D,E>& loss);
+                     explicit SplittingFunction(const SplittingFunction<T,F,D,E>& loss);
                      SplittingFunction<T,F,D,E>& operator=(const SplittingFunction<T,F,D,E>& loss);
                };
 
@@ -288,18 +288,18 @@ namespace FISTA {
                   const bool intercept = false) const = 0;
 
          private:
-            explicit Loss<T,D,E>(const Loss<T,D,E>& dict);
+            explicit Loss(const Loss<T,D,E>& dict);
             Loss<T,D,E>& operator=(const Loss<T,D,E>& dict);
       };
 
-   template <typename T> 
+   template <typename T>
       class SqLossMissing : public Loss<T> {
          public:
             SqLossMissing(const AbstractMatrixB<T>& D) : _D(&D) { };
             virtual ~SqLossMissing() { };
 
-            inline void init(const Vector<T>& x) { 
-               _x.copy(x); 
+            inline void init(const Vector<T>& x) {
+               _x.copy(x);
                _missingvalues.clear();
                for (int i = 0; i<_x.n(); ++i) {
                   if (isnan(_x[i])) {
@@ -335,8 +335,8 @@ namespace FISTA {
             virtual T fenchel(const Vector<T>& input) const {
                return 0.5*input.nrm2sq()+input.dot(_x);
             };
-            virtual void var_fenchel(const Vector<T>& x, 
-                  Vector<T>& grad1, Vector<T>& grad2, 
+            virtual void var_fenchel(const Vector<T>& x,
+                  Vector<T>& grad1, Vector<T>& grad2,
                   const bool intercept) const {
                grad1.copy(_x);
                SpVector<T> spalpha(x.n());
@@ -351,25 +351,25 @@ namespace FISTA {
             };
 
          private:
-            explicit SqLossMissing<T>(const SqLossMissing<T>& dict);
+            explicit SqLossMissing(const SqLossMissing<T>& dict);
             SqLossMissing<T>& operator=(const SqLossMissing<T>& dict);
             const AbstractMatrixB<T>* _D;
             Vector<T> _x;
             List<int> _missingvalues;
       };
 
-   template <typename T> 
+   template <typename T>
       class SqLoss : public Loss<T>, public SplittingFunction<T> {
          public:
             SqLoss(const AbstractMatrixB<T>& D) : _D(&D) { _compute_gram = false; };
             SqLoss(const AbstractMatrixB<T>& D, const Matrix<T>& G) : _D(&D), _G(&G) { _compute_gram = true; };
             virtual ~SqLoss() { };
 
-            inline void init(const Vector<T>& x) { 
-               _x.copy(x); 
+            inline void init(const Vector<T>& x) {
+               _x.copy(x);
                if (_compute_gram) {
                   _D->multTrans(x,_DtX);
-               } 
+               }
             };
 
             inline T eval(const Vector<T>& alpha) const {
@@ -414,8 +414,8 @@ namespace FISTA {
             virtual T fenchel(const Vector<T>& input) const {
                return 0.5*input.nrm2sq()+input.dot(_x);
             };
-            virtual void var_fenchel(const Vector<T>& x, 
-                  Vector<T>& grad1, Vector<T>& grad2, 
+            virtual void var_fenchel(const Vector<T>& x,
+                  Vector<T>& grad1, Vector<T>& grad2,
                   const bool intercept) const {
                grad1.copy(_x);
                SpVector<T> spalpha(x.n());
@@ -460,7 +460,7 @@ namespace FISTA {
                prim_var.resize(_D->m());
                prim_var.setZeros();
             }
-            virtual void prox_prim_var(Vector<T>& out,const Vector<T>& dual_var, 
+            virtual void prox_prim_var(Vector<T>& out,const Vector<T>& dual_var,
                   const Vector<T>& prim_var, const T c) const {
                const T gamma=T(1.0)/c;
                out.copy(dual_var);
@@ -469,8 +469,8 @@ namespace FISTA {
                out.add(_x,gamma);
                out.scal(T(1.0)/(T(1.0)+gamma));
             };
-            inline void compute_new_prim(Vector<T>& prim, const Vector<T>& prim_var, 
-                  const Vector<T>& dual_var, const T gamma, const T delta) const { 
+            inline void compute_new_prim(Vector<T>& prim, const Vector<T>& prim_var,
+                  const Vector<T>& dual_var, const T gamma, const T delta) const {
                Vector<T> tmp;
                _D->mult(prim,tmp);
                tmp.scal(-gamma);
@@ -478,13 +478,13 @@ namespace FISTA {
                tmp.add(dual_var,gamma);
                _D->multTrans(tmp,prim,T(1.0),delta);
             };
-            inline void add_mult_design_matrix(const Vector<T>& prim, 
-                  Vector<T>& out, const T fact) const { 
+            inline void add_mult_design_matrix(const Vector<T>& prim,
+                  Vector<T>& out, const T fact) const {
                _D->mult(prim,out,fact,T(1.0));
             };
 
          private:
-            explicit SqLoss<T>(const SqLoss<T>& dict);
+            explicit SqLoss(const SqLoss<T>& dict);
             SqLoss<T>& operator=(const SqLoss<T>& dict);
             const AbstractMatrixB<T>* _D;
             Vector<T> _x;
@@ -493,13 +493,13 @@ namespace FISTA {
             Vector<T> _DtX;
       };
 
-   template <typename T> 
+   template <typename T>
       class HingeLoss : public SplittingFunction<T > {
          public:
             HingeLoss(const AbstractMatrixB<T>& X) : _X(&X) {  };
             virtual ~HingeLoss() { };
 
-            inline void init(const Vector<T>& y) { 
+            inline void init(const Vector<T>& y) {
                _y.copy(y);
             };
             inline T eval(const Vector<T>& w) const {
@@ -533,7 +533,7 @@ namespace FISTA {
                prim_var.resize(_X->m());
                prim_var.setZeros();
             }
-/*            inline void prox_prim_var(Vector<T>& out,const Vector<T>& dual_var, 
+/*            inline void prox_prim_var(Vector<T>& out,const Vector<T>& dual_var,
                   const Vector<T>& prim_var, const T lambda, const T c) const {
                const T gamma=T(1.0)/c;
                out.copy(dual_var);
@@ -549,8 +549,8 @@ namespace FISTA {
                   }
                }
             }*/
-            inline void compute_new_prim(Vector<T>& prim, const Vector<T>& prim_var, 
-                  const Vector<T>& dual_var, const T gamma, const T delta) const { 
+            inline void compute_new_prim(Vector<T>& prim, const Vector<T>& prim_var,
+                  const Vector<T>& dual_var, const T gamma, const T delta) const {
                Vector<T> tmp;
                _X->mult(prim,tmp);
                tmp.scal(-gamma);
@@ -559,7 +559,7 @@ namespace FISTA {
                _X->multTrans(tmp,prim,T(1.0),delta);
             };
             inline void add_mult_design_matrix(const Vector<T>& prim, Vector<T>& out,
-                  const T fact) const { 
+                  const T fact) const {
                _X->mult(prim,out,fact,T(1.0));
             };
             inline void prox_split(Matrix<T>& splitted_w, const T lambda) const {
@@ -581,25 +581,25 @@ namespace FISTA {
             };
 
          private:
-            explicit HingeLoss<T>(const HingeLoss<T>& dict);
+            explicit HingeLoss(const HingeLoss<T>& dict);
             HingeLoss<T>& operator=(const HingeLoss<T>& dict);
 
             const AbstractMatrixB<T>* _X;
             Vector<T> _y;
       };
 
-   template <typename T, bool weighted = false> 
+   template <typename T, bool weighted = false>
       class LogLoss : public Loss<T> {
          public:
             LogLoss(const AbstractMatrixB<T>& X) : _X(&X) {  };
             virtual ~LogLoss() { };
 
-            inline void init(const Vector<T>& y) { 
+            inline void init(const Vector<T>& y) {
                _y.copy(y);
                if (weighted) {
                   int countpos=0;
                   for (int i = 0; i<y.n(); ++i)
-                     if (y[i]>0) countpos++; 
+                     if (y[i]>0) countpos++;
                   _weightpos=T(1.0)/countpos;
                   _weightneg=T(1.0)/MAX(1e-3,(y.n()-countpos));
                }
@@ -677,7 +677,7 @@ namespace FISTA {
                _X->multTrans(grad1,grad2);
             };
          private:
-            explicit LogLoss<T,weighted>(const LogLoss<T,weighted>& dict);
+            explicit LogLoss(const LogLoss<T,weighted>& dict);
             LogLoss<T,weighted>& operator=(const LogLoss<T,weighted>& dict);
 
             const AbstractMatrixB<T>* _X;
@@ -686,14 +686,14 @@ namespace FISTA {
             T _weightneg;
       };
 
-   template <typename T> 
+   template <typename T>
       class MultiLogLoss : public Loss<T, Matrix<T> > {
          public:
             MultiLogLoss(const AbstractMatrixB<T>& X) : _X(&X) {  };
 
             virtual ~MultiLogLoss() { };
 
-            inline void init(const Vector<T>& y) { 
+            inline void init(const Vector<T>& y) {
                _y.resize(y.n());
                for (int i = 0; i<y.n(); ++i)
                   _y[i] = static_cast<int>(y[i]);
@@ -788,20 +788,20 @@ namespace FISTA {
                _X->mult(grad1,grad2,true,true);
             };
          private:
-            explicit MultiLogLoss<T>(const MultiLogLoss<T>& dict);
+            explicit MultiLogLoss(const MultiLogLoss<T>& dict);
             MultiLogLoss<T>& operator=(const MultiLogLoss<T>& dict);
 
             const AbstractMatrixB<T>* _X;
             Vector<int> _y;
       };
 
-   template <typename T> 
+   template <typename T>
       class PoissonLoss : public Loss<T> {
          public:
             PoissonLoss(const AbstractMatrixB<T>& X, const T delta) : _X(&X), _delta(delta) {  };
             virtual ~PoissonLoss() { };
 
-            inline void init(const Vector<T>& y) { 
+            inline void init(const Vector<T>& y) {
                _y.copy(y);
             };
             inline T eval(const Vector<T>& w) const {
@@ -845,7 +845,7 @@ namespace FISTA {
                }
                return sum;
             };
-            virtual void var_fenchel(const Vector<T>& w, Vector<T>& grad1, 
+            virtual void var_fenchel(const Vector<T>& w, Vector<T>& grad1,
             Vector<T>& grad2, const bool intercept) const {
                grad1.resize(_X->m());
                SpVector<T> spw(w.n());
@@ -859,7 +859,7 @@ namespace FISTA {
                _X->multTrans(grad1,grad2);
             };
          private:
-            explicit PoissonLoss<T>(const PoissonLoss<T>& dict);
+            explicit PoissonLoss(const PoissonLoss<T>& dict);
             PoissonLoss<T>& operator=(const PoissonLoss<T>& dict);
 
             const AbstractMatrixB<T>* _X;
@@ -867,7 +867,7 @@ namespace FISTA {
             T _delta;
       };
 
-   template <typename T> 
+   template <typename T>
       class LossCur: public Loss<T, Matrix<T>, Matrix<T> > {
          public:
             LossCur(const AbstractMatrixB<T>& X) : _X(&X) {  };
@@ -915,25 +915,25 @@ namespace FISTA {
                _X->mult(tmp,grad2,true,false);
             };
          private:
-            explicit LossCur<T>(const LossCur<T>& dict);
+            explicit LossCur(const LossCur<T>& dict);
             LossCur<T>& operator=(const LossCur<T>& dict);
 
             const AbstractMatrixB<T>* _X;
       };
 
-   template <typename T> 
+   template <typename T>
       class SqLossMat : public Loss<T, Matrix<T> , Matrix<T> > {
          public:
             SqLossMat(const AbstractMatrixB<T>& D) : _D(&D) { _compute_gram = false; };
-            SqLossMat(const AbstractMatrixB<T>& D, const Matrix<T>& G) : _D(&D), _G(&G) { 
+            SqLossMat(const AbstractMatrixB<T>& D, const Matrix<T>& G) : _D(&D), _G(&G) {
                _compute_gram = true; };
             virtual ~SqLossMat() { };
 
-            virtual inline void init(const Matrix<T>& x) { 
-               _x.copy(x); 
+            virtual inline void init(const Matrix<T>& x) {
+               _x.copy(x);
                if (_compute_gram) {
                   _D->mult(x,_DtX,true,false);
-               } 
+               }
             };
 
             inline T eval(const Matrix<T>& alpha) const {
@@ -991,7 +991,7 @@ namespace FISTA {
             };
 
          private:
-            explicit SqLossMat<T>(const SqLossMat<T>& dict);
+            explicit SqLossMat(const SqLossMat<T>& dict);
             SqLossMat<T>& operator=(const SqLossMat<T>& dict);
             const AbstractMatrixB<T>* _D;
             Matrix<T> _x;
@@ -1005,7 +1005,7 @@ namespace FISTA {
          public:
             LossMatSup() { };
 
-            virtual ~LossMatSup() { 
+            virtual ~LossMatSup() {
                for (int i = 0; i<_N; ++i) {
                   delete(_losses[i]);
                   _losses[i]=NULL;
@@ -1062,14 +1062,14 @@ namespace FISTA {
             };
             virtual bool is_fenchel() const {
                bool ok=true;
-               for (int i = 0; i<_N; ++i) 
+               for (int i = 0; i<_N; ++i)
                   ok = ok && _losses[i]->is_fenchel();
                return ok;
             };
             virtual void dummy() = 0;
 
          private:
-            explicit LossMatSup<T,L>(const LossMatSup<T,L>& dict);
+            explicit LossMatSup(const LossMatSup<T,L>& dict);
             LossMatSup<T,L>& operator=(const LossMatSup<T,L>& dict);
             int _m;
 
@@ -1088,7 +1088,7 @@ namespace FISTA {
                this->_N=N;
                this->_losses=new LogLoss<T,weighted>*[this->_N];
                Vector<T> col;
-               for (int i = 0; i<this->_N; ++i) 
+               for (int i = 0; i<this->_N; ++i)
                   this->_losses[i]=new LogLoss<T,weighted>(X);
             }
             virtual void dummy() { };
@@ -1102,7 +1102,7 @@ namespace FISTA {
                this->_N=N;
                this->_losses=new SqLossMissing<T>*[this->_N];
                Vector<T> col;
-               for (int i = 0; i<this->_N; ++i) 
+               for (int i = 0; i<this->_N; ++i)
                   this->_losses[i]=new SqLossMissing<T>(X);
             }
             virtual void dummy() { };
@@ -1127,7 +1127,7 @@ namespace FISTA {
       class Regularizer {
          public:
             Regularizer() { };
-            Regularizer(const ParamReg<T>& param) : _id(NA) { 
+            Regularizer(const ParamReg<T>& param) : _id(NA) {
                _intercept=param.intercept;
                _pos=param.pos;
             }
@@ -1152,7 +1152,7 @@ namespace FISTA {
             virtual bool is_concave() const { return false; };
 //            virtual bool is_none() const { return false; };
 //            virtual bool is_pos() const { return _pos; };
-            
+
 
          protected:
             bool _pos;
@@ -1160,11 +1160,11 @@ namespace FISTA {
             regul_t _id;
 
          private:
-            explicit Regularizer<T,D>(const Regularizer<T,D>& reg);
+            explicit Regularizer(const Regularizer<T,D>& reg);
             Regularizer<T,D>& operator=(const Regularizer<T,D>& reg);
       };
 
-   template <typename T> 
+   template <typename T>
       class Lasso : public Regularizer<T> {
          public:
             Lasso(const ParamReg<T>& param) : Regularizer<T>(param) { this->_id = L1; };
@@ -1176,7 +1176,7 @@ namespace FISTA {
                y.softThrshold(lambda);
                if (this->_intercept) y[y.n()-1] = x[y.n()-1];
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                return (this->_intercept ? x.asum() - abs(x[x.n()-1]) : x.asum());
             };
             void inline fenchel(const Vector<T>& input, T& val, T& scal) const {
@@ -1186,10 +1186,10 @@ namespace FISTA {
                T mm = output.fmaxval();
                scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                val=0;
-               if (this->_intercept & (abs<T>(output[output.n()-1]) > EPSILON)) val=INFINITY; 
+               if (this->_intercept & (abs<T>(output[output.n()-1]) > EPSILON)) val=INFINITY;
             };
             virtual bool is_subgrad() const { return true; };
-            virtual void sub_grad(const Vector<T>& input, Vector<T>& output) const {  
+            virtual void sub_grad(const Vector<T>& input, Vector<T>& output) const {
                output.resize(input.n());
                if (!this->_pos) {
                   for (int i = 0; i<input.n(); ++i) {
@@ -1204,12 +1204,12 @@ namespace FISTA {
             }
       };
 
-   template <typename T> 
+   template <typename T>
       class LassoConstraint : public Regularizer<T> {
          public:
-            LassoConstraint(const ParamReg<T>& param) : Regularizer<T>(param) { 
+            LassoConstraint(const ParamReg<T>& param) : Regularizer<T>(param) {
                _thrs=param.lambda;
-               this->_id = L1CONSTRAINT; 
+               this->_id = L1CONSTRAINT;
             };
             virtual ~LassoConstraint() { };
 
@@ -1227,7 +1227,7 @@ namespace FISTA {
             T inline eval(const Vector<T>& x) const {
                return 0;
             };
-            void inline fenchel(const Vector<T>& input, T& val, T& scal) const { 
+            void inline fenchel(const Vector<T>& input, T& val, T& scal) const {
                scal=1.0;
                Vector<T> output;
                output.copy(input);
@@ -1239,7 +1239,7 @@ namespace FISTA {
             T _thrs;
       };
 
-   template <typename T> 
+   template <typename T>
       class Lzero : public Regularizer<T> {
          public:
             Lzero(const ParamReg<T>& param) : Regularizer<T>(param) { };
@@ -1252,13 +1252,13 @@ namespace FISTA {
                y.hardThrshold(sqrt(2*lambda));
                if (this->_intercept) y[y.n()-1] = x[y.n()-1];
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                return (this->_intercept ? x.lzero() - 1 : x.lzero());
             };
             void inline fenchel(const Vector<T>& input, T& val, T& scal) const { };
       };
 
-   template <typename T> 
+   template <typename T>
       class LogDC : public Regularizer<T> {
          public:
             LogDC(const ParamReg<T>& param) : Regularizer<T>(param), _eps(param.lambda2d1) { };
@@ -1275,7 +1275,7 @@ namespace FISTA {
                for (int i = 0; i<x.n(); ++i) _weights[i] = T(1.0)/(abs<T>(x[i])+_eps);
             };
             bool inline is_concave() const { return true; };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                T tmp=0;
                for (int i = 0; i<x.n(); ++i) tmp+= log_alt<T>(abs<T>(x[i])+_eps);
                return tmp;
@@ -1287,7 +1287,7 @@ namespace FISTA {
       };
 
 
-   template <typename T> 
+   template <typename T>
       class None: public Regularizer<T>, public SplittingFunction<T, SpMatrix<T> > {
          public:
             None() { };
@@ -1302,9 +1302,9 @@ namespace FISTA {
             void inline fenchel(const Vector<T>& input, T& val, T& scal) const {  };
             virtual bool is_fenchel() const { return false; };
             virtual bool is_subgrad() const { return true; };
-            virtual void sub_grad(const Vector<T>& input, Vector<T>& output) const { 
+            virtual void sub_grad(const Vector<T>& input, Vector<T>& output) const {
                output.setZeros();
-            } 
+            }
             virtual void reset() { };
             virtual T eval_split(const SpMatrix<T>& input) const { return 0; };
             virtual int num_components() const { return 0; };
@@ -1314,7 +1314,7 @@ namespace FISTA {
 //            virtual bool is_none() const { return true; };
       };
 
-   template <typename T> 
+   template <typename T>
       class Ridge: public Regularizer<T> {
          public:
             Ridge(const ParamReg<T>& param) : Regularizer<T>(param) { this->_id = RIDGE; };
@@ -1326,7 +1326,7 @@ namespace FISTA {
                y.scal(T(1.0/(1.0+lambda)));
                if (this->_intercept) y[y.n()-1] = x[y.n()-1];
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                return (this->_intercept ? 0.5*x.nrm2sq() - 0.5*x[x.n()-1]*x[x.n()-1] : 0.5*x.nrm2sq());
             };
             void inline fenchel(const Vector<T>& input, T& val, T& scal) const {
@@ -1335,10 +1335,10 @@ namespace FISTA {
                if (this->_pos) tmp.thrsPos();
                val=this->eval(tmp);
                scal=T(1.0);
-               if (this->_intercept & (abs<T>(tmp[tmp.n()-1]) > EPSILON)) val=INFINITY; 
+               if (this->_intercept & (abs<T>(tmp[tmp.n()-1]) > EPSILON)) val=INFINITY;
             };
             virtual bool is_subgrad() const { return true; };
-            virtual void sub_grad(const Vector<T>& input, Vector<T>& output) const {  
+            virtual void sub_grad(const Vector<T>& input, Vector<T>& output) const {
                output.resize(input.n());
                if (!this->_pos) {
                   for (int i = 0; i<input.n(); ++i) {
@@ -1352,7 +1352,7 @@ namespace FISTA {
             }
       };
 
-   template <typename T> 
+   template <typename T>
       class normL2: public Regularizer<T> {
          public:
             normL2(const ParamReg<T>& param) : Regularizer<T>(param) { };
@@ -1370,9 +1370,9 @@ namespace FISTA {
                }
                if (this->_intercept) y[y.n()-1] = x[y.n()-1];
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                Vector<T> xref(x.rawX(),this->_intercept ? x.n()-1 : x.n());
-               return xref.nrm2(); 
+               return xref.nrm2();
             };
             /// TODO add subgradient
             void inline fenchel(const Vector<T>& input, T& val, T& scal) const {
@@ -1382,11 +1382,11 @@ namespace FISTA {
                T mm = output.nrm2();
                scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                val=0;
-               if (this->_intercept & (abs<T>(output[output.n()-1]) > EPSILON)) val=INFINITY; 
+               if (this->_intercept & (abs<T>(output[output.n()-1]) > EPSILON)) val=INFINITY;
             };
       };
 
-   template <typename T> 
+   template <typename T>
       class normLINF: public Regularizer<T> {
          public:
             normLINF(const ParamReg<T>& param) : Regularizer<T>(param) { };
@@ -1402,9 +1402,9 @@ namespace FISTA {
                   y[j]=y[j]-row[j];
                if (this->_intercept) y[y.n()-1] = x[y.n()-1];
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                Vector<T> xref(x.rawX(),this->_intercept ? x.n()-1 : x.n());
-               return xref.fmaxval(); 
+               return xref.fmaxval();
             };
             /// TODO add subgradient
             void inline fenchel(const Vector<T>& input, T& val, T& scal) const {
@@ -1414,7 +1414,7 @@ namespace FISTA {
                T mm = output.asum();
                scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                val=0;
-               if (this->_intercept & (abs<T>(output[output.n()-1]) > EPSILON)) val=INFINITY; 
+               if (this->_intercept & (abs<T>(output[output.n()-1]) > EPSILON)) val=INFINITY;
             };
       };
 
@@ -1448,13 +1448,13 @@ namespace FISTA {
                   }
                }
             };
-            T inline eval(const D& x) const { 
+            T inline eval(const D& x) const {
                return _regA->eval(x) + _lambda2d1*_regB->eval(x);
             };
             virtual bool is_fenchel() const { return false; };
             void inline fenchel(const D& input, T& val, T& scal) const { };
             virtual bool is_subgrad() const { return _regA->is_subgrad() && _regB->is_subgrad(); };
-            virtual void sub_grad(const D& input, D& output) const {  
+            virtual void sub_grad(const D& input, D& output) const {
                _regA->sub_grad(input,output);
                D tmp;
                _regB->sub_grad(input,tmp);
@@ -1471,7 +1471,7 @@ namespace FISTA {
          typedef ComposeProx< T, Vector<T>, Lasso<T>, Ridge<T>, true > type;
       };
 
-   template <typename T> 
+   template <typename T>
       class FusedLasso: public Regularizer<T> {
          public:
             FusedLasso(const ParamReg<T>& param) : Regularizer<T>(param) {
@@ -1486,7 +1486,7 @@ namespace FISTA {
                copyx.copy(x);
                copyx.fusedProjectHomotopy(y,_lambda2d1*lambda,lambda,_lambda3d1*lambda,true);
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                T sum = T();
                const int maxn = this->_intercept ? x.n()-1 : x.n();
                for (int i = 0; i<maxn-1; ++i)
@@ -1502,10 +1502,10 @@ namespace FISTA {
             T _lambda3d1;
       };
 
-   template <typename T> 
+   template <typename T>
       class GraphLasso : public Regularizer<T>, public SplittingFunction<T, SpMatrix<T> > {
          public:
-            GraphLasso(const ParamReg<T>& param) : Regularizer<T>(param) { 
+            GraphLasso(const ParamReg<T>& param) : Regularizer<T>(param) {
                const bool resetflow = param.resetflow;
                const bool linf = param.linf;
                const bool clever = param.clever;
@@ -1554,7 +1554,7 @@ namespace FISTA {
                _old_lambda=lambda;
             };
 
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                Graph<T>* gr = const_cast<Graph<T>* >(&_graph);
                gr->restore_capacities();
                return gr->norm(x.rawX(),_work.rawX(),_weights.rawX(),_linf);
@@ -1577,7 +1577,7 @@ namespace FISTA {
                   gr->restore_flow();
                scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                val=0;
-               if (this->_intercept & (abs<T>(input[input.n()-1]) > EPSILON)) val=INFINITY; 
+               if (this->_intercept & (abs<T>(input[input.n()-1]) > EPSILON)) val=INFINITY;
             };
 
             virtual void init(const Vector<T>& y) { };
@@ -1588,7 +1588,7 @@ namespace FISTA {
                if (_linf) {
                   for (int i = 0; i<splitted_w.n(); ++i) {
                      splitted_w.refCol(i,col);
-                     tmp.setData(col.rawX(),col.nzmax());                     
+                     tmp.setData(col.rawX(),col.nzmax());
                      Vector<T> res;
                      res.copy(tmp);
                      vAbs<T>(res.n(),res.rawX(),res.rawX());
@@ -1598,7 +1598,7 @@ namespace FISTA {
                } else {
                   for (int i = 0; i<splitted_w.n(); ++i) {
                      splitted_w.refCol(i,col);
-                     tmp.setData(col.rawX(),col.nzmax());                     
+                     tmp.setData(col.rawX(),col.nzmax());
                      const T nrm = tmp.nrm2();
                      if (nrm > lambda*_weights[i]) {
                         tmp.scal(T(1.0)-lambda*_weights[i]/nrm);
@@ -1652,7 +1652,7 @@ namespace FISTA {
          typedef ComposeProx<T, Vector<T>, GraphLasso<T>, Ridge<T>, true> type;
       };
 
-   template <typename T> 
+   template <typename T>
       class TreeLasso : public Regularizer<T> {
          public:
             TreeLasso(const ParamReg<T>& param) : Regularizer<T>(param) {
@@ -1677,7 +1677,7 @@ namespace FISTA {
                }
                _tree.proj(yp,_linf,lambda);
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                return const_cast<Tree_Seq<T>* >(&_tree)->val_norm(x.rawX(),0,_linf);
             };
             void inline fenchel(const Vector<T>& y, T& val, T& scal) const {
@@ -1694,8 +1694,8 @@ namespace FISTA {
                   T mm = const_cast<Tree_Seq<T>* >(&_tree)->dual_norm_inf(yp2);
                   scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                   val=0;
-                  if (this->_intercept & (abs<T>(y[y.n()-1]) > EPSILON)) val=INFINITY; 
-               } 
+                  if (this->_intercept & (abs<T>(y[y.n()-1]) > EPSILON)) val=INFINITY;
+               }
             };
             virtual bool is_fenchel() const {
                return _linf;
@@ -1711,7 +1711,7 @@ namespace FISTA {
             bool _linf;
       };
 
-   template <typename T> 
+   template <typename T>
       class TreeLzero : public Regularizer<T> {
          public:
             TreeLzero(const ParamReg<T>& param) : Regularizer<T>(param) {
@@ -1734,7 +1734,7 @@ namespace FISTA {
                }
                _tree.proj_zero(yp,lambda);
             };
-            T inline eval(const Vector<T>& x) const { 
+            T inline eval(const Vector<T>& x) const {
                return const_cast<Tree_Seq<T>* >(&_tree)->val_zero(x.rawX(),0);
             };
             virtual bool is_fenchel() const { return false; };
@@ -1747,7 +1747,7 @@ namespace FISTA {
    template <typename T, typename ProxMat>
       class ProxMatToVec : public Regularizer<T> {
          public:
-            ProxMatToVec(const ParamReg<T>& param) : Regularizer<T>(param) { 
+            ProxMatToVec(const ParamReg<T>& param) : Regularizer<T>(param) {
                _size_group=param.size_group;
                ParamReg<T> param2=param;
                param2.intercept=false;
@@ -1792,12 +1792,12 @@ namespace FISTA {
                   for (int i = 0; i<param.ngroups; ++i) num_groups=MAX(num_groups,param.groups[i]);
                   _groups.resize(num_groups);
                   for (int i = 0; i<num_groups; ++i) _groups[i]=new list_int();
-                  for (int i = 0; i<param.ngroups; ++i) _groups[param.groups[i]-1]->push_back(i); 
-               } 
+                  for (int i = 0; i<param.ngroups; ++i) _groups[param.groups[i]-1]->push_back(i);
+               }
                _prox = new Reg(param2);
             }
-            virtual ~GroupProx() { 
-               delete(_prox); 
+            virtual ~GroupProx() {
+               delete(_prox);
                for (int i = 0; i<static_cast<int>(_groups.size()); ++i) delete(_groups[i]);
             };
 
@@ -1854,7 +1854,7 @@ namespace FISTA {
                return sum;
             }
             virtual bool is_fenchel() const { return _prox->is_fenchel(); };
-            void inline fenchel(const Vector<T>& x, T& val, T& scal) const { 
+            void inline fenchel(const Vector<T>& x, T& val, T& scal) const {
                const int maxn= this->_intercept ? x.n()-1 : x.n();
                T val2;
                T scal2;
@@ -1926,13 +1926,13 @@ namespace FISTA {
                for (int i = 0; i<m; ++i) {
                   if (norm[i] > lambda) {
                      T scal = (norm[i]-lambda)/norm[i];
-                     for (int j = 0; j<n; ++j) 
+                     for (int j = 0; j<n; ++j)
                         y[j*m+i] = x[j*m+i]*scal;
                   }
                }
                if (this->_pos) y.thrsPos();
                if (this->_intercept)
-                  for (int j = 0; j<n; ++j) 
+                  for (int j = 0; j<n; ++j)
                      y[j*m+m-1]=x[j*m+m-1];
             }
             T inline eval(const Matrix<T>& x) const {
@@ -1941,7 +1941,7 @@ namespace FISTA {
                return this->_intercept ? norm.asum() - norm[norm.n() -1] : norm.asum();
             }
             virtual bool is_subgrad() const { return true; };
-            virtual void sub_grad(const Matrix<T>& input, Matrix<T>& output) const { 
+            virtual void sub_grad(const Matrix<T>& input, Matrix<T>& output) const {
                Vector<T> norm;
                input.norm_2_rows(norm);
                for (int i = 0; i<norm.n(); ++i) {
@@ -1963,9 +1963,9 @@ namespace FISTA {
                   input.norm_2_rows(norm);
                }
                T mm = norm.fmaxval();
-               scal= mm > 1.0 ? T(1.0)/mm : 1.0; 
+               scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                val=0;
-               if (this->_intercept & (abs<T>(norm[norm.n()-1]) > EPSILON)) val=INFINITY; 
+               if (this->_intercept & (abs<T>(norm[norm.n()-1]) > EPSILON)) val=INFINITY;
             };
       };
 
@@ -1991,7 +1991,7 @@ namespace FISTA {
                      y(i,j) = row[j]-row2[j];
                }
             }
-            T inline eval(const Matrix<T>& x) const { 
+            T inline eval(const Matrix<T>& x) const {
                Vector<T> norm;
                x.norm_inf_rows(norm);
                return this->_intercept ? norm.asum() - norm[norm.n() -1] : norm.asum();
@@ -2008,12 +2008,12 @@ namespace FISTA {
                }
                if (this->_intercept) norm[norm.n()-1]=0;
                T mm = norm.fmaxval();
-               scal= mm > 1.0 ? T(1.0)/mm : 1.0; 
+               scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                val=0;
-               if (this->_intercept & (abs<T>(norm[norm.n()-1]) > EPSILON)) val=INFINITY; 
+               if (this->_intercept & (abs<T>(norm[norm.n()-1]) > EPSILON)) val=INFINITY;
             };
             virtual bool is_subgrad() const { return true; };
-            virtual void sub_grad(const Matrix<T>& input, Matrix<T>& output) const { 
+            virtual void sub_grad(const Matrix<T>& input, Matrix<T>& output) const {
                output.resize(input.m(),input.n());
                output.setZeros();
                const T maxm= this->_intercept ? input.m()-1 : input.m();
@@ -2024,12 +2024,12 @@ namespace FISTA {
                   if (max > 1e-15) {
                      int num_max=0;
                      for (int j = 0; j<row.n(); ++j) {
-                        if (abs<T>(max-abs<T>(row[j])) < 1e-15) 
+                        if (abs<T>(max-abs<T>(row[j])) < 1e-15)
                            num_max++;
                      }
                      T add = T(1.0)/num_max;
                      for (int j = 0; j<row.n(); ++j) {
-                        if (abs<T>(max-abs<T>(row[j])) < 1e-15) 
+                        if (abs<T>(max-abs<T>(row[j])) < 1e-15)
                            row[j] = row[j] > 0 ? add : -add;
                      }
                      output.setRow(i,row);
@@ -2041,7 +2041,7 @@ namespace FISTA {
    template <typename T>
       class TraceNorm : public Regularizer<T,Matrix<T> > {
          public:
-            TraceNorm(const ParamReg<T>& param) : Regularizer<T,Matrix<T> >(param) { 
+            TraceNorm(const ParamReg<T>& param) : Regularizer<T,Matrix<T> >(param) {
                if (param.intercept) {
                   cerr << "Trace norm implementation is not compatible with intercept, intercept deactivated" << endl;
                }
@@ -2053,7 +2053,7 @@ namespace FISTA {
             virtual ~TraceNorm() { };
 
             void inline prox(const Matrix<T>& x, Matrix<T>& y, const T lambda) {
-               //Matrix<T> tmp;               
+               //Matrix<T> tmp;
                //tmp.copy(x);
                Matrix<T> U;
                Matrix<T> V;
@@ -2113,7 +2113,7 @@ namespace FISTA {
    template <typename T>
       class Rank : public Regularizer<T,Matrix<T> > {
          public:
-            Rank(const ParamReg<T>& param) : Regularizer<T,Matrix<T> >(param) { 
+            Rank(const ParamReg<T>& param) : Regularizer<T,Matrix<T> >(param) {
                if (param.intercept) {
                   cerr << "Rank implementation is not compatible with intercept, intercept deactivated" << endl;
                }
@@ -2125,7 +2125,7 @@ namespace FISTA {
             virtual ~Rank() { };
 
             void inline prox(const Matrix<T>& x, Matrix<T>& y, const T lambda) {
-               Matrix<T> tmp;               
+               Matrix<T> tmp;
                tmp.copy(x);
                y.resize(x.m(),x.n());
                y.setZeros();
@@ -2175,9 +2175,9 @@ namespace FISTA {
           int count_col=0;
           int count=0;
           pB[0]=0;
-          for (ListIterator<Path<long long>*> it_path=paths.begin(); 
+          for (ListIterator<Path<long long>*> it_path=paths.begin();
                 it_path != paths.end(); ++it_path) {
-             for (const_iterator_int it = it_path->nodes.begin(); 
+             for (const_iterator_int it = it_path->nodes.begin();
                    it != it_path->nodes.end(); ++it) {
                 r[count]= *it;
                 v[count++]= it_path->flow;
@@ -2186,8 +2186,8 @@ namespace FISTA {
           }
           for (int i = 0; i<paths_mat.n(); ++i) sort(r,v,pB[i],pE[i]-1);
        };
- 
-    template <typename T> 
+
+    template <typename T>
        class GraphPathL0 : public Regularizer<T> {
           public:
              GraphPathL0(const ParamReg<T>& param) : Regularizer<T>(param) {
@@ -2195,33 +2195,33 @@ namespace FISTA {
                 _graph.init_graph(graph);
              }
              virtual ~GraphPathL0() { };
- 
+
              void inline prox(const Vector<T>& x, Vector<T>& y, const T lambda) {
                 // DEBUG
                 y.copy(x);
                 if (this->_pos) y.thrsPos();
                 _graph.proximal_l0(y.rawX(),lambda);
              };
-             T inline eval(const Vector<T>& x) const { 
+             T inline eval(const Vector<T>& x) const {
                 return const_cast<GraphPath<T>* >(&_graph)->eval_l0(x.rawX());
              };
-             T inline eval_paths(const Vector<T>& x, SpMatrix<T>& paths_mat) const { 
+             T inline eval_paths(const Vector<T>& x, SpMatrix<T>& paths_mat) const {
                 List<Path<long long>*> paths;
                 T val=const_cast<GraphPath<T>* >(&_graph)->eval_l0(x.rawX(),&paths);
                 convert_paths_to_mat<T>(paths,paths_mat,_graph.n());
-                for (ListIterator<Path<>*> it_path=paths.begin(); 
+                for (ListIterator<Path<>*> it_path=paths.begin();
                       it_path != paths.end(); ++it_path) delete(*it_path);
                 return val;
              };
- 
+
              virtual bool is_fenchel() const { return false; };
              void inline fenchel(const Vector<T>& input, T& val, T& scal) const {  };
- 
+
           private:
              GraphPath<T> _graph;
        };
- 
-    template <typename T> 
+
+    template <typename T>
        class GraphPathConv : public Regularizer<T> {
           public:
              GraphPathConv(const ParamReg<T>& param) : Regularizer<T>(param) {
@@ -2229,27 +2229,27 @@ namespace FISTA {
                 _graph.init_graph(graph);
              }
              virtual ~GraphPathConv() {  };
- 
+
              void inline prox(const Vector<T>& x, Vector<T>& y, const T lambda) {
                 y.copy(x);
                 if (this->_pos) y.thrsPos();
                 _graph.proximal_conv(y.rawX(),lambda);
              };
-             T inline eval(const Vector<T>& x) const { 
+             T inline eval(const Vector<T>& x) const {
                 return const_cast<GraphPath<T>* >(&_graph)->eval_conv(x.rawX());
              };
-             T inline eval_dual_norm(const Vector<T>& x) const { 
+             T inline eval_dual_norm(const Vector<T>& x) const {
                 return const_cast<GraphPath<T>* >(&_graph)->eval_dual_norm(x.rawX(),NULL);
              };
-             T inline eval_paths(const Vector<T>& x, SpMatrix<T>& paths_mat) const { 
+             T inline eval_paths(const Vector<T>& x, SpMatrix<T>& paths_mat) const {
                 List<Path<long long>*> paths;
                 T val=const_cast<GraphPath<T>* >(&_graph)->eval_conv(x.rawX(),&paths);
                 convert_paths_to_mat<T>(paths,paths_mat,_graph.n());
-                for (ListIterator<Path<long long>*> it_path=paths.begin(); 
+                for (ListIterator<Path<long long>*> it_path=paths.begin();
                       it_path != paths.end(); ++it_path) delete(*it_path);
                 return val;
              };
-             T inline eval_dual_norm_paths(const Vector<T>& x, SpMatrix<T>& paths_mat) const { 
+             T inline eval_dual_norm_paths(const Vector<T>& x, SpMatrix<T>& paths_mat) const {
                 Path<long long> path;
                 T val=const_cast<GraphPath<T>* >(&_graph)->eval_dual_norm(x.rawX(),&path.nodes);
                 List<Path<long long>*> paths;
@@ -2273,7 +2273,7 @@ namespace FISTA {
                 }
                 scal= mm > 1.0 ? T(1.0)/mm : 1.0;
                 val=0;
-                if (this->_intercept & (abs<T>(input[input.n()-1]) > EPSILON)) val=INFINITY; 
+                if (this->_intercept & (abs<T>(input[input.n()-1]) > EPSILON)) val=INFINITY;
              };
           private:
              GraphPath<T> _graph;
@@ -2283,29 +2283,29 @@ namespace FISTA {
    template <typename T,typename Reg>
       class RegMat : public Regularizer<T,Matrix<T> > {
          public:
-            RegMat(const ParamReg<T>& param) : Regularizer<T,Matrix<T> >(param) { 
+            RegMat(const ParamReg<T>& param) : Regularizer<T,Matrix<T> >(param) {
                _transpose=param.transpose;
                const int N = param.num_cols;
                _regs=new Reg*[N];
                _N=N;
-               for (int i = 0; i<N; ++i) 
+               for (int i = 0; i<N; ++i)
                   _regs[i]=new Reg(param);
             };
-            virtual ~RegMat() { 
+            virtual ~RegMat() {
                for (int i = 0; i<_N; ++i) {
                   delete(_regs[i]);
                   _regs[i]=NULL;
                }
                delete[](_regs);
             };
-            void inline reset() { 
+            void inline reset() {
                for (int i = 0; i<_N; ++i) _regs[i]->reset();
             };
             void inline prox(const Matrix<T>& x, Matrix<T>& y, const T lambda) {
                y.copy(x);
                int i;
                if (_transpose) {
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
                   for (i = 0; i<_N; ++i) {
                      Vector<T> colx, coly;
                      x.copyRow(i,colx);
@@ -2313,7 +2313,7 @@ namespace FISTA {
                      y.setRow(i,coly);
                   }
                } else {
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
                   for (i = 0; i<_N; ++i) {
                      Vector<T> colx, coly;
                      x.refCol(i,colx);
@@ -2322,9 +2322,9 @@ namespace FISTA {
                   }
                }
             };
-            virtual bool is_subgrad() const { 
+            virtual bool is_subgrad() const {
                bool ok=true;
-               for (int i = 0; i<_N; ++i) 
+               for (int i = 0; i<_N; ++i)
                   ok=ok && _regs[i]->is_subgrad();
                return ok;
             };
@@ -2345,10 +2345,10 @@ namespace FISTA {
                   }
                }
             };
-            T inline eval(const Matrix<T>& x) const { 
+            T inline eval(const Matrix<T>& x) const {
                T sum = 0;
                int i;
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
                for (i = 0; i<_N; ++i) {
                   Vector<T> col;
                   if (_transpose) {
@@ -2380,7 +2380,7 @@ namespace FISTA {
             };
             virtual bool is_fenchel() const {
                bool ok=true;
-               for (int i = 0; i<_N; ++i) 
+               for (int i = 0; i<_N; ++i)
                   ok = ok && _regs[i]->is_fenchel();
                return ok;
             };
@@ -2418,7 +2418,7 @@ namespace FISTA {
                y.toVect(yv);
                _graphlasso->prox(xv,yv,lambda);
             }
-            T inline eval(const Matrix<T>& X) const { 
+            T inline eval(const Matrix<T>& X) const {
                Vector<T> xv;
                X.toVect(xv);
                return _graphlasso->eval(xv);
@@ -2440,7 +2440,7 @@ namespace FISTA {
    template <typename T>
       class MixedL1LINFCR : public SpecGraphMat<T> {
          public:
-            MixedL1LINFCR(const int m, const ParamReg<T>& param) : SpecGraphMat<T>(param) { 
+            MixedL1LINFCR(const int m, const ParamReg<T>& param) : SpecGraphMat<T>(param) {
                const int n = param.num_cols;
                const T l2dl1 = param.lambda2d1;
                GraphStruct<T> graph_st;
@@ -2491,7 +2491,7 @@ namespace FISTA {
    template <typename T>
       class TreeMult : public SpecGraphMat<T> {
          public:
-            TreeMult(const ParamReg<T>& param) : SpecGraphMat<T>(param) { 
+            TreeMult(const ParamReg<T>& param) : SpecGraphMat<T>(param) {
                const TreeStruct<T>& tree_st=*(param.tree_st);
                const int N = param.num_cols;
                const T l1dl2 = param.lambda2d1;
@@ -2575,7 +2575,7 @@ namespace FISTA {
    template <typename T>
       class GraphMult : public SpecGraphMat<T> {
          public:
-            GraphMult(const ParamReg<T>& param) : SpecGraphMat<T>(param) { 
+            GraphMult(const ParamReg<T>& param) : SpecGraphMat<T>(param) {
                const GraphStruct<T>& graph_st=*(param.graph_st);
                const int N = param.num_cols;
                const T l1dl2 = param.lambda2d1;
@@ -2600,7 +2600,7 @@ namespace FISTA {
                   for (int j = 0; j<Ng; ++j) {
                      gv_jc[i*Ng+j]=count;
                      for (mwSize k = graph_st.gv_jc[j]; k<graph_st.gv_jc[j+1]; ++k) {
-                        gv_ir[count++] =Nv*i+graph_st.gv_ir[k]; 
+                        gv_ir[count++] =Nv*i+graph_st.gv_ir[k];
                      }
                   }
                }
@@ -2651,7 +2651,7 @@ namespace FISTA {
       };
 
    template <typename T, typename D, typename E>
-      T duality_gap(Loss<T,D,E>& loss, Regularizer<T,D>& regularizer, const D& x, 
+      T duality_gap(Loss<T,D,E>& loss, Regularizer<T,D>& regularizer, const D& x,
             const T lambda, T& best_dual, const bool verbose = false) {
          if (!regularizer.is_fenchel() || !loss.is_fenchel()) {
             cerr << "Error: no duality gap available" << endl;
@@ -2660,7 +2660,7 @@ namespace FISTA {
          T primal= loss.eval(x)+lambda*regularizer.eval(x);
          bool intercept=regularizer.is_intercept();
          D grad1, grad2;
-         loss.var_fenchel(x,grad1,grad2,intercept); 
+         loss.var_fenchel(x,grad1,grad2,intercept);
          T dual;
          grad2.scal(-T(1.0)/lambda);
          T val=0;
@@ -2680,7 +2680,7 @@ namespace FISTA {
       }
 
    template <typename T, typename D, typename E>
-      T duality_gap(Loss<T,D,E>& loss, Regularizer<T,D>& regularizer, const D& x, 
+      T duality_gap(Loss<T,D,E>& loss, Regularizer<T,D>& regularizer, const D& x,
             const T lambda, const bool verbose = false) {
          T best_dual=-INFINITY;
          return duality_gap(loss,regularizer,x,lambda,best_dual,verbose);
@@ -2688,7 +2688,7 @@ namespace FISTA {
 
    template <typename T>
       void dualityGraph(const Matrix<T>& X, const Matrix<T>& D, const Matrix<T>& alpha0,
-            Vector<T>& res, const ParamFISTA<T>& param, 
+            Vector<T>& res, const ParamFISTA<T>& param,
             const GraphStruct<T>* graph_st) {
          Regularizer<T>* regularizer=new GraphLasso<T>(*graph_st,
                param.intercept,param.resetflow,param.pos,param.clever);
@@ -2725,7 +2725,7 @@ namespace FISTA {
 
 
    template <typename T, typename D, typename E>
-      void subGradientDescent_Generic(Loss<T,D,E>& loss, Regularizer<T,D>& regularizer, const D& x0, D& x, 
+      void subGradientDescent_Generic(Loss<T,D,E>& loss, Regularizer<T,D>& regularizer, const D& x0, D& x,
             Vector<T>& optim_info,
             const ParamFISTA<T>& param) {
          D grad;
@@ -2743,13 +2743,13 @@ namespace FISTA {
             /// print loss
             if (param.verbose && ((it % it0) == 0)) {
                time.stop();
-               T los=loss.eval(x) + lambda*regularizer.eval(x);    
+               T los=loss.eval(x) + lambda*regularizer.eval(x);
                optim_info[0]=los;
                T sec=time.getElapsed();
                cout << "Iter: " << it << ", loss: " << los << ", time: " << sec << " ";
-               if (param.log) 
+               if (param.log)
                   writeLog(it,sec,los,best_dual,param.logName);
-               if (param.verbose) 
+               if (param.verbose)
                   cout << endl;
                flush(cout);
                time.start();
@@ -2771,7 +2771,7 @@ namespace FISTA {
             }
          }
          if ((it % it0) != 0 || !param.verbose) {
-            T los=loss.eval(x) + lambda*regularizer.eval(x);    
+            T los=loss.eval(x) + lambda*regularizer.eval(x);
             optim_info[0]=los;
             if (duality) {
                rel_duality_gap=duality_gap(loss,regularizer,x,lambda,best_dual,param.verbose);
@@ -2817,7 +2817,7 @@ namespace FISTA {
                T sec=time.getElapsed();
                cout << "Iter: " << it << ", loss: " << los << ", time: " << sec << ", L: " << L;
                flush(cout);
-               if (param.log) 
+               if (param.log)
                   writeLog(it,sec,los,best_dual,param.logName);
                time.start();
             }
@@ -2843,7 +2843,7 @@ namespace FISTA {
                   break;
                }
                L *= param.gamma;
-               if (param.verbose && ((it % it0) == 0)) 
+               if (param.verbose && ((it % it0) == 0))
                   cout << " " << L;
                ++iter;
             }
@@ -2860,12 +2860,12 @@ namespace FISTA {
                      regularizer.prox(prox,tmp,lambda/L);
                      break;
                   }
-                  if (param.verbose && ((it % it0) == 0)) 
+                  if (param.verbose && ((it % it0) == 0))
                   cout << " " << L;
                   ++iter;
                }
             }
-            if (param.verbose && ((it % it0) == 0)) 
+            if (param.verbose && ((it % it0) == 0))
                cout << endl;
             if (param.linesearch_mode==2) {
                sbb.copy(grad);
@@ -2887,7 +2887,7 @@ namespace FISTA {
                if (sqrt(old.nrm2sq()/MAX(EPSILON,x.nrm2sq())) < param.tol) break;
             }
          }
-         T los=loss.eval(x) + lambda*regularizer.eval(x);    
+         T los=loss.eval(x) + lambda*regularizer.eval(x);
          optim_info[0]=los;
          T sec=time.getElapsed();
          if (param.verbose) {
@@ -2931,12 +2931,12 @@ namespace FISTA {
             /// print loss
             if (param.verbose && ((it % it0) == 0)) {
                time.stop();
-               T los=loss.eval(x) + lambda*regularizer.eval(x);    
+               T los=loss.eval(x) + lambda*regularizer.eval(x);
                optim_info[0]=los;
                T sec=time.getElapsed();
                cout << "Iter: " << it << ", loss: " << los << ", time: " << sec << ", L: " << L;
                flush(cout);
-               if (param.log) 
+               if (param.log)
                   writeLog(it,sec,los,best_dual,param.logName);
                time.start();
             }
@@ -2963,11 +2963,11 @@ namespace FISTA {
                regularizer.prox(prox,tmp,lambda/L);
                if ((param.linesearch_mode==2 && it > 1) || param.fixed_step || loss.test_backtracking(y,grad,tmp,L)) break;
                L *= param.gamma;
-               if (param.verbose && ((it % it0) == 0)) 
+               if (param.verbose && ((it % it0) == 0))
                   cout << " " << L;
                ++iter;
             }
-            if (param.verbose && ((it % it0) == 0)) 
+            if (param.verbose && ((it % it0) == 0))
                cout << endl;
 
             prox.copy(x);
@@ -3007,7 +3007,7 @@ namespace FISTA {
       };
 
    template <typename T>
-      T LagrangianADMM(const SplittingFunction<T, Matrix<T> >& loss, const SplittingFunction<T, SpMatrix<T> >& reg, 
+      T LagrangianADMM(const SplittingFunction<T, Matrix<T> >& loss, const SplittingFunction<T, SpMatrix<T> >& reg,
             const T lambda, const T gamma, const Vector<T>& w,  const Matrix<T>& splitted_loss, const SpMatrix<T>& splitted_reg,
             const Matrix<T>& multi_loss, const SpMatrix<T>& multi_reg, T& los, const T* weights = NULL) {
          const int n_reg=reg.num_components();
@@ -3044,7 +3044,7 @@ namespace FISTA {
          splitted_w_loss.sum_cols(mean);
          w.copy(mean);
          multipliers_w_loss.sum_cols(mean);
-         w.add(mean,-T(1.0)/gamma);   
+         w.add(mean,-T(1.0)/gamma);
          Vector<T> number_occurences(w.n());
          number_occurences.set(splitted_w_loss.n());
          const int n_reg=splitted_w_reg.n();
@@ -3054,10 +3054,10 @@ namespace FISTA {
             for (int i = 0; i<n_reg; ++i) {
                splitted_w_reg.refCol(i,col);
                mean.add(col);
-               for (int j = 0; j<col.L(); ++j) 
+               for (int j = 0; j<col.L(); ++j)
                   number_occurences[col.r(j)]++;
             }
-            w.add(mean);   
+            w.add(mean);
             mean.setZeros();
             for (int i = 0; i<n_reg; ++i) {
                multipliers_w_reg.refCol(i,col);
@@ -3080,7 +3080,7 @@ namespace FISTA {
          splitted_w_loss.sum_cols(mean);
          w.copy(mean);
          multipliers_w_loss.sum_cols(mean);
-         w.add(mean,-T(1.0)/gamma);   
+         w.add(mean,-T(1.0)/gamma);
          Vector<T> number_occurences(w.n());
          number_occurences.set(splitted_w_loss.n());
          const int n_reg=splitted_w_reg.n();
@@ -3094,11 +3094,11 @@ namespace FISTA {
                   number_occurences[col.r(j)]+=inner_weights[j]*inner_weights[j];
                }
             }
-            w.add(mean);   
+            w.add(mean);
             mean.setZeros();
             for (int i = 0; i<n_reg; ++i) {
                multipliers_w_reg.refCol(i,col);
-               for (int j = 0; j<col.L(); ++j) 
+               for (int j = 0; j<col.L(); ++j)
                   mean[col.r(j)]+=inner_weights[j]*col.v(j);
             }
             w.add(mean,-T(1.0)/gamma);
@@ -3107,10 +3107,10 @@ namespace FISTA {
       };
 
    template <typename T>
-      void ADMM(const SplittingFunction<T, Matrix<T> >& loss, const SplittingFunction<T, SpMatrix<T> >& reg, 
+      void ADMM(const SplittingFunction<T, Matrix<T> >& loss, const SplittingFunction<T, SpMatrix<T> >& reg,
             const Vector<T>& w0, Vector<T>& w, Vector<T>& optim_info,
             const ParamFISTA<T>& param) {
-         const T gamma = param.c; 
+         const T gamma = param.c;
          const int n_reg=reg.num_components();
          const int it0 = MAX(1,param.it0);
          const T lambda=param.lambda;
@@ -3151,7 +3151,7 @@ namespace FISTA {
                if (param.verbose) {
                   cout << "Iter: " << it << ", loss: " << los << ", time: " << sec << endl;
                   flush(cout);
-                  if (param.log) 
+                  if (param.log)
                      writeLog(it,sec,los,T(0),param.logName);
                }
                time.start();
@@ -3236,7 +3236,7 @@ namespace FISTA {
             for (int i = 0; i<n_reg; ++i) {
                splitted_w_reg.refCol(i,col);
                mean.add(col);
-               for (int j = 0; j<col.L(); ++j) 
+               for (int j = 0; j<col.L(); ++j)
                   number_occurences[col.r(j)]+=gamma;
             }
             mean.scal(gamma);
@@ -3251,10 +3251,10 @@ namespace FISTA {
 
 
    template <typename T>
-      void LinADMM(const SplittingFunction<T, Matrix<T> >& loss, const SplittingFunction<T, SpMatrix<T> >& reg, 
+      void LinADMM(const SplittingFunction<T, Matrix<T> >& loss, const SplittingFunction<T, SpMatrix<T> >& reg,
             const Vector<T>& w0, Vector<T>& w, Vector<T>& optim_info,
             const ParamFISTA<T>& param) {
-         const T gamma = param.c; 
+         const T gamma = param.c;
          const int n_reg=reg.num_components();
          const int it0 = MAX(1,param.it0);
          const T lambda=param.lambda;
@@ -3294,7 +3294,7 @@ namespace FISTA {
                if (param.verbose) {
                   cout << "Iter: " << it << ", loss: " << los << ", time: " << sec << endl;
                   flush(cout);
-                  if (param.log) 
+                  if (param.log)
                      writeLog(it,sec,los,T(0),param.logName);
                }
                time.start();
@@ -3465,7 +3465,7 @@ namespace FISTA {
                   } else {
                      cout << "ADMM algorithm" << endl;
                   }
-               } 
+               }
             } else {
                if (param.ista) {
                   cout << "ISTA algorithm" << endl;
@@ -3476,7 +3476,7 @@ namespace FISTA {
                }
                if ((param.regul == GRAPH || param.regul == TREEMULT ||
                         param.regul == GRAPHMULT || param.regul==L1LINFCR) &&
-                     param.clever) 
+                     param.clever)
                   cout << "Projections with arc capacities" << endl;
                if (param.intercept) cout << "with intercept" << endl;
                if (param.pos) cout << "Non-negativity constraints" << endl;
@@ -3497,9 +3497,9 @@ namespace FISTA {
             SplittingFunction<T, Matrix<T> >** losses, const ParamFISTA<T>& param) {
          const int M = X.n();
          optim_info.resize(4,M);
-         
+
          int i1;
-#pragma omp parallel for private(i1) 
+#pragma omp parallel for private(i1)
          for (i1 = 0; i1< M; ++i1) {
 #ifdef _OPENMP
             int numT=omp_get_thread_num();
@@ -3522,7 +3522,7 @@ namespace FISTA {
                } else {
                   ADMM(*(losses[numT]),*(regularizers[numT]),alpha0i,alphai,optim_infoi,param);
                }
-            } 
+            }
          }
       }
 
@@ -3545,7 +3545,7 @@ namespace FISTA {
          optim_info.resize(4,M);
 
          int i1;
-#pragma omp parallel for private(i1) 
+#pragma omp parallel for private(i1)
          for (i1 = 0; i1< M; ++i1) {
 #ifdef _OPENMP
             int numT=omp_get_thread_num();
@@ -3586,7 +3586,7 @@ namespace FISTA {
          optim_info.resize(4,M);
 
          int i2;
-#pragma omp parallel for private(i2) 
+#pragma omp parallel for private(i2)
          for (i2 = 0; i2< M; ++i2) {
 #ifdef _OPENMP
             int numT=omp_get_thread_num();
@@ -3618,7 +3618,7 @@ namespace FISTA {
    template <typename T>
       void solver(const Matrix<T>& X, const AbstractMatrixB<T>& D, const Matrix<T>& alpha0,
             Matrix<T>& alpha, const ParamFISTA<T>& param1, Matrix<T>& optim_info,
-            const GraphStruct<T>* graph_st = NULL, 
+            const GraphStruct<T>* graph_st = NULL,
             const TreeStruct<T>* tree_st = NULL,
             const GraphPathStruct<T>* graph_path_st=NULL) {
          print_info_solver(param1);
@@ -3665,9 +3665,9 @@ namespace FISTA {
                   regularizers[i]=setRegularizerVectors(param,graph_st,tree_st,graph_path_st);
                   switch (param.loss) {
                      case SQUARE: if (param.compute_gram) {
-                                     losses[i]=new SqLoss<T>(D,G); 
+                                     losses[i]=new SqLoss<T>(D,G);
                                   } else {
-                                     losses[i]=new SqLoss<T>(D); 
+                                     losses[i]=new SqLoss<T>(D);
                                   }
                                   break;
                      case POISSON: losses[i]=new PoissonLoss<T>(D,param.delta);  break;
@@ -3717,16 +3717,16 @@ namespace FISTA {
                Regularizer<T, Matrix<T> >* regularizer;
                switch (param.loss) {
                   case SQUARE: if (param.compute_gram) {
-                                  loss=new SqLossMat<T>(D,G); 
+                                  loss=new SqLossMat<T>(D,G);
                                } else {
-                                  loss=new SqLossMat<T>(D); 
+                                  loss=new SqLossMat<T>(D);
                                }
                                break;
                   case POISSON: loss=new LossMat<T, PoissonLoss<T> >(X.n(),D,param.delta);  break;
                   case SQUARE_MISSING: loss=new LossMat<T, SqLossMissing<T> >(X.n(),D);  break;
                   case LOG:  loss = new LossMat<T, LogLoss<T,false> >(X.n(),D); break;
                   case LOGWEIGHT:  loss = new LossMat<T, LogLoss<T,true> >(X.n(),D); break;
-                  case CUR:  loss = new LossCur<T>(D); break; 
+                  case CUR:  loss = new LossCur<T>(D); break;
                   default: cerr << "Not implemented"; exit(1);
                }
                regularizer=setRegularizerMatrices(param,alpha0.m(),alpha0.n(),graph_st,tree_st,graph_path_st);
@@ -3753,16 +3753,16 @@ namespace FISTA {
 
    template <typename T>
       void PROX(const Matrix<T>& alpha0,
-            Matrix<T>& alpha, const ParamFISTA<T>& param, 
+            Matrix<T>& alpha, const ParamFISTA<T>& param,
             Vector<T>& val_loss,
-            const GraphStruct<T>* graph_st = NULL, 
+            const GraphStruct<T>* graph_st = NULL,
             const TreeStruct<T>* tree_st = NULL,
             const GraphPathStruct<T>* graph_path_st = NULL) {
          if (param.verbose) {
             print_regul(param.regul);
             if ((param.regul == GRAPH || param.regul == TREEMULT ||
                      param.regul == GRAPHMULT || param.regul==L1LINFCR) &&
-                  param.clever) 
+                  param.clever)
                cout << "Projections with arc capacities" << endl;
             if (param.intercept) cout << "with intercept" << endl;
             flush(cout);
@@ -3777,13 +3777,13 @@ namespace FISTA {
 
          if (!regul_for_matrices(param.regul)) {
             Regularizer<T>** regularizers= new Regularizer<T>*[num_threads];
-            for (int i = 0; i<num_threads; ++i) 
+            for (int i = 0; i<num_threads; ++i)
                regularizers[i]=setRegularizerVectors(param,graph_st,tree_st,graph_path_st);
 
             int i;
             if (param.eval)
                val_loss.resize(M);
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
             for (i = 0; i< M; ++i) {
 #ifdef _OPENMP
                int numT=omp_get_thread_num();
@@ -3820,7 +3820,7 @@ namespace FISTA {
 
    template <typename T>
       void EvalGraphPath(const Matrix<T>& alpha0,
-            const ParamFISTA<T>& param, 
+            const ParamFISTA<T>& param,
             Vector<T>& val_loss,
             const GraphPathStruct<T>* graph_path_st,
             SpMatrix<T>* paths = NULL) {
@@ -3836,12 +3836,12 @@ namespace FISTA {
 
          if (!regul_for_matrices(param.regul)) {
             Regularizer<T>** regularizers= new Regularizer<T>*[num_threads];
-            for (int i = 0; i<num_threads; ++i) 
+            for (int i = 0; i<num_threads; ++i)
                regularizers[i]=setRegularizerVectors<T>(param,NULL,NULL,graph_path_st);
 
             int i;
             val_loss.resize(M);
-#pragma omp parallel for private(i) 
+#pragma omp parallel for private(i)
             for (i = 0; i< M; ++i) {
 #ifdef _OPENMP
                int numT=omp_get_thread_num();

From 22486fcfd072247bc1f872a29f3dcff7c8e63e9d Mon Sep 17 00:00:00 2001
From: Samuel St-Jean <3030760+samuelstjean@users.noreply.github.com>
Date: Fri, 23 May 2025 07:44:30 -0400
Subject: [PATCH 3/3] Add changes to readme

---
 README.md | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/README.md b/README.md
index aba4d536..67bc7fe8 100644
--- a/README.md
+++ b/README.md
@@ -12,3 +12,14 @@ As such, I can probably help out with small stuff, but for technical and theoret
 You can find the original swig wrapper to re-generate these files on the branch swig_generator at https://github.com/samuelstjean/spams-python/tree/swig_generator
 
 This (unofficial) version includes pre-built wheels for python 3 on windows, mac (with openmp support) and linux and can be installed with ``pip install spams-bin``.
+
+# Improvements from the original version
+
++ Automatic build system with meson
++ Pytest for testing
++ Newer swig wrappers
++ Fixed internal naming clash for newer openmp namespace and the orthogonal matching pursuit (omp) function
++ Fixed C++20 deprecation issues with templates
++ Fixed undefined behavior with openmp and passing 0 as the number of cores (it now defaults to 1)
++ Fixed a positivity bug with spams.lassoWeighted
++ Probably other stuff