DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
linalg.h
Go to the documentation of this file.
1 /* Software SPAMS v2.1 - Copyright 2009-2011 Julien Mairal
2  *
3  * This file is part of SPAMS.
4  *
5  * SPAMS is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * SPAMS is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with SPAMS. If not, see <http://www.gnu.org/licenses/>.
17  */
18 
19 /* \file
20  * toolbox Linalg
21  *
22  * by Julien Mairal
23  * julien.mairal@inria.fr
24  *
25  * File linalg.h
26  * \brief Contains Matrix, Vector classes */
27 
28 #ifndef LINALG_H
29 #define LINALG_H
30 
31 #include "misc.h"
32 #ifdef USE_BLAS_LIB
33 #include "cblas_alt_template.h"
34 #else
35 #include "cblas_template.h" // this is obsolete
36 #endif
37 #include <fstream>
38 #ifdef WINDOWS
39 #include <string>
40 #else
41 #include <cstring>
42 #endif
43 #include <list>
44 #include <vector>
45 
46 #ifdef NEW_MATLAB
47  typedef ptrdiff_t INTT;
48 #else
49  typedef int INTT;
50 #endif
51 
52 #include <utils.h>
53 
54 #undef max
55 #undef min
56 
57 namespace spams
58 {
59 
61 template<typename T> class Matrix;
63 template<typename T> class SpMatrix;
65 template<typename T> class Vector;
67 template<typename T> class SpVector;
68 
69 
70 
71 typedef std::list< int > group;
72 typedef std::list< group > list_groups;
73 typedef std::vector< group > vector_groups;
74 
75 template <typename T>
76 static inline bool isZero(const T lambda) {
77  return static_cast<double>(abs<T>(lambda)) < 1e-99;
78 }
79 
80 template <typename T>
81 static inline bool isEqual(const T lambda1, const T lambda2) {
82  return static_cast<double>(abs<T>(lambda1-lambda2)) < 1e-99;
83 }
84 
85 
86 template <typename T>
87 static inline T softThrs(const T x, const T lambda) {
88  if (x > lambda) {
89  return x-lambda;
90  } else if (x < -lambda) {
91  return x+lambda;
92  } else {
93  return 0;
94  }
95 };
96 
97 template <typename T>
98 static inline T hardThrs(const T x, const T lambda) {
99  return (x > lambda || x < -lambda) ? x : 0;
100 };
101 
102 template <typename T>
103 static inline T alt_log(const T x);
104 template <> inline double alt_log<double>(const double x) { return log(x); };
105 template <> inline float alt_log<float>(const float x) { return logf(x); };
106 
107 template <typename T>
108 static inline T xlogx(const T x) {
109  if (x < -1e-20) {
110  return INFINITY;
111  } else if (x < 1e-20) {
112  return 0;
113  } else {
114  return x*alt_log<T>(x);
115  }
116 }
117 
118 template <typename T>
119 static inline T logexp(const T x) {
120  if (x < -30) {
121  return 0;
122  } else if (x < 30) {
123  return alt_log<T>( T(1.0) + exp_alt<T>( x ) );
124  } else {
125  return x;
126  }
127 }
128 
130 template <typename T> class Data {
131  public:
132  virtual void getData(Vector<T>& data, const int i) const = 0;
133  virtual void getGroup(Matrix<T>& data, const vector_groups& groups,
134  const int i) const = 0;
135  virtual inline T operator[](const int index) const = 0;
136  virtual int n() const = 0;
137  virtual int m() const = 0;
138  virtual int V() const = 0;
139  virtual void norm_2sq_cols(Vector<T>& norms) const { };
140  virtual ~Data() { };
141 };
142 
144 template <typename T> class AbstractMatrixB {
145  public:
146  virtual int n() const = 0;
147  virtual int m() const = 0;
148 
150  virtual void multTrans(const Vector<T>& x, Vector<T>& b,
151  const T alpha = 1.0, const T beta = 0.0) const = 0;
152 
154  virtual void mult(const SpVector<T>& x, Vector<T>& b,
155  const T alpha = 1.0, const T beta = 0.0) const = 0;
156 
157  virtual void mult(const Vector<T>& x, Vector<T>& b,
158  const T alpha = 1.0, const T beta = 0.0) const = 0;
159 
161  virtual void mult(const Matrix<T>& B, Matrix<T>& C,
162  const bool transA = false, const bool transB = false,
163  const T a = 1.0, const T b = 0.0) const = 0;
164 
165  virtual void mult(const SpMatrix<T>& B, Matrix<T>& C,
166  const bool transA = false, const bool transB = false,
167  const T a = 1.0, const T b = 0.0) const = 0;
168 
170  virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C,
171  const bool transA = false, const bool transB = false,
172  const T a = 1.0, const T b = 0.0) const = 0;
173 
175  virtual void XtX(Matrix<T>& XtX) const = 0;
176 
177  virtual void copyRow(const int i, Vector<T>& x) const = 0;
178 
179  virtual void copyTo(Matrix<T>& copy) const = 0;
180  virtual T dot(const Matrix<T>& x) const = 0;
181 
182  virtual void print(const string& name) const = 0;
183 
184  virtual ~AbstractMatrixB() { };
185 };
186 
188 template <typename T> class AbstractMatrix {
189  public:
190  virtual int n() const = 0;
191  virtual int m() const = 0;
193  virtual void copyCol(const int i, Vector<T>& Xi) const = 0;
195  virtual void add_rawCol(const int i, T* col, const T a) const = 0;
197  virtual void extract_rawCol(const int i,T* Xi) const = 0;
199  virtual void diag(Vector<T>& diag) const = 0;
201  virtual inline T operator()(const int index1, const int index2) const = 0;
202  virtual ~AbstractMatrix() { };
203 };
204 
206 template<typename T> class Matrix : public Data<T>, public AbstractMatrix<T>, public AbstractMatrixB<T> {
207  friend class SpMatrix<T>;
208  public:
209 
211  Matrix(T* X, int m, int n);
213  Matrix(int m, int n);
215  Matrix();
216 
218  virtual ~Matrix();
219 
222  inline int m() const { return _m; };
224  inline int n() const { return _n; };
226  inline T& operator()(const int i, const int j);
228  inline T operator()(const int i, const int j) const;
230  inline T& operator[](const int index) { return _X[index]; };
232  inline T operator[](const int index) const { return _X[index]; };
234  inline void copyCol(const int i, Vector<T>& x) const;
236  inline void copyRow(const int i, Vector<T>& x) const;
238  inline void extract_rawCol(const int i, T* x) const;
240  virtual void add_rawCol(const int i, T* DtXi, const T a) const;
242  inline void getData(Vector<T>& data, const int i) const;
244  virtual void getGroup(Matrix<T>& data, const vector_groups& groups,
245  const int i) const;
247  inline void refCol(int i, Vector<T>& x) const;
249  inline void refSubMat(int i, int n, Matrix<T>& mat) const;
251  inline void subMatrixSym(const Vector<int>& indices,
252  Matrix<T>& subMatrix) const;
254  inline T* rawX() const { return _X; };
256  inline const T* X() const { return _X; };
258  inline void copy(const Matrix<T>& mat);
260  inline void copyTo(Matrix<T>& mat) const { mat.copy(*this); };
262  inline void copyRef(const Matrix<T>& mat);
263 
266  inline void print(const string& name) const;
267 
268 
271  inline void clean();
273  inline void resize(int m, int n);
275  inline void setData(T* X, int m, int n);
277  inline void setm(const int m) { _m = m; }; //DANGEROUS
279  inline void setn(const int n) { _n = n; }; //DANGEROUS
281  inline void setZeros();
283  inline void set(const T a);
285  inline void clear();
287  inline void setAleat();
289  inline void eye();
291  inline void normalize();
293  inline void normalize2();
295  inline void center();
297  inline void center_rows();
299  inline void center(Vector<T>& centers);
301  inline void scal(const T a);
304  inline void fillSymmetric();
305  inline void fillSymmetric2();
307  inline void fakeSize(const int m, const int n) { _n = n; _m=m;};
309  inline void whiten(const int V);
311  inline void whiten(Vector<T>& mean, const bool pattern = false);
313  inline void whiten(Vector<T>& mean, const Vector<T>& mask);
315  inline void unwhiten(Vector<T>& mean, const bool pattern = false);
317  inline void sum_cols(Vector<T>& sum) const;
318 
321  inline bool isNormalized() const;
323  inline int fmax() const;
325  inline T fmaxval() const;
327  inline int fmin() const;
328 
329  // Algebric operations
332  inline void transpose(Matrix<T>& trans);
334  inline void neg();
336  inline void incrDiag();
337  inline void addDiag(const Vector<T>& diag);
338  inline void addDiag(const T diag);
339  inline void addToCols(const Vector<T>& diag);
340  inline void addVecToCols(const Vector<T>& diag, const T a = 1.0);
343  inline void svdRankOne(const Vector<T>& u0,
344  Vector<T>& u, Vector<T>& v) const;
345  inline void singularValues(Vector<T>& u) const;
346  inline void svd(Matrix<T>& U, Vector<T>& S, Matrix<T>&V) const;
347 
351  inline void eigLargestSymApprox(const Vector<T>& u0,
352  Vector<T>& u) const;
358  inline T eigLargestMagnSym(const Vector<T>& u0,
359  Vector<T>& u) const;
362  inline T eigLargestMagnSym() const;
364  inline void invSym();
366  inline void multTrans(const Vector<T>& x, Vector<T>& b,
367  const T alpha = 1.0, const T beta = 0.0) const;
369  inline void multTrans(const Vector<T>& x, Vector<T>& b,
370  const Vector<bool>& active) const;
372  inline void multTrans(const SpVector<T>& x, Vector<T>& b, const T alpha =1.0, const T beta = 0.0) const;
374  inline void mult(const Vector<T>& x, Vector<T>& b,
375  const T alpha = 1.0, const T beta = 0.0) const;
377  inline void mult(const SpVector<T>& x, Vector<T>& b,
378  const T alpha = 1.0, const T beta = 0.0) const;
380  inline void mult(const Matrix<T>& B, Matrix<T>& C,
381  const bool transA = false, const bool transB = false,
382  const T a = 1.0, const T b = 0.0) const;
384  inline void multSwitch(const Matrix<T>& B, Matrix<T>& C,
385  const bool transA = false, const bool transB = false,
386  const T a = 1.0, const T b = 0.0) const;
388  inline void mult(const SpMatrix<T>& B, Matrix<T>& C, const bool transA = false,
389  const bool transB = false, const T a = 1.0,
390  const T b = 0.0) const;
392  inline void multDiagLeft(const Vector<T>& diag);
394  inline void multDiagRight(const Vector<T>& diag);
396  inline void mult_elementWise(const Matrix<T>& B, Matrix<T>& C) const;
397  inline void div_elementWise(const Matrix<T>& B, Matrix<T>& C) const;
399  inline void XtX(Matrix<T>& XtX) const;
401  inline void XXt(Matrix<T>& XXt) const;
403  inline void upperTriXXt(Matrix<T>& XXt,
404  const int L) const;
406  inline void diag(Vector<T>& d) const;
408  inline void setDiag(const Vector<T>& d);
410  inline void setDiag(const T val);
412  inline void exp();
414  inline void Sqrt();
415  inline void Invsqrt();
417  inline T quad(const Vector<T>& vec1, const SpVector<T>& vec2) const;
419  inline void quad_mult(const Vector<T>& vec1, const SpVector<T>& vec2,
420  Vector<T>& y, const T a = 1.0, const T b = 0.0) const;
422  inline T quad(const SpVector<T>& vec) const;
424  inline void add(const Matrix<T>& mat, const T alpha = 1.0);
426  inline void add(const T alpha);
428  inline T dot(const Matrix<T>& mat) const;
430  inline void sub(const Matrix<T>& mat);
432  inline void inv_elem();
434  inline void inv() { this->inv_elem(); };
436  inline T trace() const;
438  inline T asum() const;
440  inline T normF() const;
442  inline T mean() const;
444  inline T normFsq() const;
446  inline T nrm2sq() const { return this->normFsq(); };
448  inline T norm_inf_2_col() const;
450  inline T norm_1_2_col() const;
452  inline void norm_2_cols(Vector<T>& norms) const;
454  inline void norm_2_rows(Vector<T>& norms) const;
456  inline void norm_inf_cols(Vector<T>& norms) const;
458  inline void norm_inf_rows(Vector<T>& norms) const;
460  inline void norm_l1_rows(Vector<T>& norms) const;
462  inline void norm_2sq_cols(Vector<T>& norms) const;
464  inline void norm_2sq_rows(Vector<T>& norms) const;
465  inline void thrsmax(const T nu);
466  inline void thrsmin(const T nu);
467  inline void thrsabsmin(const T nu);
469  inline void softThrshold(const T nu);
470  inline void hardThrshold(const T nu);
472  inline void thrsPos();
474  inline void rank1Update(const Vector<T>& vec1, const Vector<T>& vec2,
475  const T alpha = 1.0);
477  inline void rank1Update(const SpVector<T>& vec1, const Vector<T>& vec2,
478  const T alpha = 1.0);
480  inline void rank1Update(const Vector<T>& vec1, const SpVector<T>& vec2,
481  const T alpha = 1.0);
482  inline void rank1Update_mult(const Vector<T>& vec1, const Vector<T>& vec1b,
483  const SpVector<T>& vec2,
484  const T alpha = 1.0);
486  inline void rank1Update(const SpVector<T>& vec,
487  const T alpha = 1.0);
489  inline void rank1Update(const SpVector<T>& vec, const SpVector<T>& vec2,
490  const T alpha = 1.0);
492  inline void meanCol(Vector<T>& mean) const;
494  inline void meanRow(Vector<T>& mean) const;
496  inline void fillRow(const Vector<T>& row);
498  inline void extractRow(const int i, Vector<T>& row) const;
499  inline void setRow(const int i, const Vector<T>& row);
500  inline void addRow(const int i, const Vector<T>& row, const T a=1.0);
503  inline void conjugateGradient(const Vector<T>& b, Vector<T>& x,
504  const T tol = 1e-4, const int = 4) const;
507  inline void drop(char* fileName) const;
509  inline void NadarayaWatson(const Vector<int>& ind, const T sigma);
511  inline void blockThrshold(const T nu, const int sizeGroup);
513  inline void sparseProject(Matrix<T>& out, const T thrs, const int mode = 1, const T lambda1 = 0,
514  const T lambda2 = 0, const T lambda3 = 0, const bool pos = false, const int numThreads=-1);
515  inline void transformFilter();
516 
519  inline void toSparse(SpMatrix<T>& matrix) const;
521  inline void toSparseTrans(SpMatrix<T>& matrixTrans);
523  inline void toVect(Vector<T>& vec) const;
525  inline int V() const { return 1;};
527  inline void merge(const Matrix<T>& B, Matrix<T>& C) const;
529  inline void copyMask(Matrix<T>& out, Vector<bool>& mask) const;
530 
531  protected:
533  explicit Matrix<T>(const Matrix<T>& matrix);
535  Matrix<T>& operator=(const Matrix<T>& matrix);
536 
538  bool _externAlloc;
540  T* _X;
542  int _m;
544  int _n;
545 };
546 
548 template<typename T> class Vector {
549  friend class SpMatrix<T>;
550  friend class Matrix<T>;
551  friend class SpVector<T>;
552  public:
554  Vector();
556  Vector(int n);
558  Vector(T* X, int n);
560  explicit Vector<T>(const Vector<T>& vec);
561 
563  virtual ~Vector();
564 
567  inline void print(const char* name) const;
569  inline int max() const;
571  inline int min() const;
573  inline T maxval() const;
575  inline T minval() const;
577  inline int fmax() const;
579  inline int fmin() const;
581  inline T fmaxval() const;
583  inline T fminval() const;
585  inline T& operator[](const int index);
587  inline T operator[](const int index) const;
589  inline void copy(const Vector<T>& x);
591  inline int n() const { return _n; };
593  inline T* rawX() const { return _X; };
595  inline void fakeSize(const int n) { _n = n; };
597  inline void logspace(const int n, const T a, const T b);
598  inline int nnz() const;
599 
602  inline void setZeros();
604  inline void resize(const int n);
606  inline void setPointer(T* X, const int n);
607  inline void setData(T* X, const int n) { this->setPointer(X,n); };
609  inline void randperm(int n);
611  inline void setAleat();
613  inline void clear();
615  inline void softThrshold(const T nu);
616  inline void hardThrshold(const T nu);
618  inline void thrsmax(const T nu);
619  inline void thrsmin(const T nu);
620  inline void thrsabsmin(const T nu);
622  inline void thrshold(const T nu);
624  inline void thrsPos();
626  inline void set(const T val);
627  inline void setn(const int n) { _n = n; }; //DANGEROUS
628  inline bool alltrue() const;
629  inline bool allfalse() const;
630 
633  inline T nrm2() const;
635  inline T nrm2sq() const;
637  inline T dot(const Vector<T>& x) const;
639  inline T dot(const SpVector<T>& x) const;
641  inline void add(const Vector<T>& x, const T a = 1.0);
643  inline void add(const SpVector<T>& x, const T a = 1.0);
645  inline void add(const T a);
647  inline void sub(const Vector<T>& x);
649  inline void sub(const SpVector<T>& x);
651  inline void div(const Vector<T>& x);
653  inline void div(const Vector<T>& x, const Vector<T>& y);
655  inline void sqr(const Vector<T>& x);
657  inline void Sqrt(const Vector<T>& x);
659  inline void Sqrt();
661  inline void Invsqrt(const Vector<T>& x);
663  inline void Invsqrt();
665  inline void inv(const Vector<T>& x);
667  inline void inv();
669  inline void mult(const Vector<T>& x, const Vector<T>& y);
670  inline void mult_elementWise(const Vector<T>& B, Vector<T>& C) const { C.mult(*this,B); };
672  inline void normalize();
674  inline void normalize2();
676  inline void whiten(Vector<T>& mean, const bool pattern = false);
678  inline void whiten(Vector<T>& mean, const
679  Vector<T>& mask);
681  inline void whiten(const int V);
683  inline T mean();
685  inline T std();
687  inline T KL(const Vector<T>& X);
689  inline void unwhiten(Vector<T>& mean, const bool pattern = false);
691  inline void scal(const T a);
693  inline void neg();
695  inline void exp();
697  inline void log();
699  inline void logexp();
701  inline T softmax(const int y);
703  inline T asum() const;
704  inline T lzero() const;
706  inline T afused() const;
708  inline T sum() const;
710  inline void sign(Vector<T>& signs) const;
713  inline void l1project(Vector<T>& out, const T thrs, const bool simplex = false) const;
714  inline void l1project_weighted(Vector<T>& out, const Vector<T>& weights, const T thrs, const bool residual = false) const;
715  inline void l1l2projectb(Vector<T>& out, const T thrs, const T gamma, const bool pos = false,
716  const int mode = 1);
717  inline void sparseProject(Vector<T>& out, const T thrs, const int mode = 1, const T lambda1 = 0,
718  const T lambda2 = 0, const T lambda3 = 0, const bool pos = false);
719  inline void project_sft(const Vector<int>& labels, const int clas);
720  inline void project_sft_binary(const Vector<T>& labels);
724  inline void l1l2project(Vector<T>& out, const T thrs, const T gamma, const bool pos = false) const;
725  inline void fusedProject(Vector<T>& out, const T lambda1, const T lambda2, const int itermax);
726  inline void fusedProjectHomotopy(Vector<T>& out, const T lambda1,const T lambda2,const T lambda3 = 0,
727  const bool penalty = true);
730  inline void sort(Vector<T>& out, const bool mode) const;
732  inline void sort(const bool mode);
734  inline void sort2(Vector<T>& out, Vector<int>& key, const bool mode) const;
736  inline void sort2(Vector<int>& key, const bool mode);
738  inline void applyBayerPattern(const int offset);
739 
740 
743  inline void toSparse(SpVector<T>& vec) const;
745  inline void copyMask(Vector<T>& out, Vector<bool>& mask) const;
746 
747  private:
749  Vector<T>& operator=(const Vector<T>& vec);
750 
752  bool _externAlloc;
754  T* _X;
756  int _n;
757 };
758 
759 
761 template<typename T> class SpMatrix : public Data<T>, public AbstractMatrixB<T> {
762  friend class Matrix<T>;
763  friend class SpVector<T>;
764  public:
766  SpMatrix(T* v, int* r, int* pB, int* pE, int m, int n, int nzmax);
768  SpMatrix(int m, int n, int nzmax);
770  SpMatrix();
771 
773  ~SpMatrix();
774 
777  inline void refCol(int i, SpVector<T>& vec) const;
779  inline int pB(const int i) const { return _pB[i]; };
781  inline int r(const int i) const { return _r[i]; };
783  inline T v(const int i) const { return _v[i]; };
785  inline int nzmax() const { return _nzmax; };
787  inline int n() const { return _n; };
789  inline int m() const { return _m; };
791  inline int V() const { return 1; };
793  inline T operator[](const int index) const;
794  void getData(Vector<T>& data, const int index) const;
795  void getGroup(Matrix<T>& data, const vector_groups& groups,
796  const int i) const ;
798  inline void print(const string& name) const;
800  inline T asum() const;
802  inline T normFsq() const;
804  inline int* pB() const { return _pB; };
806  inline int* pE() const { return _pE; };
808  inline int* r() const { return _r; };
810  inline T* v() const { return _v; };
812  inline int nnz() const { return _pB[_n]; };
813  inline void add_direct(const SpMatrix<T>& mat, const T a);
814  inline void copy_direct(const SpMatrix<T>& mat);
815  inline T dot_direct(const SpMatrix<T>& mat) const;
816 
819  inline void clear();
821  inline void resize(const int m, const int n, const int nzmax);
823  inline void scal(const T a) const;
824 
827  inline void AAt(Matrix<T>& aat) const;
829  inline void AAt(Matrix<T>& aat, const Vector<int>& indices) const;
831  inline void wAAt(const Vector<T>& w, Matrix<T>& aat) const;
833  inline void XAt(const Matrix<T>& X, Matrix<T>& XAt) const;
835  inline void XAt(const Matrix<T>& X, Matrix<T>& XAt,
836  const Vector<int>& indices) const;
838  inline void wXAt( const Vector<T>& w, const Matrix<T>& X,
839  Matrix<T>& XAt, const int numthreads=-1) const;
840  inline void XtX(Matrix<T>& XtX) const;
841 
843  inline void multTrans(const Vector<T>& x, Vector<T>& y,
844  const T alpha = 1.0, const T beta = 0.0) const;
845  inline void multTrans(const SpVector<T>& x, Vector<T>& y,
846  const T alpha = 1.0, const T beta = 0.0) const;
848  inline void mult(const SpVector<T>& x, Vector<T>& b,
849  const T alpha = 1.0, const T beta = 0.0) const;
851  inline void mult(const Vector<T>& x, Vector<T>& b,
852  const T alpha = 1.0, const T beta = 0.0) const;
854  inline void mult(const Matrix<T>& B, Matrix<T>& C,
855  const bool transA = false, const bool transB = false,
856  const T a = 1.0, const T b = 0.0) const;
858  inline void multSwitch(const Matrix<T>& B, Matrix<T>& C,
859  const bool transA = false, const bool transB = false,
860  const T a = 1.0, const T b = 0.0) const;
862  inline void mult(const SpMatrix<T>& B, Matrix<T>& C, const bool transA = false,
863  const bool transB = false, const T a = 1.0,
864  const T b = 0.0) const;
866  inline void copyTo(Matrix<T>& mat) const { this->toFull(mat); };
868  inline T dot(const Matrix<T>& x) const;
869  inline void copyRow(const int i, Vector<T>& x) const;
870  inline void sum_cols(Vector<T>& sum) const;
871  inline void copy(const SpMatrix<T>& mat);
872 
875  inline void toFull(Matrix<T>& matrix) const;
877  inline void toFullTrans(Matrix<T>& matrix) const;
878 
880  inline void convert(const Matrix<T>&v, const Matrix<int>& r,
881  const int K);
883  inline void convert2(const Matrix<T>&v, const Vector<int>& r,
884  const int K);
886  inline void norm_2sq_cols(Vector<T>& norms) const;
888  inline void norm_0_cols(Vector<T>& norms) const;
890  inline void norm_1_cols(Vector<T>& norms) const;
891  inline void addVecToCols(const Vector<T>& diag, const T a = 1.0);
892  inline void addVecToColsWeighted(const Vector<T>& diag, const T* weights, const T a = 1.0);
893 
894  private:
896  explicit SpMatrix(const SpMatrix<T>& matrix);
897  SpMatrix<T>& operator=(const SpMatrix<T>& matrix);
898 
900  bool _externAlloc;
902  T* _v;
904  int* _r;
906  int* _pB;
908  int* _pE;
910  int _m;
912  int _n;
914  int _nzmax;
915 };
916 
918 template <typename T> class SpVector {
919  friend class Matrix<T>;
920  friend class SpMatrix<T>;
921  friend class Vector<T>;
922  public:
924  SpVector(T* v, int* r, int L, int nzmax);
926  SpVector(int nzmax);
928  SpVector();
929 
931  ~SpVector();
932 
935  inline T nzmax() const { return _nzmax; };
937  inline T length() const { return _L; };
939  inline T asum() const;
941  inline T nrm2sq() const;
943  inline T nrm2() const;
945  inline T fmaxval() const;
947  inline void print(const string& name) const;
949  inline void refIndices(Vector<int>& indices) const;
951  inline void refVal(Vector<T>& val) const;
953  inline int r(const int i) const { return _r[i]; };
955  inline T v(const int i) const { return _v[i]; };
956  inline T* rawX() const { return _v; };
958  inline int L() const { return _L; };
960  inline void setL(const int L) { _L=L; };
962  inline void sqr();
964  inline T dot(const SpVector<T>& vec) const;
965 
968  inline void clear();
970  inline void resize(const int nzmax);
971 
973  void inline toSpMatrix(SpMatrix<T>& out,
974  const int m, const int n) const;
976  void inline toFull(Vector<T>& out) const;
977 
978 
979  private:
981  explicit SpVector(const SpVector<T>& vector);
982  SpVector<T>& operator=(const SpVector<T>& vector);
983 
985  bool _externAlloc;
987  T* _v;
989  int* _r;
991  int _L;
993  int _nzmax;
994 };
995 
996 
998 template<typename T> class ProdMatrix : public AbstractMatrix<T> {
999  public:
1000  ProdMatrix();
1002  ProdMatrix(const Matrix<T>& D, const bool high_memory = true);
1004  ProdMatrix(const Matrix<T>& D, const Matrix<T>& X, const bool high_memory = true);
1006  /*ProdMatrix(const SpMatrix<T>& D, const Matrix<T>& X,
1007  const bool transD = false, const bool transX = false);*/
1008 
1010  ~ProdMatrix() { delete(_DtX);} ;
1011 
1013  inline void setMatrices(const Matrix<T>& D, const bool high_memory=true);
1014  inline void setMatrices(const Matrix<T>& D,
1015  const Matrix<T>& X, const bool high_memory=true);
1017  inline void copyCol(const int i, Vector<T>& DtXi) const;
1019  inline void extract_rawCol(const int i,T* DtXi) const;
1021  virtual void add_rawCol(const int i, T* DtXi, const T a) const;
1023  void inline addDiag(const T diag);
1025  void inline diag(Vector<T>& diag) const;
1027  inline int n() const { return _n;};
1029  inline int m() const { return _m;};
1031  inline T operator()(const int index1, const int index2) const;
1033  inline T operator[](const int index) const;
1034 
1035  private:
1037  Matrix<T>* _DtX;
1038  const Matrix<T>* _X;
1039  const Matrix<T>* _D;
1040  bool _high_memory;
1041  int _n;
1042  int _m;
1043  T _addDiag;
1044 };
1045 
1046 
1047 /* ************************************
1048  * Implementation of the class Matrix
1049  * ************************************/
1050 
1052 template <typename T> Matrix<T>::Matrix(T* X, int m, int n) :
1053  _externAlloc(true), _X(X), _m(m), _n(n) { };
1054 
1055 
1057 template <typename T> Matrix<T>::Matrix(int m, int n) :
1058  _externAlloc(false), _m(m), _n(n) {
1059 #pragma omp critical
1060  {
1061  _X= new T[_n*_m];
1062  }
1063  };
1064 
1066 template <typename T> Matrix<T>::Matrix() :
1067  _externAlloc(false), _X(NULL), _m(0), _n(0) { };
1068 
1070 template <typename T> Matrix<T>::~Matrix() {
1071  clear();
1072 };
1073 
1075 template <typename T> inline T& Matrix<T>::operator()(const int i, const int j) {
1076  return _X[j*_m+i];
1077 };
1078 
1080 template <typename T> inline T Matrix<T>::operator()(const int i, const int j) const {
1081  return _X[j*_m+i];
1082 };
1083 
1085 template <typename T> inline void Matrix<T>::print(const string& name) const {
1086  std::cerr << name << std::endl;
1087  std::cerr << _m << " x " << _n << std::endl;
1088  for (int i = 0; i<_m; ++i) {
1089  for (int j = 0; j<_n; ++j) {
1090  printf("%10.5g ",static_cast<double>(_X[j*_m+i]));
1091  // std::cerr << _X[j*_m+i] << " ";
1092  }
1093  printf("\n ");
1094  //std::cerr << std::endl;
1095  }
1096  printf("\n ");
1097 };
1098 
1100 template <typename T> inline void Matrix<T>::copyCol(const int i, Vector<T>& x) const {
1101  assert(i >= 0 && i<_n);
1102  x.resize(_m);
1103  cblas_copy<T>(_m,_X+i*_m,1,x._X,1);
1104 };
1105 
1107 template <typename T> inline void Matrix<T>::copyRow(const int i, Vector<T>& x) const {
1108  assert(i >= 0 && i<_m);
1109  x.resize(_n);
1110  cblas_copy<T>(_n,_X+i,_m,x._X,1);
1111 };
1112 
1113 
1115 template <typename T> inline void Matrix<T>::extract_rawCol(const int i, T* x) const {
1116  assert(i >= 0 && i<_n);
1117  cblas_copy<T>(_m,_X+i*_m,1,x,1);
1118 };
1119 
1121 template <typename T> inline void Matrix<T>::add_rawCol(const int i, T* x, const T a) const {
1122  assert(i >= 0 && i<_n);
1123  cblas_axpy<T>(_m,a,_X+i*_m,1,x,1);
1124 };
1125 
1127 template <typename T> inline void Matrix<T>::getData(Vector<T>& x, const int i) const {
1128  this->copyCol(i,x);
1129 };
1130 
1131 template <typename T> inline void Matrix<T>::getGroup(Matrix<T>& data,
1132  const vector_groups& groups, const int i) const {
1133  const group& gr = groups[i];
1134  const int N = gr.size();
1135  data.resize(_m,N);
1136  int count=0;
1137  for (group::const_iterator it = gr.begin(); it != gr.end(); ++it) {
1138  cblas_copy<T>(_m,_X+(*it)*_m,1,data._X+count*_m,1);
1139  ++count;
1140  }
1141 };
1142 
1144 template <typename T> inline void Matrix<T>::refCol(int i, Vector<T>& x) const {
1145  assert(i >= 0 && i<_n);
1146  x.clear();
1147  x._X=_X+i*_m;
1148  x._n=_m;
1149  x._externAlloc=true;
1150 };
1151 
1153 template <typename T> inline void Matrix<T>::refSubMat(int i, int n, Matrix<T>& mat) const {
1154  mat.setData(_X+i*_m,_m,n);
1155 }
1156 
1158 template <typename T> inline bool Matrix<T>::isNormalized() const {
1159  for (int i = 0; i<_n; ++i) {
1160  T norm=cblas_nrm2<T>(_m,_X+_m*i,1);
1161  if (fabs(norm - 1.0) > 1e-6) return false;
1162  }
1163  return true;
1164 };
1165 
1167 template <typename T>
1168 inline void Matrix<T>::clean() {
1169  this->normalize();
1170  Matrix<T> G;
1171  this->XtX(G);
1172  T* prG = G._X;
1174  for (int i = 0; i<_n; ++i) {
1175  for (int j = i+1; j<_n; ++j) {
1176  if (prG[i*_n+j] > 0.99) {
1177  // remove nasty column j and put random values inside
1178  Vector<T> col;
1179  this->refCol(j,col);
1180  col.setAleat();
1181  col.normalize();
1182  }
1183  }
1184  }
1185 };
1186 
1188 template <typename T> inline int Matrix<T>::fmax() const {
1189  return cblas_iamax<T>(_n*_m,_X,1);
1190 };
1191 
1193 template <typename T> inline T Matrix<T>::fmaxval() const {
1194  return _X[cblas_iamax<T>(_n*_m,_X,1)];
1195 };
1196 
1197 
1199 template <typename T> inline int Matrix<T>::fmin() const {
1200  return cblas_iamin<T>(_n*_m,_X,1);
1201 };
1202 
1204 template <typename T> inline void Matrix<T>::subMatrixSym(
1205  const Vector<int>& indices, Matrix<T>& subMatrix) const {
1206  int L = indices.n();
1207  subMatrix.resize(L,L);
1208  T* out = subMatrix._X;
1209  int* rawInd = indices.rawX();
1210  for (int i = 0; i<L; ++i)
1211  for (int j = 0; j<=i; ++j)
1212  out[i*L+j]=_X[rawInd[i]*_n+rawInd[j]];
1213  subMatrix.fillSymmetric();
1214 };
1215 
1217 template <typename T> inline void Matrix<T>::resize(int m, int n) {
1218  if (_n==n && _m==m) return;
1219  clear();
1220  _n=n;
1221  _m=m;
1222  _externAlloc=false;
1223 #pragma omp critical
1224  {
1225  _X=new T[_n*_m];
1226  }
1227  setZeros();
1228 };
1229 
1231 template <typename T> inline void Matrix<T>::setData(T* X, int m, int n) {
1232  clear();
1233  _X=X;
1234  _m=m;
1235  _n=n;
1236  _externAlloc=true;
1237 };
1238 
1240 template <typename T> inline void Matrix<T>::setZeros() {
1241  memset(_X,0,_n*_m*sizeof(T));
1242 };
1243 
1245 template <typename T> inline void Matrix<T>::set(const T a) {
1246  for (int i = 0; i<_n*_m; ++i) _X[i]=a;
1247 };
1248 
1250 template <typename T> inline void Matrix<T>::clear() {
1251  if (!_externAlloc) delete[](_X);
1252  _n=0;
1253  _m=0;
1254  _X=NULL;
1255  _externAlloc=true;
1256 };
1257 
1259 template <typename T> inline void Matrix<T>::setAleat() {
1260  for (int i = 0; i<_n*_m; ++i) _X[i]=normalDistrib<T>();
1261 };
1262 
1264 template <typename T> inline void Matrix<T>::eye() {
1265  this->setZeros();
1266  for (int i = 0; i<MIN(_n,_m); ++i) _X[i*_m+i] = T(1.0);
1267 };
1268 
1270 template <typename T> inline void Matrix<T>::normalize() {
1271  //T constant = 1.0/sqrt(_m);
1272  for (int i = 0; i<_n; ++i) {
1273  T norm=cblas_nrm2<T>(_m,_X+_m*i,1);
1274  if (norm > 1e-10) {
1275  T invNorm=1.0/norm;
1276  cblas_scal<T>(_m,invNorm,_X+_m*i,1);
1277  } else {
1278  // for (int j = 0; j<_m; ++j) _X[_m*i+j]=constant;
1279  Vector<T> d;
1280  this->refCol(i,d);
1281  d.setAleat();
1282  d.normalize();
1283  }
1284  }
1285 };
1286 
1288 template <typename T> inline void Matrix<T>::normalize2() {
1289  for (int i = 0; i<_n; ++i) {
1290  T norm=cblas_nrm2<T>(_m,_X+_m*i,1);
1291  if (norm > 1.0) {
1292  T invNorm=1.0/norm;
1293  cblas_scal<T>(_m,invNorm,_X+_m*i,1);
1294  }
1295  }
1296 };
1297 
1299 template <typename T> inline void Matrix<T>::center() {
1300  for (int i = 0; i<_n; ++i) {
1301  Vector<T> col;
1302  this->refCol(i,col);
1303  T sum = col.sum();
1304  col.add(-sum/static_cast<T>(_m));
1305  }
1306 };
1307 
1309 template <typename T> inline void Matrix<T>::center_rows() {
1310  Vector<T> mean_rows(_m);
1311  mean_rows.setZeros();
1312  for (int i = 0; i<_n; ++i)
1313  for (int j = 0; j<_m; ++j)
1314  mean_rows[j] += _X[i*_m+j];
1315  mean_rows.scal(T(1.0)/_n);
1316  for (int i = 0; i<_n; ++i)
1317  for (int j = 0; j<_m; ++j)
1318  _X[i*_m+j] -= mean_rows[j];
1319 };
1320 
1322 template <typename T> inline void Matrix<T>::center(Vector<T>& centers) {
1323  centers.resize(_n);
1324  for (int i = 0; i<_n; ++i) {
1325  Vector<T> col;
1326  this->refCol(i,col);
1327  T sum = col.sum()/static_cast<T>(_m);
1328  centers[i]=sum;
1329  col.add(-sum);
1330  }
1331 };
1332 
1334 template <typename T> inline void Matrix<T>::scal(const T a) {
1335  cblas_scal<T>(_n*_m,a,_X,1);
1336 };
1337 
1339 template <typename T> inline void Matrix<T>::copy(const Matrix<T>& mat) {
1340  resize(mat._m,mat._n);
1341  cblas_copy<T>(_m*_n,mat._X,1,_X,1);
1342 };
1343 
1345 template <typename T> inline void Matrix<T>::copyRef(const Matrix<T>& mat) {
1346  this->setData(mat.rawX(),mat.m(),mat.n());
1347 };
1348 
1351 template <typename T> inline void Matrix<T>::fillSymmetric() {
1352  for (int i = 0; i<_n; ++i) {
1353  for (int j =0; j<i; ++j) {
1354  _X[j*_m+i]=_X[i*_m+j];
1355  }
1356  }
1357 };
1358 template <typename T> inline void Matrix<T>::fillSymmetric2() {
1359  for (int i = 0; i<_n; ++i) {
1360  for (int j =0; j<i; ++j) {
1361  _X[i*_m+j]=_X[j*_m+i];
1362  }
1363  }
1364 };
1365 
1366 
1367 template <typename T> inline void Matrix<T>::whiten(const int V) {
1368  const int sizePatch=_m/V;
1369  for (int i = 0; i<_n; ++i) {
1370  for (int j = 0; j<V; ++j) {
1371  T mean = 0;
1372  for (int k = 0; k<sizePatch; ++k) {
1373  mean+=_X[i*_m+sizePatch*j+k];
1374  }
1375  mean /= sizePatch;
1376  for (int k = 0; k<sizePatch; ++k) {
1377  _X[i*_m+sizePatch*j+k]-=mean;
1378  }
1379  }
1380  }
1381 };
1382 
1383 template <typename T> inline void Matrix<T>::whiten(Vector<T>& mean, const bool pattern) {
1384  mean.setZeros();
1385  if (pattern) {
1386  const int n =static_cast<int>(sqrt(static_cast<T>(_m)));
1387  int count[4];
1388  for (int i = 0; i<4; ++i) count[i]=0;
1389  for (int i = 0; i<_n; ++i) {
1390  int offsetx=0;
1391  for (int j = 0; j<n; ++j) {
1392  offsetx= (offsetx+1) % 2;
1393  int offsety=0;
1394  for (int k = 0; k<n; ++k) {
1395  offsety= (offsety+1) % 2;
1396  mean[2*offsetx+offsety]+=_X[i*_m+j*n+k];
1397  count[2*offsetx+offsety]++;
1398  }
1399  }
1400  }
1401  for (int i = 0; i<4; ++i)
1402  mean[i] /= count[i];
1403  for (int i = 0; i<_n; ++i) {
1404  int offsetx=0;
1405  for (int j = 0; j<n; ++j) {
1406  offsetx= (offsetx+1) % 2;
1407  int offsety=0;
1408  for (int k = 0; k<n; ++k) {
1409  offsety= (offsety+1) % 2;
1410  _X[i*_m+j*n+k]-=mean[2*offsetx+offsety];
1411  }
1412  }
1413  }
1414  } else {
1415  const int V = mean.n();
1416  const int sizePatch=_m/V;
1417  for (int i = 0; i<_n; ++i) {
1418  for (int j = 0; j<V; ++j) {
1419  for (int k = 0; k<sizePatch; ++k) {
1420  mean[j]+=_X[i*_m+sizePatch*j+k];
1421  }
1422  }
1423  }
1424  mean.scal(T(1.0)/(_n*sizePatch));
1425  for (int i = 0; i<_n; ++i) {
1426  for (int j = 0; j<V; ++j) {
1427  for (int k = 0; k<sizePatch; ++k) {
1428  _X[i*_m+sizePatch*j+k]-=mean[j];
1429  }
1430  }
1431  }
1432  }
1433 };
1434 
1435 template <typename T> inline void Matrix<T>::whiten(Vector<T>& mean, const
1436  Vector<T>& mask) {
1437  const int V = mean.n();
1438  const int sizePatch=_m/V;
1439  mean.setZeros();
1440  for (int i = 0; i<_n; ++i) {
1441  for (int j = 0; j<V; ++j) {
1442  for (int k = 0; k<sizePatch; ++k) {
1443  mean[j]+=_X[i*_m+sizePatch*j+k];
1444  }
1445  }
1446  }
1447  for (int i = 0; i<V; ++i)
1448  mean[i] /= _n*cblas_asum(sizePatch,mask._X+i*sizePatch,1);
1449  for (int i = 0; i<_n; ++i) {
1450  for (int j = 0; j<V; ++j) {
1451  for (int k = 0; k<sizePatch; ++k) {
1452  if (mask[sizePatch*j+k])
1453  _X[i*_m+sizePatch*j+k]-=mean[j];
1454  }
1455  }
1456  }
1457 };
1458 
1459 
1460 template <typename T> inline void Matrix<T>::unwhiten(Vector<T>& mean, const bool pattern) {
1461  if (pattern) {
1462  const int n =static_cast<int>(sqrt(static_cast<T>(_m)));
1463  for (int i = 0; i<_n; ++i) {
1464  int offsetx=0;
1465  for (int j = 0; j<n; ++j) {
1466  offsetx= (offsetx+1) % 2;
1467  int offsety=0;
1468  for (int k = 0; k<n; ++k) {
1469  offsety= (offsety+1) % 2;
1470  _X[i*_m+j*n+k]+=mean[2*offsetx+offsety];
1471  }
1472  }
1473  }
1474  } else {
1475  const int V = mean.n();
1476  const int sizePatch=_m/V;
1477  for (int i = 0; i<_n; ++i) {
1478  for (int j = 0; j<V; ++j) {
1479  for (int k = 0; k<sizePatch; ++k) {
1480  _X[i*_m+sizePatch*j+k]+=mean[j];
1481  }
1482  }
1483  }
1484  }
1485 };
1486 
1489 template <typename T> inline void Matrix<T>::transpose(Matrix<T>& trans) {
1490  trans.resize(_n,_m);
1491  T* out = trans._X;
1492  for (int i = 0; i<_n; ++i)
1493  for (int j = 0; j<_m; ++j)
1494  out[j*_n+i] = _X[i*_m+j];
1495 };
1496 
1498 template <typename T> inline void Matrix<T>::neg() {
1499  for (int i = 0; i<_n*_m; ++i) _X[i]=-_X[i];
1500 };
1501 
1502 template <typename T> inline void Matrix<T>::incrDiag() {
1503  for (int i = 0; i<MIN(_n,_m); ++i) ++_X[i*_m+i];
1504 };
1505 
1506 template <typename T> inline void Matrix<T>::addDiag(
1507  const Vector<T>& diag) {
1508  T* d= diag.rawX();
1509  for (int i = 0; i<MIN(_n,_m); ++i) _X[i*_m+i] += d[i];
1510 };
1511 
1512 template <typename T> inline void Matrix<T>::addDiag(
1513  const T diag) {
1514  for (int i = 0; i<MIN(_n,_m); ++i) _X[i*_m+i] += diag;
1515 };
1516 
1517 template <typename T> inline void Matrix<T>::addToCols(
1518  const Vector<T>& cent) {
1519  Vector<T> col;
1520  for (int i = 0; i<_n; ++i) {
1521  this->refCol(i,col);
1522  col.add(cent[i]);
1523  }
1524 };
1525 
1526 template <typename T> inline void Matrix<T>::addVecToCols(
1527  const Vector<T>& vec, const T a) {
1528  Vector<T> col;
1529  for (int i = 0; i<_n; ++i) {
1530  this->refCol(i,col);
1531  col.add(vec,a);
1532  }
1533 };
1534 
1537 template <typename T> inline void Matrix<T>::svdRankOne(const Vector<T>& u0,
1538  Vector<T>& u, Vector<T>& v) const {
1539  int i;
1540  const int max_iter=MAX(_m,MAX(_n,200));
1541  const T eps=1e-10;
1542  u.resize(_m);
1543  v.resize(_n);
1544  T norm=u0.nrm2();
1545  Vector<T> up(u0);
1546  if (norm < EPSILON) up.setAleat();
1547  up.normalize();
1548  multTrans(up,v);
1549  for (i = 0; i<max_iter; ++i) {
1550  mult(v,u);
1551  norm=u.nrm2();
1552  u.scal(1.0/norm);
1553  multTrans(u,v);
1554  T theta=u.dot(up);
1555  if (i > 10 && (1 - fabs(theta)) < eps) break;
1556  up.copy(u);
1557  }
1558 };
1559 
1560 template <typename T> inline void Matrix<T>::singularValues(Vector<T>& u) const {
1561  u.resize(MIN(_m,_n));
1562  if (_m > 10*_n) {
1563  Matrix<T> XtX;
1564  this->XtX(XtX);
1565  syev<T>(no,lower,_n,XtX.rawX(),_n,u.rawX());
1566  u.thrsPos();
1567  u.Sqrt();
1568  } else if (_n > 10*_m) {
1569  Matrix<T> XXt;
1570  this->XXt(XXt);
1571  syev<T>(no,lower,_m,XXt.rawX(),_m,u.rawX());
1572  u.thrsPos();
1573  u.Sqrt();
1574  } else {
1575  T* vu, *vv;
1576  Matrix<T> copyX;
1577  copyX.copy(*this);
1578  gesvd<T>(no,no,_m,_n,copyX._X,_m,u.rawX(),vu,1,vv,1);
1579  }
1580 };
1581 
1582 template <typename T> inline void Matrix<T>::svd(Matrix<T>& U, Vector<T>& S, Matrix<T>&V) const {
1583  const int num_eig=MIN(_m,_n);
1584  S.resize(num_eig);
1585  U.resize(_m,num_eig);
1586  V.resize(num_eig,_n);
1587  if (_m > 10*_n) {
1588  Matrix<T> Vt(_n,_n);
1589  this->XtX(Vt);
1590  syev<T>(allV,lower,_n,Vt.rawX(),_n,S.rawX());
1591  S.thrsPos();
1592  S.Sqrt();
1593  this->mult(Vt,U);
1594  Vt.transpose(V);
1595  Vector<T> inveigs;
1596  inveigs.copy(S);
1597  for (int i = 0; i<num_eig; ++i)
1598  if (S[i] > 1e-10) {
1599  inveigs[i]=T(1.0)/S[i];
1600  } else {
1601  inveigs[i]=T(1.0);
1602  }
1603  U.multDiagRight(inveigs);
1604  } else if (_n > 10*_m) {
1605  this->XXt(U);
1606  syev<T>(allV,lower,_m,U.rawX(),_m,S.rawX());
1607  S.thrsPos();
1608  S.Sqrt();
1609  U.mult(*this,V,true,false);
1610  Vector<T> inveigs;
1611  inveigs.copy(S);
1612  for (int i = 0; i<num_eig; ++i)
1613  if (S[i] > 1e-10) {
1614  inveigs[i]=T(1.0)/S[i];
1615  } else {
1616  inveigs[i]=T(1.0);
1617  }
1618  V.multDiagLeft(inveigs);
1619  } else {
1620  Matrix<T> copyX;
1621  copyX.copy(*this);
1622  gesvd<T>(reduced,reduced,_m,_n,copyX._X,_m,S.rawX(),U.rawX(),_m,V.rawX(),num_eig);
1623  }
1624 };
1625 
1629 template <typename T> inline void Matrix<T>::eigLargestSymApprox(
1630  const Vector<T>& u0, Vector<T>& u) const {
1631  int i,j;
1632  const int max_iter=100;
1633  const T eps=10e-6;
1634  u.copy(u0);
1635  T norm = u.nrm2();
1636  T theta;
1637  u.scal(1.0/norm);
1638  Vector<T> up(u);
1639  Vector<T> uor(u);
1640  T lambda=T();
1641 
1642  for (j = 0; j<2;++j) {
1643  up.copy(u);
1644  for (i = 0; i<max_iter; ++i) {
1645  mult(up,u);
1646  norm = u.nrm2();
1647  u.scal(1.0/norm);
1648  theta=u.dot(up);
1649  if ((1 - fabs(theta)) < eps) break;
1650  up.copy(u);
1651  }
1652  lambda+=theta*norm;
1653  if (isnan(lambda)) {
1654  std::cerr << "eigLargestSymApprox failed" << std::endl;
1655  exit(1);
1656  }
1657  if (j == 1 && lambda < eps) {
1658  u.copy(uor);
1659  break;
1660  }
1661  if (theta >= 0) break;
1662  u.copy(uor);
1663  for (i = 0; i<_m; ++i) _X[i*_m+i]-=lambda;
1664  }
1665 };
1666 
1672 template <typename T> inline T Matrix<T>::eigLargestMagnSym(
1673  const Vector<T>& u0, Vector<T>& u) const {
1674  const int max_iter=1000;
1675  const T eps=10e-6;
1676  u.copy(u0);
1677  T norm = u.nrm2();
1678  u.scal(1.0/norm);
1679  Vector<T> up(u);
1680  T lambda=T();
1681 
1682  for (int i = 0; i<max_iter; ++i) {
1683  mult(u,up);
1684  u.copy(up);
1685  norm=u.nrm2();
1686  if (norm > 0) u.scal(1.0/norm);
1687  if (norm == 0 || fabs(norm-lambda)/norm < eps) break;
1688  lambda=norm;
1689  }
1690  return norm;
1691 };
1692 
1695 template <typename T> inline T Matrix<T>::eigLargestMagnSym() const {
1696  const int max_iter=1000;
1697  const T eps=10e-6;
1698  Vector<T> u(_m);
1699  u.setAleat();
1700  T norm = u.nrm2();
1701  u.scal(1.0/norm);
1702  Vector<T> up(u);
1703  T lambda=T();
1704  for (int i = 0; i<max_iter; ++i) {
1705  mult(u,up);
1706  u.copy(up);
1707  norm=u.nrm2();
1708  if (fabs(norm-lambda) < eps) break;
1709  lambda=norm;
1710  u.scal(1.0/norm);
1711  }
1712  return norm;
1713 };
1714 
1716 template <typename T> inline void Matrix<T>::invSym() {
1717  // int lwork=2*_n;
1718  // T* work;
1719 //#ifdef USE_BLAS_LIB
1720 // INTT* ipiv;
1721 //#else
1722 // int* ipiv;
1723 //#endif
1724 //#pragma omp critical
1725 // {
1726 // work= new T[lwork];
1727 //#ifdef USE_BLAS_LIB
1729 //#else
1730 // ipiv= new int[lwork];
1731 //#endif
1732 // }
1733 // sytrf<T>(upper,_n,_X,_n,ipiv,work,lwork);
1734 // sytri<T>(upper,_n,_X,_n,ipiv,work);
1735 // sytrf<T>(upper,_n,_X,_n);
1736  sytri<T>(upper,_n,_X,_n);
1737  this->fillSymmetric();
1738 // delete[](work);
1739 // delete[](ipiv);
1740 };
1741 
1743 template <typename T> inline void Matrix<T>::multTrans(const Vector<T>& x,
1744  Vector<T>& b, const T a, const T c) const {
1745  b.resize(_n);
1746  // assert(x._n == _m && b._n == _n);
1747  cblas_gemv<T>(CblasColMajor,CblasTrans,_m,_n,a,_X,_m,x._X,1,c,b._X,1);
1748 };
1749 
1751 template <typename T> inline void Matrix<T>::multTrans(const SpVector<T>& x,
1752  Vector<T>& b, const T alpha, const T beta) const {
1753  b.resize(_n);
1754  Vector<T> col;
1755  if (beta) {
1756  for (int i = 0; i<_n; ++i) {
1757  refCol(i,col);
1758  b._X[i] = alpha*col.dot(x);
1759  }
1760  } else {
1761 
1762  for (int i = 0; i<_n; ++i) {
1763  refCol(i,col);
1764  b._X[i] = beta*b._X[i]+alpha*col.dot(x);
1765  }
1766  }
1767 };
1768 
1769 template <typename T> inline void Matrix<T>::multTrans(
1770  const Vector<T>& x, Vector<T>& b, const Vector<bool>& active) const {
1771  b.setZeros();
1772  Vector<T> col;
1773  bool* pr_active=active.rawX();
1774  for (int i = 0; i<_n; ++i) {
1775  if (pr_active[i]) {
1776  this->refCol(i,col);
1777  b._X[i]=col.dot(x);
1778  }
1779  }
1780 };
1781 
1783 template <typename T> inline void Matrix<T>::mult(const Vector<T>& x,
1784  Vector<T>& b, const T a, const T c) const {
1785  // assert(x._n == _n && b._n == _m);
1786  b.resize(_m);
1787  cblas_gemv<T>(CblasColMajor,CblasNoTrans,_m,_n,a,_X,_m,x._X,1,c,b._X,1);
1788 };
1789 
1791 template <typename T> inline void Matrix<T>::mult(const SpVector<T>& x,
1792  Vector<T>& b, const T a, const T a2) const {
1793  if (!a2) {
1794  b.setZeros();
1795  } else if (a2 != 1.0) {
1796  b.scal(a2);
1797  }
1798  if (a == 1.0) {
1799  for (int i = 0; i<x._L; ++i) {
1800  cblas_axpy<T>(_m,x._v[i],_X+x._r[i]*_m,1,b._X,1);
1801  }
1802  } else {
1803  for (int i = 0; i<x._L; ++i) {
1804  cblas_axpy<T>(_m,a*x._v[i],_X+x._r[i]*_m,1,b._X,1);
1805  }
1806  }
1807 };
1808 
1810 template <typename T> inline void Matrix<T>::mult(const Matrix<T>& B,
1811  Matrix<T>& C, const bool transA, const bool transB,
1812  const T a, const T b) const {
1813  CBLAS_TRANSPOSE trA,trB;
1814  int m,k,n;
1815  if (transA) {
1816  trA = CblasTrans;
1817  m = _n;
1818  k = _m;
1819  } else {
1820  trA= CblasNoTrans;
1821  m = _m;
1822  k = _n;
1823  }
1824  if (transB) {
1825  trB = CblasTrans;
1826  n = B._m;
1827  // assert(B._n == k);
1828  } else {
1829  trB = CblasNoTrans;
1830  n = B._n;
1831  // assert(B._m == k);
1832  }
1833  C.resize(m,n);
1834  cblas_gemm<T>(CblasColMajor,trA,trB,m,n,k,a,_X,_m,B._X,B._m,
1835  b,C._X,C._m);
1836 };
1837 
1839 template <typename T>
1840 inline void Matrix<T>::multSwitch(const Matrix<T>& B, Matrix<T>& C,
1841  const bool transA, const bool transB,
1842  const T a, const T b) const {
1843  B.mult(*this,C,transB,transA,a,b);
1844 };
1845 
1847 template <typename T>
1848 inline void Matrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
1849  const bool transA, const bool transB,
1850  const T a, const T b) const {
1851  if (transA) {
1852  if (transB) {
1853  C.resize(_n,B.m());
1854  if (b) {
1855  C.scal(b);
1856  } else {
1857  C.setZeros();
1858  }
1859  Vector<T> rowC(B.m());
1860  Vector<T> colA;
1861  for (int i = 0; i<_n; ++i) {
1862  this->refCol(i,colA);
1863  B.mult(colA,rowC,a);
1864  C.addRow(i,rowC,a);
1865  }
1866  } else {
1867  C.resize(_n,B.n());
1868  if (b) {
1869  C.scal(b);
1870  } else {
1871  C.setZeros();
1872  }
1873  Vector<T> colC;
1874  SpVector<T> colB;
1875  for (int i = 0; i<B.n(); ++i) {
1876  C.refCol(i,colC);
1877  B.refCol(i,colB);
1878  this->multTrans(colB,colC,a,T(1.0));
1879  }
1880  }
1881  } else {
1882  if (transB) {
1883  C.resize(_m,B.m());
1884  if (b) {
1885  C.scal(b);
1886  } else {
1887  C.setZeros();
1888  }
1889  Vector<T> colA;
1890  SpVector<T> colB;
1891  for (int i = 0; i<_n; ++i) {
1892  this->refCol(i,colA);
1893  B.refCol(i,colB);
1894  C.rank1Update(colA,colB,a);
1895  }
1896  } else {
1897  C.resize(_m,B.n());
1898  if (b) {
1899  C.scal(b);
1900  } else {
1901  C.setZeros();
1902  }
1903  Vector<T> colC;
1904  SpVector<T> colB;
1905  for (int i = 0; i<B.n(); ++i) {
1906  C.refCol(i,colC);
1907  B.refCol(i,colB);
1908  this->mult(colB,colC,a,T(1.0));
1909  }
1910  }
1911  };
1912 }
1913 
1914 
1916 template <typename T>
1917  inline void Matrix<T>::multDiagLeft(const Vector<T>& diag) {
1918  if (diag.n() != _m)
1919  return;
1920  T* d = diag.rawX();
1921  for (int i = 0; i< _n; ++i) {
1922  for (int j = 0; j<_m; ++j) {
1923  _X[i*_m+j] *= d[j];
1924  }
1925  }
1926  };
1927 
1929 template <typename T> inline void Matrix<T>::multDiagRight(
1930  const Vector<T>& diag) {
1931  if (diag.n() != _n)
1932  return;
1933  T* d = diag.rawX();
1934  for (int i = 0; i< _n; ++i) {
1935  for (int j = 0; j<_m; ++j) {
1936  _X[i*_m+j] *= d[i];
1937  }
1938  }
1939 };
1940 
1942 template <typename T> inline void Matrix<T>::mult_elementWise(
1943  const Matrix<T>& B, Matrix<T>& C) const {
1944  assert(_n == B._n && _m == B._m);
1945  C.resize(_m,_n);
1946  vMul<T>(_n*_m,_X,B._X,C._X);
1947 };
1948 
1950 template <typename T> inline void Matrix<T>::div_elementWise(
1951  const Matrix<T>& B, Matrix<T>& C) const {
1952  assert(_n == B._n && _m == B._m);
1953  C.resize(_m,_n);
1954  vDiv<T>(_n*_m,_X,B._X,C._X);
1955 };
1956 
1957 
1959 template <typename T> inline void Matrix<T>::XtX(Matrix<T>& xtx) const {
1960  xtx.resize(_n,_n);
1961  cblas_syrk<T>(CblasColMajor,CblasUpper,CblasTrans,_n,_m,T(1.0),
1962  _X,_m,T(),xtx._X,_n);
1963  xtx.fillSymmetric();
1964 };
1965 
1967 template <typename T> inline void Matrix<T>::XXt(Matrix<T>& xxt) const {
1968  xxt.resize(_m,_m);
1969  cblas_syrk<T>(CblasColMajor,CblasUpper,CblasNoTrans,_m,_n,T(1.0),
1970  _X,_m,T(),xxt._X,_m);
1971  xxt.fillSymmetric();
1972 };
1973 
1975 template <typename T> inline void Matrix<T>::upperTriXXt(Matrix<T>& XXt, const int L) const {
1976  XXt.resize(L,L);
1977  for (int i = 0; i<L; ++i) {
1978  cblas_syr<T>(CblasColMajor,CblasUpper,i+1,T(1.0),_X+i*_m,1,XXt._X,L);
1979  }
1980  XXt.fillSymmetric();
1981 }
1982 
1983 
1985 template <typename T> inline void Matrix<T>::diag(Vector<T>& dv) const {
1986  int size_diag=MIN(_n,_m);
1987  dv.resize(size_diag);
1988  T* const d = dv.rawX();
1989  for (int i = 0; i<size_diag; ++i)
1990  d[i]=_X[i*_m+i];
1991 };
1992 
1994 template <typename T> inline void Matrix<T>::setDiag(const Vector<T>& dv) {
1995  int size_diag=MIN(_n,_m);
1996  T* const d = dv.rawX();
1997  for (int i = 0; i<size_diag; ++i)
1998  _X[i*_m+i]=d[i];
1999 };
2000 
2002 template <typename T> inline void Matrix<T>::setDiag(const T val) {
2003  int size_diag=MIN(_n,_m);
2004  for (int i = 0; i<size_diag; ++i)
2005  _X[i*_m+i]=val;
2006 };
2007 
2008 
2010 template <typename T> inline void Matrix<T>::exp() {
2011  vExp<T>(_n*_m,_X,_X);
2012 };
2013 template <typename T> inline void Matrix<T>::Sqrt() {
2014  vSqrt<T>(_n*_m,_X,_X);
2015 };
2016 
2017 template <typename T> inline void Matrix<T>::Invsqrt() {
2018  vInvSqrt<T>(_n*_m,_X,_X);
2019 };
2021 template <typename T> inline T Matrix<T>::quad(
2022  const SpVector<T>& vec) const {
2023  T sum = T();
2024  int L = vec._L;
2025  int* r = vec._r;
2026  T* v = vec._v;
2027  for (int i = 0; i<L; ++i)
2028  for (int j = 0; j<L; ++j)
2029  sum += _X[r[i]*_m+r[j]]*v[i]*v[j];
2030  return sum;
2031 };
2032 
2033 template <typename T> inline void Matrix<T>::quad_mult(const Vector<T>& vec1,
2034  const SpVector<T>& vec2, Vector<T>& y, const T a, const T b) const {
2035  const int size_y= y.n();
2036  const int nn = _n/size_y;
2037  //y.resize(size_y);
2038  //y.setZeros();
2039  Matrix<T> tmp;
2040  for (int i = 0; i<size_y; ++i) {
2041  tmp.setData(_X+(i*nn)*_m,_m,nn);
2042  y[i]=b*y[i]+a*tmp.quad(vec1,vec2);
2043  }
2044 }
2045 
2047 template <typename T> inline T Matrix<T>::quad(
2048  const Vector<T>& vec1, const SpVector<T>& vec) const {
2049  T sum = T();
2050  int L = vec._L;
2051  int* r = vec._r;
2052  T* v = vec._v;
2053  Vector<T> col;
2054  for (int i = 0; i<L; ++i) {
2055  this->refCol(r[i],col);
2056  sum += v[i]*col.dot(vec1);
2057  }
2058  return sum;
2059 };
2060 
2062 template <typename T> inline void Matrix<T>::add(const Matrix<T>& mat, const T alpha) {
2063  assert(mat._m == _m && mat._n == _n);
2064  cblas_axpy<T>(_n*_m,alpha,mat._X,1,_X,1);
2065 };
2066 
2068 template <typename T> inline T Matrix<T>::dot(const Matrix<T>& mat) const {
2069  assert(mat._m == _m && mat._n == _n);
2070  return cblas_dot<T>(_n*_m,mat._X,1,_X,1);
2071 };
2072 
2073 
2075 template <typename T> inline void Matrix<T>::add(const T alpha) {
2076  for (int i = 0; i<_n*_m; ++i) _X[i]+=alpha;
2077 };
2078 
2080 template <typename T> inline void Matrix<T>::sub(const Matrix<T>& mat) {
2081  vSub<T>(_n*_m,_X,mat._X,_X);
2082 };
2083 
2085 template <typename T> inline T Matrix<T>::asum() const {
2086  return cblas_asum<T>(_n*_m,_X,1);
2087 };
2088 
2090 template <typename T> inline T Matrix<T>::trace() const {
2091  T sum=T();
2092  int m = MIN(_n,_m);
2093  for (int i = 0; i<m; ++i)
2094  sum += _X[i*_m+i];
2095  return sum;
2096 };
2097 
2099 template <typename T> inline T Matrix<T>::normF() const {
2100  return cblas_nrm2<T>(_n*_m,_X,1);
2101 };
2102 
2103 template <typename T> inline T Matrix<T>::mean() const {
2104  Vector<T> vec;
2105  this->toVect(vec);
2106  return vec.mean();
2107 };
2108 
2110 template <typename T> inline T Matrix<T>::normFsq() const {
2111  return cblas_dot<T>(_n*_m,_X,1,_X,1);
2112 };
2113 
2115 template <typename T> inline T Matrix<T>::norm_inf_2_col() const {
2116  Vector<T> col;
2117  T max = -1.0;
2118  for (int i = 0; i<_n; ++i) {
2119  refCol(i,col);
2120  T norm_col = col.nrm2();
2121  if (norm_col > max)
2122  max = norm_col;
2123  }
2124  return max;
2125 };
2126 
2128 template <typename T> inline T Matrix<T>::norm_1_2_col() const {
2129  Vector<T> col;
2130  T sum = 0.0;
2131  for (int i = 0; i<_n; ++i) {
2132  refCol(i,col);
2133  sum += col.nrm2();
2134  }
2135  return sum;
2136 };
2137 
2139 template <typename T> inline void Matrix<T>::norm_2_rows(
2140  Vector<T>& norms) const {
2141  norms.resize(_m);
2142  norms.setZeros();
2143  for (int i = 0; i<_n; ++i)
2144  for (int j = 0; j<_m; ++j)
2145  norms[j] += _X[i*_m+j]*_X[i*_m+j];
2146  for (int j = 0; j<_m; ++j)
2147  norms[j]=sqrt(norms[j]);
2148 };
2149 
2151 template <typename T> inline void Matrix<T>::norm_2sq_rows(
2152  Vector<T>& norms) const {
2153  norms.resize(_m);
2154  norms.setZeros();
2155  for (int i = 0; i<_n; ++i)
2156  for (int j = 0; j<_m; ++j)
2157  norms[j] += _X[i*_m+j]*_X[i*_m+j];
2158 };
2159 
2160 
2162 template <typename T> inline void Matrix<T>::norm_2_cols(
2163  Vector<T>& norms) const {
2164  norms.resize(_n);
2165  Vector<T> col;
2166  for (int i = 0; i<_n; ++i) {
2167  refCol(i,col);
2168  norms[i] = col.nrm2();
2169  }
2170 };
2171 
2172 
2174 template <typename T> inline void Matrix<T>::norm_inf_cols(Vector<T>& norms) const {
2175  norms.resize(_n);
2176  Vector<T> col;
2177  for (int i = 0; i<_n; ++i) {
2178  refCol(i,col);
2179  norms[i] = col.fmaxval();
2180  }
2181 };
2182 
2184 template <typename T> inline void Matrix<T>::norm_inf_rows(Vector<T>& norms) const {
2185  norms.resize(_m);
2186  norms.setZeros();
2187  for (int i = 0; i<_n; ++i)
2188  for (int j = 0; j<_m; ++j)
2189  norms[j] = MAX(abs<T>(_X[i*_m+j]),norms[j]);
2190 };
2191 
2193 template <typename T> inline void Matrix<T>::norm_l1_rows(Vector<T>& norms) const {
2194  norms.resize(_m);
2195  norms.setZeros();
2196  for (int i = 0; i<_n; ++i)
2197  for (int j = 0; j<_m; ++j)
2198  norms[j] += abs<T>(_X[i*_m+j]);
2199 };
2200 
2201 
2202 
2204 template <typename T> inline void Matrix<T>::norm_2sq_cols(
2205  Vector<T>& norms) const {
2206  norms.resize(_n);
2207  Vector<T> col;
2208  for (int i = 0; i<_n; ++i) {
2209  refCol(i,col);
2210  norms[i] = col.nrm2sq();
2211  }
2212 };
2213 
2214 template <typename T>
2215 inline void Matrix<T>::sum_cols(Vector<T>& sum) const {
2216  sum.resize(_m);
2217  sum.setZeros();
2218  Vector<T> tmp;
2219  for (int i = 0; i<_n; ++i) {
2220  this->refCol(i,tmp);
2221  sum.add(tmp);
2222  }
2223 };
2224 
2226 template <typename T> inline void Matrix<T>::meanCol(Vector<T>& mean) const {
2227  Vector<T> ones(_n);
2228  ones.set(T(1.0/_n));
2229  this->mult(ones,mean,1.0,0.0);
2230 };
2231 
2233 template <typename T> inline void Matrix<T>::meanRow(Vector<T>& mean) const {
2234  Vector<T> ones(_m);
2235  ones.set(T(1.0/_m));
2236  this->multTrans(ones,mean,1.0,0.0);
2237 };
2238 
2239 
2241 template <typename T> inline void Matrix<T>::fillRow(const Vector<T>& row) {
2242  for (int i = 0; i<_n; ++i) {
2243  T val = row[i];
2244  for (int j = 0; j<_m; ++j) {
2245  _X[i*_m+j]=val;
2246  }
2247  }
2248 };
2249 
2251 template <typename T> inline void Matrix<T>::extractRow(const int j,
2252  Vector<T>& row) const {
2253  row.resize(_n);
2254  for (int i = 0; i<_n; ++i) {
2255  row[i]=_X[i*_m+j];
2256  }
2257 };
2258 
2260 template <typename T> inline void Matrix<T>::setRow(const int j,
2261  const Vector<T>& row) {
2262  for (int i = 0; i<_n; ++i) {
2263  _X[i*_m+j]=row[i];
2264  }
2265 };
2266 
2268 template <typename T> inline void Matrix<T>::addRow(const int j,
2269  const Vector<T>& row, const T a) {
2270  if (a==1.0) {
2271  for (int i = 0; i<_n; ++i) {
2272  _X[i*_m+j]+=row[i];
2273  }
2274  } else {
2275  for (int i = 0; i<_n; ++i) {
2276  _X[i*_m+j]+=a*row[i];
2277  }
2278  }
2279 };
2280 
2281 
2283 template <typename T> inline void Matrix<T>::softThrshold(const T nu) {
2284  Vector<T> vec;
2285  toVect(vec);
2286  vec.softThrshold(nu);
2287 };
2288 
2290 template <typename T> inline void Matrix<T>::hardThrshold(const T nu) {
2291  Vector<T> vec;
2292  toVect(vec);
2293  vec.hardThrshold(nu);
2294 };
2295 
2296 
2298 template <typename T> inline void Matrix<T>::thrsmax(const T nu) {
2299  Vector<T> vec;
2300  toVect(vec);
2301  vec.thrsmax(nu);
2302 };
2303 
2305 template <typename T> inline void Matrix<T>::thrsmin(const T nu) {
2306  Vector<T> vec;
2307  toVect(vec);
2308  vec.thrsmin(nu);
2309 };
2310 
2311 
2313 template <typename T> inline void Matrix<T>::inv_elem() {
2314  Vector<T> vec;
2315  toVect(vec);
2316  vec.inv();
2317 };
2318 
2320 template <typename T> inline void Matrix<T>::blockThrshold(const T nu,
2321  const int sizeGroup) {
2322  for (int i = 0; i<_n; ++i) {
2323  int j;
2324  for (j = 0; j<_m-sizeGroup+1; j+=sizeGroup) {
2325  T nrm=0;
2326  for (int k = 0; k<sizeGroup; ++k)
2327  nrm += _X[i*_m +j+k]*_X[i*_m +j+k];
2328  nrm=sqrt(nrm);
2329  if (nrm < nu) {
2330  for (int k = 0; k<sizeGroup; ++k)
2331  _X[i*_m +j+k]=0;
2332  } else {
2333  T scal = (nrm-nu)/nrm;
2334  for (int k = 0; k<sizeGroup; ++k)
2335  _X[i*_m +j+k]*=scal;
2336  }
2337  }
2338  j -= sizeGroup;
2339  for ( ; j<_m; ++j)
2340  _X[j]=softThrs<T>(_X[j],nu);
2341  }
2342 }
2343 
2344 template <typename T> inline void Matrix<T>::sparseProject(Matrix<T>& Y,
2345  const T thrs, const int mode, const T lambda1,
2346  const T lambda2, const T lambda3, const bool pos,
2347  const int numThreads) {
2348 
2349  int NUM_THREADS=init_omp(numThreads);
2350  Vector<T>* XXT= new Vector<T>[NUM_THREADS];
2351  for (int i = 0; i<NUM_THREADS; ++i) {
2352  XXT[i].resize(_m);
2353  }
2354 
2355  int i;
2356 #pragma omp parallel for private(i)
2357  for (i = 0; i< _n; ++i) {
2358 #ifdef _OPENMP
2359  int numT=omp_get_thread_num();
2360 #else
2361  int numT=0;
2362 #endif
2363  Vector<T> Xi;
2364  this->refCol(i,Xi);
2365  Vector<T> Yi;
2366  Y.refCol(i,Yi);
2367  Vector<T>& XX = XXT[numT];
2368  XX.copy(Xi);
2369  XX.sparseProject(Yi,thrs,mode,lambda1,lambda2,lambda3,pos);
2370  }
2371  delete[](XXT);
2372 };
2373 
2374 
2376 template <typename T> inline void Matrix<T>::thrsPos() {
2377  Vector<T> vec;
2378  toVect(vec);
2379  vec.thrsPos();
2380 };
2381 
2382 
2384 template <typename T> inline void Matrix<T>::rank1Update(
2385  const Vector<T>& vec1, const Vector<T>& vec2, const T alpha) {
2386  cblas_ger<T>(CblasColMajor,_m,_n,alpha,vec1._X,1,vec2._X,1,_X,_m);
2387 };
2388 
2390 template <typename T> inline void Matrix<T>::rank1Update(
2391  const SpVector<T>& vec1, const Vector<T>& vec2, const T alpha) {
2392  int* r = vec1._r;
2393  T* v = vec1._v;
2394  T* X2 = vec2._X;
2395  assert(vec2._n == _n);
2396  if (alpha == 1.0) {
2397  for (int i = 0; i<_n; ++i) {
2398  for (int j = 0; j<vec1._L; ++j) {
2399  _X[i*_m+r[j]] += v[j]*X2[i];
2400  }
2401  }
2402  } else {
2403  for (int i = 0; i<_n; ++i) {
2404  for (int j = 0; j<vec1._L; ++j) {
2405  _X[i*_m+r[j]] += alpha*v[j]*X2[i];
2406  }
2407  }
2408  }
2409 };
2410 
2411 template <typename T>
2412 inline void Matrix<T>::rank1Update_mult(const Vector<T>& vec1,
2413  const Vector<T>& vec1b,
2414  const SpVector<T>& vec2,
2415  const T alpha) {
2416  const int nn = vec1b.n();
2417  const int size_A = _n/nn;
2418  Matrix<T> tmp;
2419  for (int i = 0; i<nn; ++i) {
2420  tmp.setData(_X+i*size_A*_m,_m,size_A);
2421  tmp.rank1Update(vec1,vec2,alpha*vec1b[i]);
2422  }
2423 };
2424 
2426 template <typename T> inline void Matrix<T>::rank1Update(
2427  const SpVector<T>& vec1, const SpVector<T>& vec2, const T alpha) {
2428  int* r = vec1._r;
2429  T* v = vec1._v;
2430  T* v2 = vec2._v;
2431  int* r2 = vec2._r;
2432  if (alpha == 1.0) {
2433  for (int i = 0; i<vec2._L; ++i) {
2434  for (int j = 0; j<vec1._L; ++j) {
2435  _X[r2[i]*_m+r[j]] += v[j]*v2[i];
2436  }
2437  }
2438  } else {
2439  for (int i = 0; i<vec2._L; ++i) {
2440  for (int j = 0; j<vec1._L; ++j) {
2441  _X[r[i]*_m+r[j]] += alpha*v[j]*v2[i];
2442  }
2443  }
2444  }
2445 };
2446 
2447 
2449 template <typename T> inline void Matrix<T>::rank1Update(
2450  const Vector<T>& vec1, const SpVector<T>& vec2, const T alpha) {
2451  int* r = vec2._r;
2452  T* v = vec2._v;
2453  Vector<T> Xi;
2454  for (int i = 0; i<vec2._L; ++i) {
2455  this->refCol(r[i],Xi);
2456  Xi.add(vec1,v[i]*alpha);
2457  }
2458 };
2459 
2461 template <typename T> inline void Matrix<T>::rank1Update(
2462  const SpVector<T>& vec1, const T alpha) {
2463  int* r = vec1._r;
2464  T* v = vec1._v;
2465  if (alpha == 1.0) {
2466  for (int i = 0; i<vec1._L; ++i) {
2467  for (int j = 0; j<vec1._L; ++j) {
2468  _X[r[i]*_m+r[j]] += v[j]*v[i];
2469  }
2470  }
2471  } else {
2472  for (int i = 0; i<vec1._L; ++i) {
2473  for (int j = 0; j<vec1._L; ++j) {
2474  _X[_m*r[i]+r[j]] += alpha*v[j]*v[i];
2475  }
2476  }
2477  }
2478 };
2479 
2480 
2482 template <typename T> inline void Matrix<T>::conjugateGradient(
2483  const Vector<T>& b, Vector<T>& x, const T tol, const int itermax) const {
2484  Vector<T> R,P,AP;
2485  R.copy(b);
2486  this->mult(x,R,T(-1.0),T(1.0));
2487  P.copy(R);
2488  int k = 0;
2489  T normR = R.nrm2sq();
2490  T alpha;
2491  while (normR > tol && k < itermax) {
2492  this->mult(P,AP);
2493  alpha = normR/P.dot(AP);
2494  x.add(P,alpha);
2495  R.add(AP,-alpha);
2496  T tmp = R.nrm2sq();
2497  P.scal(tmp/normR);
2498  normR = tmp;
2499  P.add(R,T(1.0));
2500  ++k;
2501  };
2502 };
2503 
2504 template <typename T> inline void Matrix<T>::drop(char* fileName) const {
2505  std::ofstream f;
2506  f.precision(12);
2507  f.flags(std::ios_base::scientific);
2508  f.open(fileName, ofstream::trunc);
2509  std::cout << "Matrix written in " << fileName << std::endl;
2510  for (int i = 0; i<_n; ++i) {
2511  for (int j = 0; j<_m; ++j)
2512  f << _X[i*_m+j] << " ";
2513  f << std::endl;
2514  }
2515  f.close();
2516 };
2517 
2519 template <typename T> inline void Matrix<T>::NadarayaWatson(
2520  const Vector<int>& ind, const T sigma) {
2521  if (ind.n() != _n) return;
2522 
2523  init_omp(MAX_THREADS);
2524 
2525  const int Ngroups=ind.maxval();
2526  int i;
2527 #pragma omp parallel for private(i)
2528  for (i = 1; i<=Ngroups; ++i) {
2529  Vector<int> indicesGroup(_n);
2530  int count = 0;
2531  for (int j = 0; j<_n; ++j)
2532  if (ind[j] == i) indicesGroup[count++]=j;
2533  Matrix<T> Xm(_m,count);
2534  Vector<T> col, col2;
2535  for (int j= 0; j<count; ++j) {
2536  this->refCol(indicesGroup[j],col);
2537  Xm.refCol(j,col2);
2538  col2.copy(col);
2539  }
2540  Vector<T> norms;
2541  Xm.norm_2sq_cols(norms);
2542  Matrix<T> weights;
2543  Xm.XtX(weights);
2544  weights.scal(T(-2.0));
2545  Vector<T> ones(Xm.n());
2546  ones.set(T(1.0));
2547  weights.rank1Update(ones,norms);
2548  weights.rank1Update(norms,ones);
2549  weights.scal(-sigma);
2550  weights.exp();
2551  Vector<T> den;
2552  weights.mult(ones,den);
2553  den.inv();
2554  weights.multDiagRight(den);
2555  Matrix<T> num;
2556  Xm.mult(weights,num);
2557  for (int j= 0; j<count; ++j) {
2558  this->refCol(indicesGroup[j],col);
2559  num.refCol(j,col2);
2560  col.copy(col2);
2561  }
2562  }
2563 };
2564 
2566 template <typename T> inline void Matrix<T>::toSparse(SpMatrix<T>& out) const {
2567  out.clear();
2568  int count=0;
2569  int* pB;
2570 #pragma omp critical
2571  {
2572  pB=new int[_n+1];
2573  }
2574  int* pE=pB+1;
2575  for (int i = 0; i<_n*_m; ++i)
2576  if (_X[i] != 0) ++count;
2577  int* r;
2578  T* v;
2579 #pragma omp critical
2580  {
2581  r=new int[count];
2582  v=new T[count];
2583  }
2584  count=0;
2585  for (int i = 0; i<_n; ++i) {
2586  pB[i]=count;
2587  for (int j = 0; j<_m; ++j) {
2588  if (_X[i*_m+j] != 0) {
2589  v[count]=_X[i*_m+j];
2590  r[count++]=j;
2591  }
2592  }
2593  pE[i]=count;
2594  }
2595  out._v=v;
2596  out._r=r;
2597  out._pB=pB;
2598  out._pE=pE;
2599  out._m=_m;
2600  out._n=_n;
2601  out._nzmax=count;
2602  out._externAlloc=false;
2603 };
2604 
2606 template <typename T> inline void Matrix<T>::toSparseTrans(
2607  SpMatrix<T>& out) {
2608  out.clear();
2609  int count=0;
2610  int* pB;
2611 #pragma omp critical
2612  {
2613  pB=new int[_m+1];
2614  }
2615  int* pE=pB+1;
2616  for (int i = 0; i<_n*_m; ++i)
2617  if (_X[i] != 0) ++count;
2618  int* r;
2619  T* v;
2620 #pragma omp critical
2621  {
2622  r=new int[count];
2623  v=new T[count];
2624  }
2625  count=0;
2626  for (int i = 0; i<_m; ++i) {
2627  pB[i]=count;
2628  for (int j = 0; j<_n; ++j) {
2629  if (_X[i+j*_m] != 0) {
2630  v[count]=_X[j*_m+i];
2631  r[count++]=j;
2632  }
2633  }
2634  pE[i]=count;
2635  }
2636  out._v=v;
2637  out._r=r;
2638  out._pB=pB;
2639  out._pE=pE;
2640  out._m=_n;
2641  out._n=_m;
2642  out._nzmax=count;
2643  out._externAlloc=false;
2644 };
2645 
2647 template <typename T> inline void Matrix<T>::toVect(
2648  Vector<T>& vec) const {
2649  vec.clear();
2650  vec._externAlloc=true;
2651  vec._n=_n*_m;
2652  vec._X=_X;
2653 };
2654 
2656 template <typename T> inline void Matrix<T>::merge(const Matrix<T>& B,
2657  Matrix<T>& C) const {
2658  const int K =_n;
2659  Matrix<T> G;
2660  this->mult(B,G,true,false);
2661  std::list<int> list;
2662  for (int i = 0; i<G.n(); ++i) {
2663  Vector<T> g;
2664  G.refCol(i,g);
2665  T fmax=g.fmaxval();
2666  if (fmax < 0.995) list.push_back(i);
2667  }
2668  C.resize(_m,K+list.size());
2669 
2670  for (int i = 0; i<K; ++i) {
2671  Vector<T> d, d2;
2672  C.refCol(i,d);
2673  this->refCol(i,d2);
2674  d.copy(d2);
2675  }
2676  int count=0;
2677  for (std::list<int>::const_iterator it = list.begin();
2678  it != list.end(); ++it) {
2679  Vector<T> d, d2;
2680  C.refCol(K+count,d);
2681  B.refCol(*it,d2);
2682  d.copy(d2);
2683  ++count;
2684  }
2685 };
2686 
2687 /* ***********************************
2688  * Implementation of the class Vector
2689  * ***********************************/
2690 
2691 
2693 template <typename T> Vector<T>::Vector() :
2694  _externAlloc(true), _X(NULL), _n(0) { };
2695 
2697 template <typename T> Vector<T>::Vector(int n) :
2698  _externAlloc(false), _n(n) {
2699 #pragma omp critical
2700  {
2701  _X=new T[_n];
2702  }
2703  };
2704 
2706 template <typename T> Vector<T>::Vector(T* X, int n) :
2707  _externAlloc(true), _X(X), _n(n) { };
2708 
2710 template <typename T> Vector<T>::Vector(const Vector<T>& vec) :
2711  _externAlloc(false), _n(vec._n) {
2712 #pragma omp critical
2713  {
2714  _X=new T[_n];
2715  }
2716  cblas_copy<T>(_n,vec._X,1,_X,1);
2717  };
2718 
2720 template <typename T> Vector<T>::~Vector() {
2721  clear();
2722 };
2723 
2725 template <> inline void Vector<double>::print(const char* name) const {
2726  printf("%s, %d\n",name,_n);
2727  for (int i = 0; i<_n; ++i) {
2728  printf("%g ",_X[i]);
2729  }
2730  printf("\n");
2731 };
2732 
2734 template <> inline void Vector<float>::print(const char* name) const {
2735  printf("%s, %d\n",name,_n);
2736  for (int i = 0; i<_n; ++i) {
2737  printf("%g ",_X[i]);
2738  }
2739  printf("\n");
2740 };
2741 
2743 template <> inline void Vector<int>::print(const char* name) const {
2744  printf("%s, %d\n",name,_n);
2745  for (int i = 0; i<_n; ++i) {
2746  printf("%d ",_X[i]);
2747  }
2748  printf("\n");
2749 };
2750 
2752 template <> inline void Vector<bool>::print(const char* name) const {
2753  printf("%s, %d\n",name,_n);
2754  for (int i = 0; i<_n; ++i) {
2755  printf("%d ",_X[i] ? 1 : 0);
2756  }
2757  printf("\n");
2758 };
2759 
2761 template <typename T> inline int Vector<T>::max() const {
2762  int imax=0;
2763  T max=_X[0];
2764  for (int j = 1; j<_n; ++j) {
2765  T cur = _X[j];
2766  if (cur > max) {
2767  imax=j;
2768  max = cur;
2769  }
2770  }
2771  return imax;
2772 };
2773 
2775 template <typename T> inline int Vector<T>::min() const {
2776  int imin=0;
2777  T min=_X[0];
2778  for (int j = 1; j<_n; ++j) {
2779  T cur = _X[j];
2780  if (cur < min) {
2781  imin=j;
2782  min = cur;
2783  }
2784  }
2785  return imin;
2786 };
2787 
2789 template <typename T> inline T Vector<T>::maxval() const {
2790  return _X[this->max()];
2791 };
2792 
2794 template <typename T> inline T Vector<T>::minval() const {
2795  return _X[this->min()];
2796 };
2797 
2799 template <typename T> inline T Vector<T>::fmaxval() const {
2800  return fabs(_X[this->fmax()]);
2801 };
2802 
2804 template <typename T> inline T Vector<T>::fminval() const {
2805  return fabs(_X[this->fmin()]);
2806 };
2807 
2808 template <typename T>
2809 inline void Vector<T>::logspace(const int n, const T a, const T b) {
2810  T first=log10(a);
2811  T last=log10(b);
2812  T step = (last-first)/(n-1);
2813  this->resize(n);
2814  _X[0]=first;
2815  for (int i = 1; i<_n; ++i)
2816  _X[i]=_X[i-1]+step;
2817  for (int i = 0; i<_n; ++i)
2818  _X[i]=pow(T(10.0),_X[i]);
2819 }
2820 
2821 template <typename T>
2822 inline int Vector<T>::nnz() const {
2823  int sum=0;
2824  for (int i = 0; i<_n; ++i)
2825  if (_X[i] != T()) ++sum;
2826  return sum;
2827 };
2829 template <>
2830 inline void Vector<int>::logspace(const int n, const int a, const int b) {
2831  Vector<double> tmp(n);
2832  tmp.logspace(n,double(a),double(b));
2833  this->resize(n);
2834  _X[0]=a;
2835  _X[n-1]=b;
2836  for (int i = 1; i<_n-1; ++i) {
2837  int candidate=static_cast<int>(floor(static_cast<double>(tmp[i])));
2838  _X[i]= candidate > _X[i-1] ? candidate : _X[i-1]+1;
2839  }
2840 }
2841 
2843 template <typename T> inline int Vector<T>::fmax() const {
2844  return cblas_iamax<T>(_n,_X,1);
2845 };
2846 
2848 template <typename T> inline int Vector<T>::fmin() const {
2849  return cblas_iamin<T>(_n,_X,1);
2850 };
2851 
2853 template <typename T> inline T& Vector<T>::operator[] (const int i) {
2854  assert(i>=0 && i<_n);
2855  return _X[i];
2856 };
2857 
2859 template <typename T> inline T Vector<T>::operator[] (const int i) const {
2860  assert(i>=0 && i<_n);
2861  return _X[i];
2862 };
2863 
2865 template <typename T> inline void Vector<T>::copy(const Vector<T>& x) {
2866  this->resize(x.n());
2867  cblas_copy<T>(_n,x._X,1,_X,1);
2868 };
2869 
2871 template <typename T> inline void Vector<T>::setZeros() {
2872  memset(_X,0,_n*sizeof(T));
2873 };
2874 
2876 template <typename T> inline void Vector<T>::resize(const int n) {
2877  if (_n == n) return;
2878  clear();
2879 #pragma omp critical
2880  {
2881  _X=new T[n];
2882  }
2883  _n=n;
2884  _externAlloc=false;
2885  this->setZeros();
2886 };
2887 
2889 template <typename T> inline void Vector<T>::setPointer(T* X, const int n) {
2890  clear();
2891  _externAlloc=true;
2892  _X=X;
2893  _n=n;
2894 };
2895 
2897 template <> inline void Vector<int>::randperm(int n) {
2898  resize(n);
2899  Vector<int> table(n);
2900  for (int i = 0; i<n; ++i)
2901  table[i]=i;
2902  int size=n;
2903  for (int i = 0; i<n; ++i) {
2904  const int ind=random() % size;
2905  _X[i]=table[ind];
2906  table[ind]=table[size-1];
2907  --size;
2908  }
2909 };
2910 
2912 template <typename T> inline void Vector<T>::setAleat() {
2913  for (int i = 0; i<_n; ++i) _X[i]=normalDistrib<T>();
2914 };
2915 
2917 template <typename T> inline void Vector<T>::clear() {
2918  if (!_externAlloc) delete[](_X);
2919  _n=0;
2920  _X=NULL;
2921  _externAlloc=true;
2922 };
2923 
2925 template <typename T> inline void Vector<T>::softThrshold(const T nu) {
2926  for (int i = 0; i<_n; ++i) {
2927  if (_X[i] > nu) {
2928  _X[i] -= nu;
2929  } else if (_X[i] < -nu) {
2930  _X[i] += nu;
2931  } else {
2932  _X[i] = T();
2933  }
2934  }
2935 };
2936 
2938 template <typename T> inline void Vector<T>::hardThrshold(const T nu) {
2939  for (int i = 0; i<_n; ++i) {
2940  if (!(_X[i] > nu || _X[i] < -nu)) {
2941  _X[i] = 0;
2942  }
2943  }
2944 };
2945 
2946 
2948 template <typename T> inline void Vector<T>::thrsmax(const T nu) {
2949  for (int i = 0; i<_n; ++i)
2950  _X[i]=MAX(_X[i],nu);
2951 }
2952 
2954 template <typename T> inline void Vector<T>::thrsmin(const T nu) {
2955  for (int i = 0; i<_n; ++i)
2956  _X[i]=MIN(_X[i],nu);
2957 }
2958 
2960 template <typename T> inline void Vector<T>::thrsabsmin(const T nu) {
2961  for (int i = 0; i<_n; ++i)
2962  _X[i]=MAX(MIN(_X[i],nu),-nu);
2963 }
2964 
2965 
2967 template <typename T> inline void Vector<T>::thrshold(const T nu) {
2968  for (int i = 0; i<_n; ++i)
2969  if (abs<T>(_X[i]) < nu)
2970  _X[i]=0;
2971 }
2973 template <typename T> inline void Vector<T>::thrsPos() {
2974  for (int i = 0; i<_n; ++i) {
2975  if (_X[i] < 0) _X[i]=0;
2976  }
2977 };
2978 
2979 template <>
2980 inline bool Vector<bool>::alltrue() const {
2981  for (int i = 0; i<_n; ++i) {
2982  if (!_X[i]) return false;
2983  }
2984  return true;
2985 };
2986 
2987 template <>
2988 inline bool Vector<bool>::allfalse() const {
2989  for (int i = 0; i<_n; ++i) {
2990  if (_X[i]) return false;
2991  }
2992  return true;
2993 };
2994 
2995 
2997 template <typename T> inline void Vector<T>::set(const T val) {
2998  for (int i = 0; i<_n; ++i) _X[i]=val;
2999 };
3000 
3002 template <typename T> inline T Vector<T>::nrm2() const {
3003  return cblas_nrm2<T>(_n,_X,1);
3004 };
3005 
3007 template <typename T> inline T Vector<T>::nrm2sq() const {
3008  return cblas_dot<T>(_n,_X,1,_X,1);
3009 };
3010 
3012 template <typename T> inline T Vector<T>::dot(const Vector<T>& x) const {
3013  assert(_n == x._n);
3014  return cblas_dot<T>(_n,_X,1,x._X,1);
3015 };
3016 
3018 template <typename T> inline T Vector<T>::dot(const SpVector<T>& x) const {
3019  T sum=0;
3020  const T* v = x._v;
3021  const int* r = x._r;
3022  for (int i = 0; i<x._L; ++i) {
3023  sum += _X[r[i]]*v[i];
3024  }
3025  return sum;
3026 };
3027 
3029 template <typename T> inline void Vector<T>::add(const Vector<T>& x, const T a) {
3030  assert(_n == x._n);
3031  cblas_axpy<T>(_n,a,x._X,1,_X,1);
3032 };
3033 
3035 template <typename T> inline void Vector<T>::add(const SpVector<T>& x,
3036  const T a) {
3037  if (a == 1.0) {
3038  for (int i = 0; i<x._L; ++i)
3039  _X[x._r[i]]+=x._v[i];
3040  } else {
3041  for (int i = 0; i<x._L; ++i)
3042  _X[x._r[i]]+=a*x._v[i];
3043  }
3044 };
3045 
3047 template <typename T> inline void Vector<T>::add(const T a) {
3048  for (int i = 0; i<_n; ++i) _X[i]+=a;
3049 };
3050 
3052 template <typename T> inline void Vector<T>::sub(const Vector<T>& x) {
3053  assert(_n == x._n);
3054  vSub<T>(_n,_X,x._X,_X);
3055 };
3056 
3058 template <typename T> inline void Vector<T>::sub(const SpVector<T>& x) {
3059  for (int i = 0; i<x._L; ++i)
3060  _X[x._r[i]]-=x._v[i];
3061 };
3062 
3064 template <typename T> inline void Vector<T>::div(const Vector<T>& x) {
3065  assert(_n == x._n);
3066  vDiv<T>(_n,_X,x._X,_X);
3067 };
3068 
3070 template <typename T> inline void Vector<T>::div(const Vector<T>& x, const Vector<T>& y) {
3071  assert(_n == x._n);
3072  vDiv<T>(_n,x._X,y._X,_X);
3073 };
3074 
3075 
3077 template <typename T> inline void Vector<T>::sqr(const Vector<T>& x) {
3078  this->resize(x._n);
3079  vSqr<T>(_n,x._X,_X);
3080 }
3081 
3083 template <typename T> inline void Vector<T>::Invsqrt(const Vector<T>& x) {
3084  this->resize(x._n);
3085  vInvSqrt<T>(_n,x._X,_X);
3086 }
3088 template <typename T> inline void Vector<T>::Sqrt(const Vector<T>& x) {
3089  this->resize(x._n);
3090  vSqrt<T>(_n,x._X,_X);
3091 }
3093 template <typename T> inline void Vector<T>::Invsqrt() {
3094  vInvSqrt<T>(_n,_X,_X);
3095 }
3097 template <typename T> inline void Vector<T>::Sqrt() {
3098  vSqrt<T>(_n,_X,_X);
3099 }
3100 
3101 
3103 template <typename T> inline void Vector<T>::inv(const Vector<T>& x) {
3104  this->resize(x.n());
3105  vInv<T>(_n,x._X,_X);
3106 };
3107 
3109 template <typename T> inline void Vector<T>::inv() {
3110  vInv<T>(_n,_X,_X);
3111 };
3112 
3114 template <typename T> inline void Vector<T>::mult(const Vector<T>& x,
3115  const Vector<T>& y) {
3116  this->resize(x.n());
3117  vMul<T>(_n,x._X,y._X,_X);
3118 };
3119 ;
3120 
3122 template <typename T> inline void Vector<T>::normalize() {
3123  T norm=nrm2();
3124  if (norm > EPSILON) scal(1.0/norm);
3125 };
3126 
3128 template <typename T> inline void Vector<T>::normalize2() {
3129  T norm=nrm2();
3130  if (norm > T(1.0)) scal(1.0/norm);
3131 };
3132 
3134 template <typename T> inline void Vector<T>::whiten(
3135  Vector<T>& meanv, const bool pattern) {
3136  if (pattern) {
3137  const int n =static_cast<int>(sqrt(static_cast<T>(_n)));
3138  int count[4];
3139  for (int i = 0; i<4; ++i) count[i]=0;
3140  int offsetx=0;
3141  for (int j = 0; j<n; ++j) {
3142  offsetx= (offsetx+1) % 2;
3143  int offsety=0;
3144  for (int k = 0; k<n; ++k) {
3145  offsety= (offsety+1) % 2;
3146  meanv[2*offsetx+offsety]+=_X[j*n+k];
3147  count[2*offsetx+offsety]++;
3148  }
3149  }
3150  for (int i = 0; i<4; ++i)
3151  meanv[i] /= count[i];
3152  offsetx=0;
3153  for (int j = 0; j<n; ++j) {
3154  offsetx= (offsetx+1) % 2;
3155  int offsety=0;
3156  for (int k = 0; k<n; ++k) {
3157  offsety= (offsety+1) % 2;
3158  _X[j*n+k]-=meanv[2*offsetx+offsety];
3159  }
3160  }
3161  } else {
3162  const int V = meanv.n();
3163  const int sizePatch=_n/V;
3164  for (int j = 0; j<V; ++j) {
3165  T mean = 0;
3166  for (int k = 0; k<sizePatch; ++k) {
3167  mean+=_X[sizePatch*j+k];
3168  }
3169  mean /= sizePatch;
3170  for (int k = 0; k<sizePatch; ++k) {
3171  _X[sizePatch*j+k]-=mean;
3172  }
3173  meanv[j]=mean;
3174  }
3175  }
3176 };
3177 
3179 template <typename T> inline void Vector<T>::whiten(
3180  Vector<T>& meanv, const Vector<T>& mask) {
3181  const int V = meanv.n();
3182  const int sizePatch=_n/V;
3183  for (int j = 0; j<V; ++j) {
3184  T mean = 0;
3185  for (int k = 0; k<sizePatch; ++k) {
3186  mean+=_X[sizePatch*j+k];
3187  }
3188  mean /= cblas_asum(sizePatch,mask._X+j*sizePatch,1);
3189  for (int k = 0; k<sizePatch; ++k) {
3190  if (mask[sizePatch*j+k])
3191  _X[sizePatch*j+k]-=mean;
3192  }
3193  meanv[j]=mean;
3194  }
3195 };
3196 
3198 template <typename T> inline void Vector<T>::whiten(const int V) {
3199  const int sizePatch=_n/V;
3200  for (int j = 0; j<V; ++j) {
3201  T mean = 0;
3202  for (int k = 0; k<sizePatch; ++k) {
3203  mean+=_X[sizePatch*j+k];
3204  }
3205  mean /= sizePatch;
3206  for (int k = 0; k<sizePatch; ++k) {
3207  _X[sizePatch*j+k]-=mean;
3208  }
3209  }
3210 };
3211 
3212 template <typename T> inline T Vector<T>::KL(const Vector<T>& Y) {
3213  T sum = 0;
3214  T* prY = Y.rawX();
3215  // Y.print("Y");
3216  // this->print("X");
3217  // stop();
3218  for (int i = 0; i<_n; ++i) {
3219  if (_X[i] > 1e-20) {
3220  if (prY[i] < 1e-60) {
3221  sum += 1e200;
3222  } else {
3223  sum += _X[i]*log_alt<T>(_X[i]/prY[i]);
3224  }
3225  //sum += _X[i]*log_alt<T>(_X[i]/(prY[i]+1e-100));
3226  }
3227  }
3228  sum += T(-1.0) + Y.sum();
3229  return sum;
3230 };
3231 
3233 template <typename T> inline void Vector<T>::unwhiten(
3234  Vector<T>& meanv, const bool pattern) {
3235  if (pattern) {
3236  const int n =static_cast<int>(sqrt(static_cast<T>(_n)));
3237  int offsetx=0;
3238  for (int j = 0; j<n; ++j) {
3239  offsetx= (offsetx+1) % 2;
3240  int offsety=0;
3241  for (int k = 0; k<n; ++k) {
3242  offsety= (offsety+1) % 2;
3243  _X[j*n+k]+=meanv[2*offsetx+offsety];
3244  }
3245  }
3246  } else {
3247  const int V = meanv.n();
3248  const int sizePatch=_n/V;
3249  for (int j = 0; j<V; ++j) {
3250  T mean = meanv[j];
3251  for (int k = 0; k<sizePatch; ++k) {
3252  _X[sizePatch*j+k]+=mean;
3253  }
3254  }
3255  }
3256 };
3257 
3258 
3260 template <typename T> inline T Vector<T>::mean() {
3261  return this->sum()/_n;
3262 }
3263 
3265 template <typename T> inline T Vector<T>::std() {
3266  T E = this->mean();
3267  T std=0;
3268  for (int i = 0; i<_n; ++i) {
3269  T tmp=_X[i]-E;
3270  std += tmp*tmp;
3271  }
3272  std /= _n;
3273  return sqr_alt<T>(std);
3274 }
3275 
3277 template <typename T> inline void Vector<T>::scal(const T a) {
3278  return cblas_scal<T>(_n,a,_X,1);
3279 };
3280 
3282 template <typename T> inline void Vector<T>::neg() {
3283  for (int i = 0; i<_n; ++i) _X[i]=-_X[i];
3284 };
3285 
3287 template <typename T> inline void Vector<T>::exp() {
3288  vExp<T>(_n,_X,_X);
3289 };
3290 
3292 template <typename T> inline void Vector<T>::log() {
3293  for (int i=0; i<_n; ++i) _X[i]=alt_log<T>(_X[i]);
3294 };
3295 
3297 template <typename T> inline void Vector<T>::logexp() {
3298  for (int i = 0; i<_n; ++i) {
3299  if (_X[i] < -30) {
3300  _X[i]=0;
3301  } else if (_X[i] < 30) {
3302  _X[i]= alt_log<T>( T(1.0) + exp_alt<T>( _X[i] ) );
3303  }
3304  }
3305 };
3306 
3308 template <typename T> inline T Vector<T>::softmax(const int y) {
3309  this->add(-_X[y]);
3310  _X[y]=-INFINITY;
3311  T max=this->maxval();
3312  if (max > 30) {
3313  return max;
3314  } else if (max < -30) {
3315  return 0;
3316  } else {
3317  _X[y]=T(0.0);
3318  this->exp();
3319  return alt_log<T>(this->sum());
3320  }
3321 };
3322 
3324 template <typename T> inline T Vector<T>::asum() const {
3325  return cblas_asum<T>(_n,_X,1);
3326 };
3327 
3328 template <typename T> inline T Vector<T>::lzero() const {
3329  int count=0;
3330  for (int i = 0; i<_n; ++i)
3331  if (_X[i] != 0) ++count;
3332  return count;
3333 };
3334 
3335 
3336 template <typename T> inline T Vector<T>::afused() const {
3337  T sum = 0;
3338  for (int i = 1; i<_n; ++i) {
3339  sum += abs<T>(_X[i]-_X[i-1]);
3340  }
3341  return sum;
3342 }
3344 template <typename T> inline T Vector<T>::sum() const {
3345  T sum=T();
3346  for (int i = 0; i<_n; ++i) sum +=_X[i];
3347  return sum;
3348 };
3349 
3351 template <typename T> inline void Vector<T>::sign(Vector<T>& signs) const {
3352  T* prSign=signs.rawX();
3353  for (int i = 0; i<_n; ++i) {
3354  if (_X[i] == 0) {
3355  prSign[i]=0.0;
3356  } else {
3357  prSign[i] = _X[i] > 0 ? 1.0 : -1.0;
3358  }
3359  }
3360 };
3361 
3364 template <typename T> inline void Vector<T>::l1project(Vector<T>& out,
3365  const T thrs, const bool simplex) const {
3366  out.copy(*this);
3367  if (simplex) {
3368  out.thrsPos();
3369  } else {
3370  vAbs<T>(_n,out._X,out._X);
3371  }
3372  T norm1 = out.sum();
3373  if (norm1 <= thrs) {
3374  if (!simplex) out.copy(*this);
3375  return;
3376  }
3377  T* prU = out._X;
3378  int sizeU = _n;
3379 
3380  T sum = T();
3381  int sum_card = 0;
3382 
3383  while (sizeU > 0) {
3384  // put the pivot in prU[0]
3385  swap(prU[0],prU[sizeU/2]);
3386  T pivot = prU[0];
3387  int sizeG=1;
3388  T sumG=pivot;
3389 
3390  for (int i = 1; i<sizeU; ++i) {
3391  if (prU[i] >= pivot) {
3392  sumG += prU[i];
3393  swap(prU[sizeG++],prU[i]);
3394  }
3395  }
3396 
3397  if (sum + sumG - pivot*(sum_card + sizeG) <= thrs) {
3398  sum_card += sizeG;
3399  sum += sumG;
3400  prU +=sizeG;
3401  sizeU -= sizeG;
3402  } else {
3403  ++prU;
3404  sizeU = sizeG-1;
3405  }
3406  }
3407  T lambda = (sum-thrs)/sum_card;
3408  out.copy(*this);
3409  if (simplex) {
3410  out.thrsPos();
3411  }
3412  out.softThrshold(lambda);
3413 };
3414 
3417 template <typename T> inline void Vector<T>::l1project_weighted(Vector<T>& out, const Vector<T>& weights,
3418  const T thrs, const bool residual) const {
3419  out.copy(*this);
3420  if (thrs==0) {
3421  out.setZeros();
3422  return;
3423  }
3424  vAbs<T>(_n,out._X,out._X);
3425  out.div(weights);
3426  Vector<int> keys(_n);
3427  for (int i = 0; i<_n; ++i) keys[i]=i;
3428  out.sort2(keys,false);
3429  T sum1=0;
3430  T sum2=0;
3431  T lambda=0;
3432  for (int i = 0; i<_n; ++i) {
3433  const T lambda_old=lambda;
3434  const T fact=weights[keys[i]]*weights[keys[i]];
3435  lambda=out[i];
3436  sum2 += fact;
3437  sum1 += fact*lambda;
3438  if (sum1 - lambda*sum2 >= thrs) {
3439  sum2-=fact;
3440  sum1-=fact*lambda;
3441  lambda=lambda_old;
3442  break;
3443  }
3444  }
3445  lambda=MAX(0,(sum1-thrs)/sum2);
3446 
3447  if (residual) {
3448  for (int i = 0; i<_n; ++i) {
3449  out._X[i]=_X[i] > 0 ? MIN(_X[i],lambda*weights[i]) : MAX(_X[i],-lambda*weights[i]);
3450  }
3451  } else {
3452  for (int i = 0; i<_n; ++i) {
3453  out._X[i]=_X[i] > 0 ? MAX(0,_X[i]-lambda*weights[i]) : MIN(0,_X[i]+lambda*weights[i]);
3454  }
3455  }
3456 };
3457 
3458 
3459 template <typename T>
3460 inline void Vector<T>::project_sft_binary(const Vector<T>& y) {
3461  T mean = this->mean();
3462  T thrs=mean;
3463  while (abs(mean) > EPSILON) {
3464  int n_seuils=0;
3465  for (int i = 0; i< _n; ++i) {
3466  _X[i] = _X[i]-thrs;
3467  const T val = y[i]*_X[i];
3468  if (val > 0) {
3469  ++n_seuils;
3470  _X[i]=0;
3471  } else if (val < -1.0) {
3472  ++n_seuils;
3473  _X[i] = -y[i];
3474  }
3475  }
3476  mean = this->mean();
3477  thrs= mean * _n/(_n-n_seuils);
3478  }
3479 };
3480 
3481 template <typename T>
3482 inline void Vector<T>::project_sft(const Vector<int>& labels, const int clas) {
3483  T mean = this->mean();
3484  T thrs=mean;
3485 
3486  while (abs(mean) > EPSILON) {
3487  int n_seuils=0;
3488  for (int i = 0; i< _n; ++i) {
3489  _X[i] = _X[i]-thrs;
3490  if (labels[i]==clas) {
3491  if (_X[i] < -1.0) {
3492  _X[i]=-1.0;
3493  ++n_seuils;
3494  }
3495  } else {
3496  if (_X[i] < 0) {
3497  ++n_seuils;
3498  _X[i]=0;
3499  }
3500  }
3501  }
3502  mean = this->mean();
3503  thrs= mean * _n/(_n-n_seuils);
3504  }
3505 };
3506 
3507 template <typename T>
3508 inline void Vector<T>::sparseProject(Vector<T>& out, const T thrs, const int mode, const T lambda1,
3509  const T lambda2, const T lambda3, const bool pos) {
3510  if (mode == 1) {
3512  this->l1project(out,thrs,pos);
3513  } else if (mode == 2) {
3515  if (lambda1 > 1e-10) {
3516  this->scal(lambda1);
3517  this->l1l2project(out,thrs,2.0/(lambda1*lambda1),pos);
3518  this->scal(T(1.0/lambda1));
3519  out.scal(T(1.0/lambda1));
3520  } else {
3521  out.copy(*this);
3522  out.normalize2();
3523  out.scal(sqrt(thrs));
3524  }
3525  } else if (mode == 3) {
3527  this->l1l2project(out,thrs,lambda1,pos);
3528  } else if (mode == 4) {
3530  out.copy(*this);
3531  if (pos)
3532  out.thrsPos();
3533  out.softThrshold(lambda1);
3534  T nrm=out.nrm2sq();
3535  if (nrm > thrs)
3536  out.scal(sqr_alt<T>(thrs/nrm));
3537  } else if (mode == 5) {
3539  // this->fusedProject(out,lambda1,lambda2,100);
3540  // T nrm=out.nrm2sq();
3541  // if (nrm > thrs)
3542  // out.scal(sqr_alt<T>(thrs/nrm));
3543  // } else if (mode == 6) {
3545  this->fusedProjectHomotopy(out,lambda1,lambda2,lambda3,true);
3546 } else if (mode==6) {
3548  this->fusedProjectHomotopy(out,lambda1/thrs,lambda2/thrs,lambda3/thrs,false);
3549 } else {
3551  if (lambda1 < 1e-10) {
3552  out.copy(*this);
3553  if (pos)
3554  out.thrsPos();
3555  out.normalize2();
3556  out.scal(sqrt(thrs));
3557  } else if (lambda1 > 0.999999) {
3558  this->l1project(out,thrs,pos);
3559  } else {
3560  this->sparseProject(out,thrs/(1.0-lambda1),2,lambda1/(1-lambda1),0,0,pos);
3561  }
3562 }
3563 };
3564 
3566 template <typename T>
3567 inline void Vector<T>::l1l2projectb(Vector<T>& out, const T thrs, const T gamma, const bool pos,
3568  const int mode) {
3569  if (mode == 1) {
3571  this->scal(gamma);
3572  this->l1l2project(out,thrs,2.0/(gamma*gamma),pos);
3573  this->scal(T(1.0/gamma));
3574  out.scal(T(1.0/gamma));
3575  } else if (mode == 2) {
3577  this->l1l2project(out,thrs,gamma,pos);
3578  } else if (mode == 3) {
3580  out.copy(*this);
3581  if (pos)
3582  out.thrsPos();
3583  out.softThrshold(gamma);
3584  T nrm=out.nrm2();
3585  if (nrm > thrs)
3586  out.scal(thrs/nrm);
3587  }
3588 }
3589 
3592 template <typename T>
3593  inline void Vector<T>::l1l2project(Vector<T>& out, const T thrs, const T gamma, const bool pos) const {
3594  if (gamma == 0)
3595  return this->l1project(out,thrs,pos);
3596  out.copy(*this);
3597  if (pos) {
3598  out.thrsPos();
3599  } else {
3600  vAbs<T>(_n,out._X,out._X);
3601  }
3602  T norm = out.sum() + gamma*out.nrm2sq();
3603  if (norm <= thrs) {
3604  if (!pos) out.copy(*this);
3605  return;
3606  }
3607 
3609  T* prU = out._X;
3610  int sizeU = _n;
3611 
3612  T sum = 0;
3613  int sum_card = 0;
3614 
3615  while (sizeU > 0) {
3616  // put the pivot in prU[0]
3617  swap(prU[0],prU[sizeU/2]);
3618  T pivot = prU[0];
3619  int sizeG=1;
3620  T sumG=pivot+0.5*gamma*pivot*pivot;
3621 
3622  for (int i = 1; i<sizeU; ++i) {
3623  if (prU[i] >= pivot) {
3624  sumG += prU[i]+0.5*gamma*prU[i]*prU[i];
3625  swap(prU[sizeG++],prU[i]);
3626  }
3627  }
3628  if (sum + sumG - pivot*(1+0.5*gamma*pivot)*(sum_card + sizeG) <
3629  thrs*(1+gamma*pivot)*(1+gamma*pivot)) {
3630  sum_card += sizeG;
3631  sum += sumG;
3632  prU +=sizeG;
3633  sizeU -= sizeG;
3634  } else {
3635  ++prU;
3636  sizeU = sizeG-1;
3637  }
3638  }
3639  T a = gamma*gamma*thrs+0.5*gamma*sum_card;
3640  T b = 2*gamma*thrs+sum_card;
3641  T c=thrs-sum;
3642  T delta = b*b-4*a*c;
3643  T lambda = (-b+sqrt(delta))/(2*a);
3644 
3645  out.copy(*this);
3646  if (pos) {
3647  out.thrsPos();
3648  }
3649  out.softThrshold(lambda);
3650  out.scal(T(1.0/(1+lambda*gamma)));
3651  };
3652 
3653 template <typename T>
3654 static inline T fusedHomotopyAux(const bool& sign1,
3655  const bool& sign2,
3656  const bool& sign3,
3657  const T& c1,
3658  const T& c2) {
3659  if (sign1) {
3660  if (sign2) {
3661  return sign3 ? 0 : c2;
3662  } else {
3663  return sign3 ? -c2-c1 : -c1;
3664  }
3665  } else {
3666  if (sign2) {
3667  return sign3 ? c1 : c1+c2;
3668  } else {
3669  return sign3 ? -c2 : 0;
3670  }
3671  }
3672 };
3673 
3674 template <typename T>
3675 inline void Vector<T>::fusedProjectHomotopy(Vector<T>& alpha,
3676  const T lambda1,const T lambda2,const T lambda3,
3677  const bool penalty) {
3678  T* pr_DtR=_X;
3679  const int K = _n;
3680  alpha.setZeros();
3681  Vector<T> u(K); // regularization path for gamma
3682  Vector<T> Du(K); // regularization path for alpha
3683  Vector<T> DDu(K); // regularization path for alpha
3684  Vector<T> gamma(K); // auxiliary variable
3685  Vector<T> c(K); // auxiliary variables
3686  Vector<T> scores(K); // auxiliary variables
3687  gamma.setZeros();
3688  T* pr_gamma = gamma.rawX();
3689  T* pr_u = u.rawX();
3690  T* pr_Du = Du.rawX();
3691  T* pr_DDu = DDu.rawX();
3692  T* pr_c = c.rawX();
3693  T* pr_scores = scores.rawX();
3694  Vector<int> ind(K+1);
3695  Vector<bool> signs(K);
3696  ind.set(K);
3697  int* pr_ind = ind.rawX();
3698  bool* pr_signs = signs.rawX();
3699 
3701  T sumBeta = this->sum();
3702 
3704  pr_gamma[0]=sumBeta/K;
3706  alpha.set(pr_gamma[0]);
3708  this->sub(alpha);
3709  for (int j = K-2; j>=0; --j)
3710  pr_DtR[j] += pr_DtR[j+1];
3711 
3712  pr_DtR[0]=0;
3713  pr_ind[0]=0;
3714  pr_signs[0] = pr_DtR[0] > 0;
3715  pr_c[0]=T(1.0)/K;
3716  int currentInd=this->fmax();
3717  T currentLambda=abs<T>(pr_DtR[currentInd]);
3718  bool newAtom = true;
3719 
3721  for (int i = 1; i<K; ++i) {
3724  if (penalty && currentLambda <= lambda2) break;
3725  if (!penalty) {
3727  scores.copy(alpha);
3728  scores.softThrshold(lambda1*currentLambda/lambda2);
3729  scores.scal(T(1.0/(1.0+lambda3*currentLambda/lambda2)));
3730  if (lambda1*scores.asum()+lambda2*scores.afused()+0.5*
3731  lambda3*scores.nrm2sq() >= T(1.0)) break;
3732  }
3733 
3735  if (newAtom) {
3736  int j;
3737  for (j = 1; j<i; ++j)
3738  if (pr_ind[j] > currentInd) break;
3739  for (int k = i; k>j; --k) {
3740  pr_ind[k]=pr_ind[k-1];
3741  pr_c[k]=pr_c[k-1];
3742  pr_signs[k]=pr_signs[k-1];
3743  }
3744  pr_ind[j]=currentInd;
3745  pr_signs[j]=pr_DtR[currentInd] > 0;
3746  pr_c[j-1]=T(1.0)/(pr_ind[j]-pr_ind[j-1]);
3747  pr_c[j]=T(1.0)/(pr_ind[j+1]-pr_ind[j]);
3748  }
3749 
3750  // Compute u
3751  pr_u[0]= pr_signs[1] ? -pr_c[0] : pr_c[0];
3752  if (i == 1) {
3753  pr_u[1]=pr_signs[1] ? pr_c[0]+pr_c[1] : -pr_c[0]-pr_c[1];
3754  } else {
3755  pr_u[1]=pr_signs[1] ? pr_c[0]+pr_c[1] : -pr_c[0]-pr_c[1];
3756  pr_u[1]+=pr_signs[2] ? -pr_c[1] : pr_c[1];
3757  for (int j = 2; j<i; ++j) {
3758  pr_u[j]=2*fusedHomotopyAux<T>(pr_signs[j-1],
3759  pr_signs[j],pr_signs[j+1], pr_c[j-1],pr_c[j]);
3760  }
3761  pr_u[i] = pr_signs[i-1] ? -pr_c[i-1] : pr_c[i-1];
3762  pr_u[i] += pr_signs[i] ? pr_c[i-1]+pr_c[i] : -pr_c[i-1]-pr_c[i];
3763  }
3764 
3765  // Compute Du
3766  pr_Du[0]=pr_u[0];
3767  for (int k = 1; k<pr_ind[1]; ++k)
3768  pr_Du[k]=pr_Du[0];
3769  for (int j = 1; j<=i; ++j) {
3770  pr_Du[pr_ind[j]]=pr_Du[pr_ind[j]-1]+pr_u[j];
3771  for (int k = pr_ind[j]+1; k<pr_ind[j+1]; ++k)
3772  pr_Du[k]=pr_Du[pr_ind[j]];
3773  }
3774 
3776  DDu.copy(Du);
3777  for (int j = K-2; j>=0; --j)
3778  pr_DDu[j] += pr_DDu[j+1];
3779 
3781  T max_step1 = INFINITY;
3782  if (penalty) {
3783  max_step1 = currentLambda-lambda2;
3784  }
3785 
3787  T max_step2 = INFINITY;
3788  int step_out = -1;
3789  for (int j = 1; j<=i; ++j) {
3790  T ratio = -pr_gamma[pr_ind[j]]/pr_u[j];
3791  if (ratio > 0 && ratio <= max_step2) {
3792  max_step2=ratio;
3793  step_out=j;
3794  }
3795  }
3796  T max_step3 = INFINITY;
3798  for (int j = 1; j<K; ++j) {
3799  T sc1 = (currentLambda-pr_DtR[j])/(T(1.0)-pr_DDu[j]);
3800  T sc2 = (currentLambda+pr_DtR[j])/(T(1.0)+pr_DDu[j]);
3801  if (sc1 <= 1e-10) sc1=INFINITY;
3802  if (sc2 <= 1e-10) sc2=INFINITY;
3803  pr_scores[j]= MIN(sc1,sc2);
3804  }
3805  for (int j = 0; j<=i; ++j) {
3806  pr_scores[pr_ind[j]]=INFINITY;
3807  }
3808  currentInd = scores.fmin();
3809  max_step3 = pr_scores[currentInd];
3810  T step = MIN(max_step1,MIN(max_step3,max_step2));
3811  if (step == 0 || step == INFINITY) break;
3812 
3814  for (int j = 0; j<=i; ++j) {
3815  pr_gamma[pr_ind[j]]+=step*pr_u[j];
3816  }
3817  alpha.add(Du,step);
3818  this->add(DDu,-step);
3819  currentLambda -= step;
3820  if (step == max_step2) {
3822  for (int k = step_out; k<=i; ++k)
3823  pr_ind[k]=pr_ind[k+1];
3824  pr_ind[i]=K;
3825  for (int k = step_out; k<=i; ++k)
3826  pr_signs[k]=pr_signs[k+1];
3827  pr_c[step_out-1]=T(1.0)/(pr_ind[step_out]-pr_ind[step_out-1]);
3828  pr_c[step_out]=T(1.0)/(pr_ind[step_out+1]-pr_ind[step_out]);
3829  i-=2;
3830  newAtom=false;
3831  } else {
3832  newAtom=true;
3833  }
3834  }
3835 
3836  if (penalty) {
3837  alpha.softThrshold(lambda1);
3838  alpha.scal(T(1.0/(1.0+lambda3)));
3839  } else {
3840  alpha.softThrshold(lambda1*currentLambda/lambda2);
3841  alpha.scal(T(1.0/(1.0+lambda3*currentLambda/lambda2)));
3842  }
3843 };
3844 
3845 template <typename T>
3846 inline void Vector<T>::fusedProject(Vector<T>& alpha, const T lambda1, const T lambda2,
3847  const int itermax) {
3848  T* pr_alpha= alpha.rawX();
3849  T* pr_beta=_X;
3850  const int K = alpha.n();
3851 
3852  T total_alpha =alpha.sum();
3854  for (int i = K-2; i>=0; --i)
3855  pr_beta[i]+=pr_beta[i+1];
3856 
3857  for (int i = 0; i<itermax; ++i) {
3858  T sum_alpha=0;
3859  T sum_diff = 0;
3861  T gamma_old=pr_alpha[0];
3862  pr_alpha[0]=(K*gamma_old+pr_beta[0]-
3863  total_alpha)/K;
3864  T diff = pr_alpha[0]-gamma_old;
3865  sum_diff += diff;
3866  sum_alpha += pr_alpha[0];
3867  total_alpha +=K*diff;
3868 
3870  for (int j = 1; j<K; ++j) {
3871  pr_alpha[j]+=sum_diff;
3872  T gamma_old=pr_alpha[j]-pr_alpha[j-1];
3873  T gamma_new=softThrs((K-j)*gamma_old+pr_beta[j]-
3874  (total_alpha-sum_alpha),lambda2)/(K-j);
3875  pr_alpha[j]=pr_alpha[j-1]+gamma_new;
3876  T diff = gamma_new-gamma_old;
3877  sum_diff += diff;
3878  sum_alpha+=pr_alpha[j];
3879  total_alpha +=(K-j)*diff;
3880  }
3881  }
3882  alpha.softThrshold(lambda1);
3883 
3884 };
3885 
3887 template <typename T>
3888 inline void Vector<T>::sort(const bool mode) {
3889  if (mode) {
3890  lasrt<T>(incr,_n,_X);
3891  } else {
3892  lasrt<T>(decr,_n,_X);
3893  }
3894 };
3895 
3896 
3898 template <typename T>
3899 inline void Vector<T>::sort(Vector<T>& out, const bool mode) const {
3900  out.copy(*this);
3901  out.sort(mode);
3902 };
3903 
3904 template <typename T>
3905 inline void Vector<T>::sort2(Vector<int>& key, const bool mode) {
3906  quick_sort(key.rawX(),_X,0,_n-1,mode);
3907 };
3908 
3909 
3910 template <typename T>
3911 inline void Vector<T>::sort2(Vector<T>& out, Vector<int>& key, const bool mode) const {
3912  out.copy(*this);
3913  out.sort2(key,mode);
3914 }
3915 
3916 template <typename T>
3917 inline void Vector<T>::applyBayerPattern(const int offset) {
3918  int sizePatch=_n/3;
3919  int n = static_cast<int>(sqrt(static_cast<T>(sizePatch)));
3920  if (offset == 0) {
3921  // R
3922  for (int i = 0; i<n; ++i) {
3923  const int step = (i % 2) ? 1 : 2;
3924  const int off = (i % 2) ? 0 : 1;
3925  for (int j = off; j<n; j+=step) {
3926  _X[i*n+j]=0;
3927  }
3928  }
3929  // G
3930  for (int i = 0; i<n; ++i) {
3931  const int step = 2;
3932  const int off = (i % 2) ? 1 : 0;
3933  for (int j = off; j<n; j+=step) {
3934  _X[sizePatch+i*n+j]=0;
3935  }
3936  }
3937  // B
3938  for (int i = 0; i<n; ++i) {
3939  const int step = (i % 2) ? 2 : 1;
3940  const int off = 0;
3941  for (int j = off; j<n; j+=step) {
3942  _X[2*sizePatch+i*n+j]=0;
3943  }
3944  }
3945  } else if (offset == 1) {
3946  // R
3947  for (int i = 0; i<n; ++i) {
3948  const int step = (i % 2) ? 2 : 1;
3949  const int off = (i % 2) ? 1 : 0;
3950  for (int j = off; j<n; j+=step) {
3951  _X[i*n+j]=0;
3952  }
3953  }
3954  // G
3955  for (int i = 0; i<n; ++i) {
3956  const int step = 2;
3957  const int off = (i % 2) ? 0 : 1;
3958  for (int j = off; j<n; j+=step) {
3959  _X[sizePatch+i*n+j]=0;
3960  }
3961  }
3962  // B
3963  for (int i = 0; i<n; ++i) {
3964  const int step = (i % 2) ? 1 : 2;
3965  const int off = 0;
3966  for (int j = off; j<n; j+=step) {
3967  _X[2*sizePatch+i*n+j]=0;
3968  }
3969  }
3970  } else if (offset == 2) {
3971  // R
3972  for (int i = 0; i<n; ++i) {
3973  const int step = (i % 2) ? 1 : 2;
3974  const int off = 0;
3975  for (int j = off; j<n; j+=step) {
3976  _X[i*n+j]=0;
3977  }
3978  }
3979  // G
3980  for (int i = 0; i<n; ++i) {
3981  const int step = 2;
3982  const int off = (i % 2) ? 0 : 1;
3983  for (int j = off; j<n; j+=step) {
3984  _X[sizePatch+i*n+j]=0;
3985  }
3986  }
3987  // B
3988  for (int i = 0; i<n; ++i) {
3989  const int step = (i % 2) ? 2 : 1;
3990  const int off = (i % 2) ? 1 : 0;
3991  for (int j = off; j<n; j+=step) {
3992  _X[2*sizePatch+i*n+j]=0;
3993  }
3994  }
3995  } else if (offset == 3) {
3996  // R
3997  for (int i = 0; i<n; ++i) {
3998  const int step = (i % 2) ? 2 : 1;
3999  const int off = 0;
4000  for (int j = off; j<n; j+=step) {
4001  _X[i*n+j]=0;
4002  }
4003  }
4004  // G
4005  for (int i = 0; i<n; ++i) {
4006  const int step = 2;
4007  const int off = (i % 2) ? 1 : 0;
4008  for (int j = off; j<n; j+=step) {
4009  _X[sizePatch+i*n+j]=0;
4010  }
4011  }
4012  // B
4013  for (int i = 0; i<n; ++i) {
4014  const int step = (i % 2) ? 1 : 2;
4015  const int off = (i % 2) ? 0 : 1;
4016  for (int j = off; j<n; j+=step) {
4017  _X[2*sizePatch+i*n+j]=0;
4018  }
4019  }
4020  }
4021 };
4022 
4023 
4025 template <typename T> inline void Vector<T>::toSparse(
4026  SpVector<T>& vec) const {
4027  int L=0;
4028  T* v = vec._v;
4029  int* r = vec._r;
4030  for (int i = 0; i<_n; ++i) {
4031  if (_X[i] != T()) {
4032  v[L]=_X[i];
4033  r[L++]=i;
4034  }
4035  }
4036  vec._L=L;
4037 };
4038 
4039 
4040 template <typename T>
4041 inline void Vector<T>::copyMask(Vector<T>& out, Vector<bool>& mask) const {
4042  out.resize(_n);
4043  int pointer=0;
4044  for (int i = 0; i<_n; ++i) {
4045  if (mask[i])
4046  out[pointer++]=_X[i];
4047  }
4048  out.setn(pointer);
4049 };
4050 
4051 template <typename T>
4052 inline void Matrix<T>::copyMask(Matrix<T>& out, Vector<bool>& mask) const {
4053  out.resize(_m,_n);
4054  int count=0;
4055  for (int i = 0; i<mask.n(); ++i)
4056  if (mask[i])
4057  ++count;
4058  out.setm(count);
4059  for (int i = 0; i<_n; ++i) {
4060  int pointer=0;
4061  for (int j = 0; j<_m; ++j) {
4062  if (mask[j]) {
4063  out[i*count+pointer]=_X[i*_m+j];
4064  ++pointer;
4065  }
4066  }
4067  }
4068 };
4069 
4070 
4071 
4072 /* ****************************
4073  * Implementation of SpMatrix
4074  * ****************************/
4075 
4076 
4078 template <typename T> SpMatrix<T>::SpMatrix(T* v, int* r, int* pB, int* pE,
4079  int m, int n, int nzmax) :
4080  _externAlloc(true), _v(v), _r(r), _pB(pB), _pE(pE), _m(m), _n(n), _nzmax(nzmax)
4081 { };
4082 
4084 template <typename T> SpMatrix<T>::SpMatrix(int m, int n, int nzmax) :
4085  _externAlloc(false), _m(m), _n(n), _nzmax(nzmax) {
4086 #pragma omp critical
4087  {
4088  _v=new T[nzmax];
4089  _r=new int[nzmax];
4090  _pB=new int[_n+1];
4091  }
4092  _pE=_pB+1;
4093  };
4094 
4096 template <typename T> SpMatrix<T>::SpMatrix() :
4097  _externAlloc(true), _v(NULL), _r(NULL), _pB(NULL), _pE(NULL),
4098  _m(0),_n(0),_nzmax(0) { };
4099 
4100 
4101 template <typename T>
4102 inline void SpMatrix<T>::copy(const SpMatrix<T>& mat) {
4103  this->resize(mat._m,mat._n,mat._nzmax);
4104  memcpy(_v,mat._v,_nzmax*sizeof(T));
4105  memcpy(_r,mat._r,_nzmax*sizeof(int));
4106  memcpy(_pB,mat._pB,(_n+1)*sizeof(int));
4107 }
4108 
4109 
4111 template <typename T> SpMatrix<T>::~SpMatrix() {
4112  clear();
4113 };
4114 
4116 template <typename T> inline void SpMatrix<T>::refCol(int i,
4117  SpVector<T>& vec) const {
4118  if (vec._nzmax > 0) vec.clear();
4119  vec._v=_v+_pB[i];
4120  vec._r=_r+_pB[i];
4121  vec._externAlloc=true;
4122  vec._L=_pE[i]-_pB[i];
4123  vec._nzmax=vec._L;
4124 };
4125 
4127 template<typename T> inline void SpMatrix<T>::print(const string& name) const {
4128  cerr << name << endl;
4129  cerr << _m << " x " << _n << " , " << _nzmax << endl;
4130  for (int i = 0; i<_n; ++i) {
4131  for (int j = _pB[i]; j<_pE[i]; ++j) {
4132  cerr << "(" <<_r[j] << "," << i << ") = " << _v[j] << endl;
4133  }
4134  }
4135 };
4136 
4137 template<typename T>
4138 inline T SpMatrix<T>::operator[](const int index) const {
4139  const int num_col=(index/_m);
4140  const int num_row=index -num_col*_m;
4141  T val = 0;
4142  for (int j = _pB[num_col]; j<_pB[num_col+1]; ++j) {
4143  if (_r[j]==num_row) {
4144  val=_v[j];
4145  break;
4146  }
4147  }
4148  return val;
4149 };
4150 template<typename T>
4151 void SpMatrix<T>::getData(Vector<T>& data, const int index) const {
4152  data.resize(_m);
4153  data.setZeros();
4154  for (int i = _pB[index]; i< _pB[index+1]; ++i)
4155  data[_r[i]]=_v[i];
4156 };
4157 
4158 template<typename T>
4159 void SpMatrix<T>::getGroup(Matrix<T>& data, const vector_groups& groups, const int i) const {
4160  const group& gr = groups[i];
4161  const int N = gr.size();
4162  data.resize(_m,N);
4163  int count=0;
4164  Vector<T> col;
4165  for (group::const_iterator it = gr.begin(); it != gr.end(); ++it) {
4166  data.refCol(count,col);
4167  this->getData(col,*it);
4168  ++count;
4169  }
4170 };
4171 
4173 template <typename T> inline T SpMatrix<T>::asum() const {
4174  return cblas_asum<T>(_pB[_n],_v,1);
4175 };
4176 
4178 template <typename T> inline T SpMatrix<T>::normFsq() const {
4179  return cblas_dot<T>(_pB[_n],_v,1,_v,1);
4180 };
4181 
4182 template <typename T>
4183 inline void SpMatrix<T>::add_direct(const SpMatrix<T>& mat, const T a) {
4184  Vector<T> v2(mat._v,mat._nzmax);
4185  Vector<T> v1(_v,_nzmax);
4186  v1.add(v2,a);
4187 }
4188 
4189 template <typename T>
4190 inline void SpMatrix<T>::copy_direct(const SpMatrix<T>& mat) {
4191  Vector<T> v2(mat._v,_pB[_n]);
4192  Vector<T> v1(_v,_pB[_n]);
4193  v1.copy(v2);
4194 }
4195 
4196 template <typename T>
4197 inline T SpMatrix<T>::dot_direct(const SpMatrix<T>& mat) const {
4198  Vector<T> v2(mat._v,_pB[_n]);
4199  Vector<T> v1(_v,_pB[_n]);
4200  return v1.dot(v2);
4201 }
4202 
4204 template <typename T> inline void SpMatrix<T>::clear() {
4205  if (!_externAlloc) {
4206  delete[](_r);
4207  delete[](_v);
4208  delete[](_pB);
4209  }
4210  _n=0;
4211  _m=0;
4212  _nzmax=0;
4213  _v=NULL;
4214  _r=NULL;
4215  _pB=NULL;
4216  _pE=NULL;
4217  _externAlloc=true;
4218 };
4219 
4221 template <typename T> inline void SpMatrix<T>::resize(const int m,
4222  const int n, const int nzmax) {
4223  if (n == _n && m == _m && nzmax == _nzmax) return;
4224  this->clear();
4225  _n=n;
4226  _m=m;
4227  _nzmax=nzmax;
4228  _externAlloc=false;
4229 #pragma omp critical
4230  {
4231  _v = new T[nzmax];
4232  _r = new int[nzmax];
4233  _pB = new int[_n+1];
4234  }
4235  _pE = _pB+1;
4236  for (int i = 0; i<=_n; ++i) _pB[i]=0;
4237 };
4238 
4240 template <typename T> inline void SpMatrix<T>::scal(const T a) const {
4241  cblas_scal<T>(_pB[_n],a,_v,1);
4242 };
4243 
4245 template <typename T>
4246 inline void SpMatrix<T>::multTrans(const Vector<T>& x, Vector<T>& y,
4247  const T alpha, const T beta) const {
4248  y.resize(_n);
4249  if (beta) {
4250  y.scal(beta);
4251  } else {
4252  y.setZeros();
4253  }
4254  const T* prX = x.rawX();
4255  for (int i = 0; i<_n; ++i) {
4256  T sum=T();
4257  for (int j = _pB[i]; j<_pE[i]; ++j) {
4258  sum+=_v[j]*prX[_r[j]];
4259  }
4260  y[i] += alpha*sum;
4261  }
4262 };
4263 
4265 template <typename T>
4266 inline void SpMatrix<T>::multTrans(const SpVector<T>& x, Vector<T>& y,
4267  const T alpha, const T beta) const {
4268  y.resize(_n);
4269  if (beta) {
4270  y.scal(beta);
4271  } else {
4272  y.setZeros();
4273  }
4274  T* prY = y.rawX();
4275  SpVector<T> col;
4276  for (int i = 0; i<_n; ++i) {
4277  this->refCol(i,col);
4278  prY[i] += alpha*x.dot(col);
4279  }
4280 };
4281 
4282 
4284 template <typename T>
4285 inline void SpMatrix<T>::mult(const Vector<T>& x, Vector<T>& y,
4286  const T alpha, const T beta) const {
4287  y.resize(_m);
4288  if (beta) {
4289  y.scal(beta);
4290  } else {
4291  y.setZeros();
4292  }
4293  const T* prX = x.rawX();
4294  for (int i = 0; i<_n; ++i) {
4295  T sca=alpha* prX[i];
4296  for (int j = _pB[i]; j<_pE[i]; ++j) {
4297  y[_r[j]] += sca*_v[j];
4298  }
4299  }
4300 };
4301 
4302 
4304 template <typename T>
4305 inline void SpMatrix<T>::mult(const SpVector<T>& x, Vector<T>& y,
4306  const T alpha, const T beta) const {
4307  y.resize(_m);
4308  if (beta) {
4309  y.scal(beta);
4310  } else {
4311  y.setZeros();
4312  }
4313  T* prY = y.rawX();
4314  for (int i = 0; i<x.L(); ++i) {
4315  int ind=x.r(i);
4316  T val = alpha * x.v(i);
4317  for (int j = _pB[ind]; j<_pE[ind]; ++j) {
4318  prY[_r[j]] += val *_v[j];
4319  }
4320  }
4321 };
4322 
4324 template <typename T>
4325 inline void SpMatrix<T>::mult(const Matrix<T>& B, Matrix<T>& C,
4326  const bool transA, const bool transB,
4327  const T a, const T b) const {
4328  if (transA) {
4329  if (transB) {
4330  C.resize(_n,B.m());
4331  if (b) {
4332  C.scal(b);
4333  } else {
4334  C.setZeros();
4335  }
4336  SpVector<T> tmp;
4337  Vector<T> row(B.m());
4338  for (int i = 0; i<_n; ++i) {
4339  this->refCol(i,tmp);
4340  B.mult(tmp,row);
4341  C.addRow(i,row,a);
4342  }
4343  } else {
4344  C.resize(_n,B.n());
4345  if (b) {
4346  C.scal(b);
4347  } else {
4348  C.setZeros();
4349  }
4350  SpVector<T> tmp;
4351  Vector<T> row(B.n());
4352  for (int i = 0; i<_n; ++i) {
4353  this->refCol(i,tmp);
4354  B.multTrans(tmp,row);
4355  C.addRow(i,row,a);
4356  }
4357  }
4358  } else {
4359  if (transB) {
4360  C.resize(_m,B.m());
4361  if (b) {
4362  C.scal(b);
4363  } else {
4364  C.setZeros();
4365  }
4366  Vector<T> row(B.n());
4367  Vector<T> col;
4368  for (int i = 0; i<B.m(); ++i) {
4369  B.copyRow(i,row);
4370  C.refCol(i,col);
4371  this->mult(row,col,a,T(1.0));
4372  }
4373  } else {
4374  C.resize(_m,B.n());
4375  if (b) {
4376  C.scal(b);
4377  } else {
4378  C.setZeros();
4379  }
4380  Vector<T> colB;
4381  Vector<T> colC;
4382  for (int i = 0; i<B.n(); ++i) {
4383  B.refCol(i,colB);
4384  C.refCol(i,colC);
4385  this->mult(colB,colC,a,T(1.0));
4386  }
4387  }
4388  }
4389 };
4390 
4392 template <typename T>
4393 inline void SpMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
4394  const bool transA, const bool transB,
4395  const T a, const T b) const {
4396  if (transA) {
4397  if (transB) {
4398  C.resize(_n,B.m());
4399  if (b) {
4400  C.scal(b);
4401  } else {
4402  C.setZeros();
4403  }
4404  SpVector<T> tmp;
4405  Vector<T> row(B.m());
4406  for (int i = 0; i<_n; ++i) {
4407  this->refCol(i,tmp);
4408  B.mult(tmp,row);
4409  C.addRow(i,row,a);
4410  }
4411  } else {
4412  C.resize(_n,B.n());
4413  if (b) {
4414  C.scal(b);
4415  } else {
4416  C.setZeros();
4417  }
4418  SpVector<T> tmp;
4419  Vector<T> row(B.n());
4420  for (int i = 0; i<_n; ++i) {
4421  this->refCol(i,tmp);
4422  B.multTrans(tmp,row);
4423  C.addRow(i,row,a);
4424  }
4425  }
4426  } else {
4427  if (transB) {
4428  C.resize(_m,B.m());
4429  if (b) {
4430  C.scal(b);
4431  } else {
4432  C.setZeros();
4433  }
4434  SpVector<T> colB;
4435  SpVector<T> colA;
4436  for (int i = 0; i<_n; ++i) {
4437  this->refCol(i,colA);
4438  B.refCol(i,colB);
4439  C.rank1Update(colA,colB,a);
4440  }
4441  } else {
4442  C.resize(_m,B.n());
4443  if (b) {
4444  C.scal(b);
4445  } else {
4446  C.setZeros();
4447  }
4448  SpVector<T> colB;
4449  Vector<T> colC;
4450  for (int i = 0; i<B.n(); ++i) {
4451  B.refCol(i,colB);
4452  C.refCol(i,colC);
4453  this->mult(colB,colC,a);
4454  }
4455  }
4456  }
4457 };
4458 
4460 template <typename T>
4461 inline void SpMatrix<T>::multSwitch(const Matrix<T>& B, Matrix<T>& C,
4462  const bool transA, const bool transB,
4463  const T a, const T b) const {
4464  B.mult(*this,C,transB,transA,a,b);
4465 };
4466 
4467 template <typename T>
4468 inline T SpMatrix<T>::dot(const Matrix<T>& x) const {
4469  T sum=0;
4470  for (int i = 0; i<_n; ++i)
4471  for (int j = _pB[i]; j<_pE[i]; ++j) {
4472  sum+=_v[j]*x(_r[j],j);
4473  }
4474  return sum;
4475 };
4476 
4477 
4478 template <typename T>
4479 inline void SpMatrix<T>::copyRow(const int ind, Vector<T>& x) const {
4480  x.resize(_n);
4481  x.setZeros();
4482  for (int i = 0; i<_n; ++i) {
4483  for (int j = _pB[i]; j<_pE[i]; ++j) {
4484  if (_r[j]==ind) {
4485  x[i]=_v[j];
4486  } else if (_r[j] > ind) {
4487  break;
4488  }
4489  }
4490  }
4491 };
4492 
4493 template <typename T>
4494 inline void SpMatrix<T>::addVecToCols(
4495  const Vector<T>& vec, const T a) {
4496  const T* pr_vec = vec.rawX();
4497  if (isEqual(a,T(1.0))) {
4498  for (int i = 0; i<_n; ++i)
4499  for (int j = _pB[i]; j<_pE[i]; ++j)
4500  _v[j] += pr_vec[_r[j]];
4501  } else {
4502  for (int i = 0; i<_n; ++i)
4503  for (int j = _pB[i]; j<_pE[i]; ++j)
4504  _v[j] += a*pr_vec[_r[j]];
4505  }
4506 };
4507 
4508 template <typename T>
4509 inline void SpMatrix<T>::addVecToColsWeighted(
4510  const Vector<T>& vec, const T* weights, const T a) {
4511  const T* pr_vec = vec.rawX();
4512  if (isEqual(a,T(1.0))) {
4513  for (int i = 0; i<_n; ++i)
4514  for (int j = _pB[i]; j<_pE[i]; ++j)
4515  _v[j] += pr_vec[_r[j]]*weights[j-_pB[i]];
4516  } else {
4517  for (int i = 0; i<_n; ++i)
4518  for (int j = _pB[i]; j<_pE[i]; ++j)
4519  _v[j] += a*pr_vec[_r[j]]*weights[j-_pB[i]];
4520  }
4521 };
4522 
4523 template <typename T>
4524 inline void SpMatrix<T>::sum_cols(Vector<T>& sum) const {
4525  sum.resize(_m);
4526  sum.setZeros();
4527  SpVector<T> tmp;
4528  for (int i = 0; i<_n; ++i) {
4529  this->refCol(i,tmp);
4530  sum.add(tmp);
4531  }
4532 };
4533 
4535 template <typename T> inline void SpMatrix<T>::AAt(Matrix<T>& aat) const {
4536  int i,j,k;
4537  int K=_m;
4538  int M=_n;
4539 
4540  /* compute alpha alpha^T */
4541  aat.resize(K,K);
4542  int NUM_THREADS=init_omp(MAX_THREADS);
4543  T* aatT=new T[NUM_THREADS*K*K];
4544  for (j = 0; j<NUM_THREADS*K*K; ++j) aatT[j]=T();
4545 
4546 #pragma omp parallel for private(i,j,k)
4547  for (i = 0; i<M; ++i) {
4548 #ifdef _OPENMP
4549  int numT=omp_get_thread_num();
4550 #else
4551  int numT=0;
4552 #endif
4553  T* write_area=aatT+numT*K*K;
4554  for (j = _pB[i]; j<_pE[i]; ++j) {
4555  for (k = _pB[i]; k<=j; ++k) {
4556  write_area[_r[j]*K+_r[k]]+=_v[j]*_v[k];
4557  }
4558  }
4559  }
4560 
4561  cblas_copy<T>(K*K,aatT,1,aat._X,1);
4562  for (i = 1; i<NUM_THREADS; ++i)
4563  cblas_axpy<T>(K*K,1.0,aatT+K*K*i,1,aat._X,1);
4564  aat.fillSymmetric();
4565  delete[](aatT);
4566 }
4567 
4568 template <typename T>
4569 inline void SpMatrix<T>::XtX(Matrix<T>& XtX) const {
4570  XtX.resize(_n,_n);
4571  XtX.setZeros();
4572  SpVector<T> col;
4573  Vector<T> col_out;
4574  for (int i = 0; i<_n; ++i) {
4575  this->refCol(i,col);
4576  XtX.refCol(i,col_out);
4577  this->multTrans(col,col_out);
4578  }
4579 };
4580 
4581 
4583 template <typename T> inline void SpMatrix<T>::AAt(Matrix<T>& aat,
4584  const Vector<int>& indices) const {
4585  int i,j,k;
4586  int K=_m;
4587  int M=indices.n();
4588 
4589  /* compute alpha alpha^T */
4590  aat.resize(K,K);
4591  int NUM_THREADS=init_omp(MAX_THREADS);
4592  T* aatT=new T[NUM_THREADS*K*K];
4593  for (j = 0; j<NUM_THREADS*K*K; ++j) aatT[j]=T();
4594 
4595 #pragma omp parallel for private(i,j,k)
4596  for (i = 0; i<M; ++i) {
4597  int ii = indices[i];
4598 #ifdef _OPENMP
4599  int numT=omp_get_thread_num();
4600 #else
4601  int numT=0;
4602 #endif
4603  T* write_area=aatT+numT*K*K;
4604  for (j = _pB[ii]; j<_pE[ii]; ++j) {
4605  for (k = _pB[ii]; k<=j; ++k) {
4606  write_area[_r[j]*K+_r[k]]+=_v[j]*_v[k];
4607  }
4608  }
4609  }
4610 
4611  cblas_copy<T>(K*K,aatT,1,aat._X,1);
4612  for (i = 1; i<NUM_THREADS; ++i)
4613  cblas_axpy<T>(K*K,1.0,aatT+K*K*i,1,aat._X,1);
4614  aat.fillSymmetric();
4615  delete[](aatT);
4616 }
4617 
4619 template <typename T> inline void SpMatrix<T>::wAAt(const Vector<T>& w,
4620  Matrix<T>& aat) const {
4621  int i,j,k;
4622  int K=_m;
4623  int M=_n;
4624 
4625  /* compute alpha alpha^T */
4626  aat.resize(K,K);
4627  int NUM_THREADS=init_omp(MAX_THREADS);
4628  T* aatT=new T[NUM_THREADS*K*K];
4629  for (j = 0; j<NUM_THREADS*K*K; ++j) aatT[j]=T();
4630 
4631 #pragma omp parallel for private(i,j,k)
4632  for (i = 0; i<M; ++i) {
4633 #ifdef _OPENMP
4634  int numT=omp_get_thread_num();
4635 #else
4636  int numT=0;
4637 #endif
4638  T* write_area=aatT+numT*K*K;
4639  for (j = _pB[i]; j<_pE[i]; ++j) {
4640  for (k = _pB[i]; k<=j; ++k) {
4641  write_area[_r[j]*K+_r[k]]+=w._X[i]*_v[j]*_v[k];
4642  }
4643  }
4644  }
4645 
4646  cblas_copy<T>(K*K,aatT,1,aat._X,1);
4647  for (i = 1; i<NUM_THREADS; ++i)
4648  cblas_axpy<T>(K*K,1.0,aatT+K*K*i,1,aat._X,1);
4649  aat.fillSymmetric();
4650  delete[](aatT);
4651 }
4652 
4654 template <typename T> inline void SpMatrix<T>::XAt(const Matrix<T>& X,
4655  Matrix<T>& XAt) const {
4656  int j,i;
4657  int n=X._m;
4658  int K=_m;
4659  int M=_n;
4660 
4661  XAt.resize(n,K);
4662  /* compute X alpha^T */
4663  int NUM_THREADS=init_omp(MAX_THREADS);
4664  T* XatT=new T[NUM_THREADS*n*K];
4665  for (j = 0; j<NUM_THREADS*n*K; ++j) XatT[j]=T();
4666 
4667 #pragma omp parallel for private(i,j)
4668  for (i = 0; i<M; ++i) {
4669 #ifdef _OPENMP
4670  int numT=omp_get_thread_num();
4671 #else
4672  int numT=0;
4673 #endif
4674  T* write_area=XatT+numT*n*K;
4675  for (j = _pB[i]; j<_pE[i]; ++j) {
4676  cblas_axpy<T>(n,_v[j],X._X+i*n,1,write_area+_r[j]*n,1);
4677  }
4678  }
4679 
4680  cblas_copy<T>(n*K,XatT,1,XAt._X,1);
4681  for (i = 1; i<NUM_THREADS; ++i)
4682  cblas_axpy<T>(n*K,1.0,XatT+n*K*i,1,XAt._X,1);
4683  delete[](XatT);
4684 };
4685 
4687 template <typename T> inline void SpMatrix<T>::XAt(const Matrix<T>& X,
4688  Matrix<T>& XAt, const Vector<int>& indices) const {
4689  int j,i;
4690  int n=X._m;
4691  int K=_m;
4692  int M=indices.n();
4693 
4694  XAt.resize(n,K);
4695  /* compute X alpha^T */
4696  int NUM_THREADS=init_omp(MAX_THREADS);
4697  T* XatT=new T[NUM_THREADS*n*K];
4698  for (j = 0; j<NUM_THREADS*n*K; ++j) XatT[j]=T();
4699 
4700 #pragma omp parallel for private(i,j)
4701  for (i = 0; i<M; ++i) {
4702  int ii = indices[i];
4703 #ifdef _OPENMP
4704  int numT=omp_get_thread_num();
4705 #else
4706  int numT=0;
4707 #endif
4708  T* write_area=XatT+numT*n*K;
4709  for (j = _pB[ii]; j<_pE[ii]; ++j) {
4710  cblas_axpy<T>(n,_v[j],X._X+i*n,1,write_area+_r[j]*n,1);
4711  }
4712  }
4713 
4714  cblas_copy<T>(n*K,XatT,1,XAt._X,1);
4715  for (i = 1; i<NUM_THREADS; ++i)
4716  cblas_axpy<T>(n*K,1.0,XatT+n*K*i,1,XAt._X,1);
4717  delete[](XatT);
4718 };
4719 
4721 template <typename T> inline void SpMatrix<T>::wXAt(const Vector<T>& w,
4722  const Matrix<T>& X, Matrix<T>& XAt, const int numThreads) const {
4723  int j,l,i;
4724  int n=X._m;
4725  int K=_m;
4726  int M=_n;
4727  int Mx = X._n;
4728  int numRepX= M/Mx;
4729  assert(numRepX*Mx == M);
4730  XAt.resize(n,K);
4731  /* compute X alpha^T */
4732  int NUM_THREADS=init_omp(numThreads);
4733  T* XatT=new T[NUM_THREADS*n*K];
4734  for (j = 0; j<NUM_THREADS*n*K; ++j) XatT[j]=T();
4735 
4736 #pragma omp parallel for private(i,j,l)
4737  for (i = 0; i<Mx; ++i) {
4738 #ifdef _OPENMP
4739  int numT=omp_get_thread_num();
4740 #else
4741  int numT=0;
4742 #endif
4743  T * write_area=XatT+numT*n*K;
4744  for (l = 0; l<numRepX; ++l) {
4745  int ind=numRepX*i+l;
4746  if (w._X[ind] != 0)
4747  for (j = _pB[ind]; j<_pE[ind]; ++j) {
4748  cblas_axpy<T>(n,w._X[ind]*_v[j],X._X+i*n,1,write_area+_r[j]*n,1);
4749  }
4750  }
4751  }
4752 
4753  cblas_copy<T>(n*K,XatT,1,XAt._X,1);
4754  for (i = 1; i<NUM_THREADS; ++i)
4755  cblas_axpy<T>(n*K,1.0,XatT+n*K*i,1,XAt._X,1);
4756  delete[](XatT);
4757 };
4758 
4760 template<typename T> inline void SpMatrix<T>::toFull(Matrix<T>& matrix) const {
4761  matrix.resize(_m,_n);
4762  matrix.setZeros();
4763  T* out = matrix._X;
4764  for (int i=0; i<_n; ++i) {
4765  for (int j = _pB[i]; j<_pE[i]; ++j) {
4766  out[i*_m+_r[j]]=_v[j];
4767  }
4768  }
4769 };
4770 
4772 template <typename T> inline void SpMatrix<T>::toFullTrans(
4773  Matrix<T>& matrix) const {
4774  matrix.resize(_n,_m);
4775  matrix.setZeros();
4776  T* out = matrix._X;
4777  for (int i=0; i<_n; ++i) {
4778  for (int j = _pB[i]; j<_pE[i]; ++j) {
4779  out[i+_r[j]*_n]=_v[j];
4780  }
4781  }
4782 };
4783 
4784 
4786 template <typename T> inline void SpMatrix<T>::convert(const Matrix<T>&vM,
4787  const Matrix<int>& rM, const int K) {
4788  const int M = rM.n();
4789  const int L = rM.m();
4790  const int* r = rM.X();
4791  const T* v = vM.X();
4792  int count=0;
4793  for (int i = 0; i<M*L; ++i) if (r[i] != -1) ++count;
4794  resize(K,M,count);
4795  count=0;
4796  for (int i = 0; i<M; ++i) {
4797  _pB[i]=count;
4798  for (int j = 0; j<L; ++j) {
4799  if (r[i*L+j] == -1) break;
4800  _v[count]=v[i*L+j];
4801  _r[count++]=r[i*L+j];
4802  }
4803  _pE[i]=count;
4804  }
4805  for (int i = 0; i<M; ++i) sort(_r,_v,_pB[i],_pE[i]-1);
4806 };
4807 
4809 template <typename T> inline void SpMatrix<T>::convert2(
4810  const Matrix<T>&vM, const Vector<int>& rv, const int K) {
4811  const int M = vM.n();
4812  const int L = vM.m();
4813  int* r = rv.rawX();
4814  const T* v = vM.X();
4815  int LL=0;
4816  for (int i = 0; i<L; ++i) if (r[i] != -1) ++LL;
4817  this->resize(K,M,LL*M);
4818  int count=0;
4819  for (int i = 0; i<M; ++i) {
4820  _pB[i]=count;
4821  for (int j = 0; j<LL; ++j) {
4822  _v[count]=v[i*L+j];
4823  _r[count++]=r[j];
4824  }
4825  _pE[i]=count;
4826  }
4827  for (int i = 0; i<M; ++i) sort(_r,_v,_pB[i],_pE[i]-1);
4828 };
4829 
4831 template <typename T>
4832 inline void SpMatrix<T>::norm_2sq_cols(Vector<T>& norms) const {
4833  norms.resize(_n);
4834  SpVector<T> col;
4835  for (int i = 0; i<_n; ++i) {
4836  this->refCol(i,col);
4837  norms[i] = col.nrm2sq();
4838  }
4839 };
4840 
4841 template <typename T>
4842 inline void SpMatrix<T>::norm_0_cols(Vector<T>& norms) const {
4843  norms.resize(_n);
4844  SpVector<T> col;
4845  for (int i = 0; i<_n; ++i) {
4846  this->refCol(i,col);
4847  norms[i] = static_cast<T>(col.length());
4848  }
4849 };
4850 
4851 template <typename T>
4852 inline void SpMatrix<T>::norm_1_cols(Vector<T>& norms) const {
4853  norms.resize(_n);
4854  SpVector<T> col;
4855  for (int i = 0; i<_n; ++i) {
4856  this->refCol(i,col);
4857  norms[i] =col.asum();
4858  }
4859 };
4860 
4861 
4862 /* ***************************
4863  * Implementation of SpVector
4864  * ***************************/
4865 
4866 
4868 template <typename T> SpVector<T>::SpVector(T* v, int* r, int L, int nzmax) :
4869  _externAlloc(true), _v(v), _r(r), _L(L), _nzmax(nzmax) { };
4870 
4872 template <typename T> SpVector<T>::SpVector(int nzmax) :
4873  _externAlloc(false), _L(0), _nzmax(nzmax) {
4874 #pragma omp critical
4875  {
4876  _v = new T[nzmax];
4877  _r = new int[nzmax];
4878  }
4879  };
4880 
4882 template <typename T> SpVector<T>::SpVector() : _externAlloc(true), _v(NULL), _r(NULL), _L(0),
4883  _nzmax(0) { };
4884 
4885 
4887 template <typename T> SpVector<T>::~SpVector() { clear(); };
4888 
4889 
4891 template <typename T> inline T SpVector<T>::asum() const {
4892  return cblas_asum<T>(_L,_v,1);
4893 };
4894 
4896 template <typename T> inline T SpVector<T>::nrm2sq() const {
4897  return cblas_dot<T>(_L,_v,1,_v,1);
4898 };
4899 
4901 template <typename T> inline T SpVector<T>::nrm2() const {
4902  return cblas_nrm2<T>(_L,_v,1);
4903 };
4904 
4906 template <typename T> inline T SpVector<T>::fmaxval() const {
4907  Vector<T> tmp(_v,_L);
4908  return tmp.fmaxval();
4909 };
4910 
4912 template <typename T> inline void SpVector<T>::print(const string& name) const {
4913  std::cerr << name << std::endl;
4914  std::cerr << _nzmax << std::endl;
4915  for (int i = 0; i<_L; ++i)
4916  cerr << "(" <<_r[i] << ", " << _v[i] << ")" << endl;
4917 };
4918 
4920 template <typename T> inline void SpVector<T>::refIndices(
4921  Vector<int>& indices) const {
4922  indices.setPointer(_r,_L);
4923 };
4924 
4926 template <typename T> inline void SpVector<T>::refVal(
4927  Vector<T>& val) const {
4928  val.setPointer(_v,_L);
4929 };
4930 
4932 template <typename T> inline void SpVector<T>::sqr() {
4933  vSqr<T>(_L,_v,_v);
4934 };
4935 
4936 template <typename T>
4937 inline T SpVector<T>::dot(const SpVector<T>& vec) const {
4938  T sum=T();
4939  int countI = 0;
4940  int countJ = 0;
4941  while (countI < _L && countJ < vec._L) {
4942  const int rI = _r[countI];
4943  const int rJ = vec._r[countJ];
4944  if (rI > rJ) {
4945  ++countJ;
4946  } else if (rJ > rI) {
4947  ++countI;
4948  } else {
4949  sum+=_v[countI]*vec._v[countJ];
4950  ++countI;
4951  ++countJ;
4952  }
4953  }
4954  return sum;
4955 };
4956 
4958 template <typename T> inline void SpVector<T>::clear() {
4959  if (!_externAlloc) {
4960  delete[](_v);
4961  delete[](_r);
4962  }
4963  _v=NULL;
4964  _r=NULL;
4965  _L=0;
4966  _nzmax=0;
4967  _externAlloc=true;
4968 };
4969 
4971 template <typename T> inline void SpVector<T>::resize(const int nzmax) {
4972  if (_nzmax != nzmax) {
4973  clear();
4974  _nzmax=nzmax;
4975  _L=0;
4976  _externAlloc=false;
4977 #pragma omp critical
4978  {
4979  _v=new T[nzmax];
4980  _r=new int[nzmax];
4981  }
4982  }
4983 };
4984 
4985 template <typename T> void inline SpVector<T>::toSpMatrix(
4986  SpMatrix<T>& out, const int m, const int n) const {
4987  out.resize(m,n,_L);
4988  cblas_copy<T>(_L,_v,1,out._v,1);
4989  int current_col=0;
4990  T* out_v=out._v;
4991  int* out_r=out._r;
4992  int* out_pB=out._pB;
4993  out_pB[0]=current_col;
4994  for (int i = 0; i<_L; ++i) {
4995  int col=_r[i]/m;
4996  if (col > current_col) {
4997  out_pB[current_col+1]=i;
4998  current_col++;
4999  i--;
5000  } else {
5001  out_r[i]=_r[i]-col*m;
5002  }
5003  }
5004  for (current_col++ ; current_col < n+1; ++current_col)
5005  out_pB[current_col]=_L;
5006 };
5007 
5008 template <typename T> void inline SpVector<T>::toFull(Vector<T>& out)
5009  const {
5010  out.setZeros();
5011  T* X = out.rawX();
5012  for (int i = 0; i<_L; ++i)
5013  X[_r[i]]=_v[i];
5014  };
5015 
5016 /* ****************************
5017  * Implementaton of ProdMatrix
5018  * ****************************/
5019 
5020 template <typename T> ProdMatrix<T>::ProdMatrix() {
5021  _DtX= NULL;
5022  _X=NULL;
5023  _D=NULL;
5024  _high_memory=true;
5025  _n=0;
5026  _m=0;
5027  _addDiag=0;
5028 };
5029 
5031 template <typename T> ProdMatrix<T>::ProdMatrix(const Matrix<T>& D,
5032  const bool high_memory) {
5033  if (high_memory) _DtX = new Matrix<T>();
5034  this->setMatrices(D,high_memory);
5035 };
5036 
5038 template <typename T> ProdMatrix<T>::ProdMatrix(const Matrix<T>& D, const Matrix<T>& X,
5039  const bool high_memory) {
5040  if (high_memory) _DtX = new Matrix<T>();
5041  this->setMatrices(D,X,high_memory);
5042 };
5043 
5044 template <typename T> inline void ProdMatrix<T>::setMatrices(const Matrix<T>& D, const Matrix<T>& X,
5045  const bool high_memory) {
5046  _high_memory=high_memory;
5047  _m = D.n();
5048  _n = X.n();
5049  if (high_memory) {
5050  D.mult(X,*_DtX,true,false);
5051  } else {
5052  _X=&X;
5053  _D=&D;
5054  _DtX=NULL;
5055  }
5056  _addDiag=0;
5057 };
5058 
5059 template <typename T> inline void ProdMatrix<T>::setMatrices(
5060  const Matrix<T>& D, const bool high_memory) {
5061  _high_memory=high_memory;
5062  _m = D.n();
5063  _n = D.n();
5064  if (high_memory) {
5065  D.XtX(*_DtX);
5066  } else {
5067  _X=&D;
5068  _D=&D;
5069  _DtX=NULL;
5070  }
5071  _addDiag=0;
5072 };
5073 
5075 template <typename T> inline void ProdMatrix<T>::copyCol(const int i, Vector<T>& DtXi) const {
5076  if (_high_memory) {
5077  _DtX->copyCol(i,DtXi);
5078  } else {
5079  Vector<T> Xi;
5080  _X->refCol(i,Xi);
5081  _D->multTrans(Xi,DtXi);
5082  if (_addDiag && _m == _n) DtXi[i] += _addDiag;
5083  }
5084 };
5085 
5087 template <typename T> inline void ProdMatrix<T>::extract_rawCol(const int i,T* DtXi) const {
5088  if (_high_memory) {
5089  _DtX->extract_rawCol(i,DtXi);
5090  } else {
5091  Vector<T> Xi;
5092  Vector<T> vDtXi(DtXi,_m);
5093  _X->refCol(i,Xi);
5094  _D->multTrans(Xi,vDtXi);
5095  if (_addDiag && _m == _n) DtXi[i] += _addDiag;
5096  }
5097 };
5098 
5099 template <typename T> inline void ProdMatrix<T>::add_rawCol(const int i,T* DtXi,
5100  const T a) const {
5101  if (_high_memory) {
5102  _DtX->add_rawCol(i,DtXi,a);
5103  } else {
5104  Vector<T> Xi;
5105  Vector<T> vDtXi(DtXi,_m);
5106  _X->refCol(i,Xi);
5107  _D->multTrans(Xi,vDtXi,a,T(1.0));
5108  if (_addDiag && _m == _n) DtXi[i] += a*_addDiag;
5109  }
5110 };
5111 
5112 template <typename T> void inline ProdMatrix<T>::addDiag(const T diag) {
5113  if (_m == _n) {
5114  if (_high_memory) {
5115  _DtX->addDiag(diag);
5116  } else {
5117  _addDiag=diag;
5118  }
5119  }
5120 };
5121 
5122 template <typename T> inline T ProdMatrix<T>::operator[](const int index) const {
5123  if (_high_memory) {
5124  return (*_DtX)[index];
5125  } else {
5126  const int index2=index/this->_m;
5127  const int index1=index-this->_m*index2;
5128  Vector<T> col1, col2;
5129  _D->refCol(index1,col1);
5130  _X->refCol(index2,col2);
5131  return col1.dot(col2);
5132  }
5133 };
5134 
5135 
5136 template <typename T> inline T ProdMatrix<T>::operator()(const int index1,
5137  const int index2) const {
5138  if (_high_memory) {
5139  return (*_DtX)(index1,index2);
5140  } else {
5141  Vector<T> col1, col2;
5142  _D->refCol(index1,col1);
5143  _X->refCol(index2,col2);
5144  return col1.dot(col2);
5145  }
5146 };
5147 
5148 template <typename T> void inline ProdMatrix<T>::diag(Vector<T>& diag) const {
5149  if (_m == _n) {
5150  if (_high_memory) {
5151  _DtX->diag(diag);
5152  } else {
5153  Vector<T> col1, col2;
5154  for (int i = 0; i <_m; ++i) {
5155  _D->refCol(i,col1);
5156  _X->refCol(i,col2);
5157  diag[i] = col1.dot(col2);
5158  }
5159  }
5160  }
5161 };
5162 
5163 template <typename T> class SubMatrix : public AbstractMatrix<T> {
5164 
5165  public:
5166  SubMatrix(AbstractMatrix<T>& G, Vector<int>& indI, Vector<int>& indJ);
5167 
5168  void inline convertIndicesI(Vector<int>& ind) const;
5169  void inline convertIndicesJ(Vector<int>& ind) const;
5170  int inline n() const { return _indicesJ.n(); };
5171  int inline m() const { return _indicesI.n(); };
5172  void inline extract_rawCol(const int i, T* pr) const;
5174  inline void copyCol(const int i, Vector<T>& DtXi) const;
5176  inline void add_rawCol(const int i, T* DtXi, const T a) const;
5178  inline void diag(Vector<T>& diag) const;
5179  inline T operator()(const int index1, const int index2) const;
5180 
5181  private:
5182  Vector<int> _indicesI;
5183  Vector<int> _indicesJ;
5184  AbstractMatrix<T>* _matrix;
5185 };
5186 
5187 template <typename T>
5188 SubMatrix<T>::SubMatrix(AbstractMatrix<T>& G, Vector<int>& indI, Vector<int>& indJ) {
5189  _matrix = &G;
5190  _indicesI.copy(indI);
5191  _indicesJ.copy(indJ);
5192 };
5193 
5194 template <typename T> void inline SubMatrix<T>::convertIndicesI(
5195  Vector<int>& ind) const {
5196  int* pr_ind = ind.rawX();
5197  for (int i = 0; i<ind.n(); ++i) {
5198  if (pr_ind[i] == -1) break;
5199  pr_ind[i]=_indicesI[pr_ind[i]];
5200  }
5201 };
5202 
5203 template <typename T> void inline SubMatrix<T>::convertIndicesJ(
5204  Vector<int>& ind) const {
5205  int* pr_ind = ind.rawX();
5206  for (int i = 0; i<ind.n(); ++i) {
5207  if (pr_ind[i] == -1) break;
5208  pr_ind[i]=_indicesJ[pr_ind[i]];
5209  }
5210 };
5211 
5212 template <typename T> void inline SubMatrix<T>::extract_rawCol(const int i, T* pr) const {
5213  int* pr_ind=_indicesI.rawX();
5214  int* pr_ind2=_indicesJ.rawX();
5215  for (int j = 0; j<_indicesI.n(); ++j) {
5216  pr[j]=(*_matrix)(pr_ind[j],pr_ind2[i]);
5217  }
5218 };
5219 
5220 template <typename T> inline void SubMatrix<T>::copyCol(const int i,
5221  Vector<T>& DtXi) const {
5222  this->extract_rawCol(i,DtXi.rawX());
5223 };
5224 
5225 template <typename T> void inline SubMatrix<T>::add_rawCol(const int i, T* pr,
5226  const T a) const {
5227  int* pr_ind=_indicesI.rawX();
5228  int* pr_ind2=_indicesJ.rawX();
5229  for (int j = 0; j<_indicesI.n(); ++j) {
5230  pr[j]+=a*(*_matrix)(pr_ind[j],pr_ind2[i]);
5231  }
5232 };
5233 
5234 template <typename T> void inline SubMatrix<T>::diag(Vector<T>& diag) const {
5235  T* pr = diag.rawX();
5236  int* pr_ind=_indicesI.rawX();
5237  for (int j = 0; j<_indicesI.n(); ++j) {
5238  pr[j]=(*_matrix)(pr_ind[j],pr_ind[j]);
5239  }
5240 };
5241 
5242 template <typename T> inline T SubMatrix<T>::operator()(const int index1,
5243  const int index2) const {
5244  return (*_matrix)(_indicesI[index1],_indicesJ[index2]);
5245 }
5246 
5248 template <typename T> class ShiftMatrix : public AbstractMatrixB<T> {
5249  public:
5250  ShiftMatrix(const AbstractMatrixB<T>& inputmatrix, const int shifts, const bool center = false) : _shifts(shifts), _inputmatrix(&inputmatrix), _centered(false) {
5251  _m=_inputmatrix->m()-shifts+1;
5252  _n=_inputmatrix->n()*shifts;
5253  if (center) this->center();
5254  };
5255  int n() const { return _n; };
5256  int m() const { return _m; };
5257 
5259  void multTrans(const Vector<T>& x, Vector<T>& b,
5260  const T alpha = 1.0, const T beta = 0.0) const;
5261 
5263  virtual void mult(const SpVector<T>& x, Vector<T>& b,
5264  const T alpha = 1.0, const T beta = 0.0) const;
5265 
5266  virtual void mult(const Vector<T>& x, Vector<T>& b,
5267  const T alpha = 1.0, const T beta = 0.0) const;
5268 
5270  virtual void mult(const Matrix<T>& B, Matrix<T>& C,
5271  const bool transA = false, const bool transB = false,
5272  const T a = 1.0, const T b = 0.0) const;
5273 
5274  virtual void mult(const SpMatrix<T>& B, Matrix<T>& C,
5275  const bool transA = false, const bool transB = false,
5276  const T a = 1.0, const T b = 0.0) const;
5277 
5279  virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C,
5280  const bool transA = false, const bool transB = false,
5281  const T a = 1.0, const T b = 0.0) const;
5282 
5284  virtual void XtX(Matrix<T>& XtX) const;
5285 
5286  virtual void copyRow(const int i, Vector<T>& x) const;
5287 
5288  virtual void copyTo(Matrix<T>& copy) const;
5289  virtual T dot(const Matrix<T>& x) const;
5290 
5291  virtual void print(const string& name) const;
5292 
5293  virtual ~ShiftMatrix() { };
5294 
5295  private:
5296  void center() {
5297  Vector<T> ones(_m);
5298  ones.set(T(1.0)/_m);
5299  this->multTrans(ones,_means);
5300  _centered=true; };
5301 
5302  int _m;
5303  int _n;
5304  int _shifts;
5305  bool _centered;
5306  Vector<T> _means;
5307  const AbstractMatrixB<T>* _inputmatrix;
5308 };
5309 
5310 template <typename T> void ShiftMatrix<T>::multTrans(const
5311  Vector<T>& x, Vector<T>& b, const T alpha, const T beta) const {
5312  b.resize(_n);
5313  if (beta==0) b.setZeros();
5314  Vector<T> tmp(_inputmatrix->m());
5315  Vector<T> subvec;
5316  Vector<T> subvec2;
5317  const int nn=_inputmatrix->n();
5318  for (int i = 0; i<_shifts; ++i) {
5319  tmp.setZeros();
5320  subvec2.setData(tmp.rawX()+i,_m);
5321  subvec2.copy(x);
5322  subvec.setData(b.rawX()+i*nn,nn);
5323  _inputmatrix->multTrans(tmp,subvec,alpha,beta);
5324  }
5325  if (_centered) {
5326  b.add(_means,-alpha*x.sum());
5327  }
5328 };
5329 
5330 
5332 template <typename T> void ShiftMatrix<T>::mult(const
5333  SpVector<T>& x, Vector<T>& b, const T alpha, const T beta) const {
5334  b.resize(_m);
5335  if (beta==0) {
5336  b.setZeros();
5337  } else {
5338  b.scal(beta);
5339  }
5340  const int nn=_inputmatrix->n();
5341  const int mm=_inputmatrix->m();
5342  Vector<T> fullx(_n);
5343  x.toFull(fullx);
5344  SpVector<T> sptmp(nn);
5345  Vector<T> tmp;
5346  Vector<T> tmp2(mm);
5347  for (int i = 0; i<_shifts; ++i) {
5348  tmp.setData(fullx.rawX()+i*nn,nn);
5349  tmp.toSparse(sptmp);
5350  _inputmatrix->mult(sptmp,tmp2,alpha,0);
5351  tmp.setData(tmp2.rawX()+i,_m);
5352  b.add(tmp);
5353  }
5354  if (_centered) {
5355  b.add(-alpha*_means.dot(x));
5356  }
5357 };
5358 
5360 template <typename T> void ShiftMatrix<T>::mult(const
5361  Vector<T>& x, Vector<T>& b, const T alpha, const T beta) const {
5362  b.resize(_m);
5363  const int nn=_inputmatrix->n();
5364  const int mm=_inputmatrix->m();
5365  Vector<T> tmp;
5366  Vector<T> tmp2(mm);
5367  if (beta==0) {
5368  b.setZeros();
5369  } else {
5370  b.scal(beta);
5371  }
5372  for (int i = 0; i<_shifts; ++i) {
5373  tmp.setData(x.rawX()+i*nn,nn);
5374  _inputmatrix->mult(tmp,tmp2,alpha,0);
5375  tmp.setData(tmp2.rawX()+i,_m);
5376  b.add(tmp);
5377  }
5378  if (_centered) {
5379  b.add(-alpha*_means.dot(x));
5380  }
5381 };
5382 
5384 template <typename T> void ShiftMatrix<T>::mult(const Matrix<T>&
5385  B, Matrix<T>& C, const bool transA, const bool transB, const T a, const T
5386  b) const {
5387  cerr << "Shift Matrix is used in inadequate setting" << endl;
5388 }
5389 
5390 template <typename T> void ShiftMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
5391  const bool transA, const bool transB, const T a, const T b) const {
5392  cerr << "Shift Matrix is used in inadequate setting" << endl;
5393 }
5394 
5396 template <typename T> void ShiftMatrix<T>::multSwitch(const
5397  Matrix<T>& B, Matrix<T>& C, const bool transA, const bool transB,
5398  const T a, const T b) const {
5399  cerr << "Shift Matrix is used in inadequate setting" << endl;
5400 }
5401 
5402 template <typename T> void ShiftMatrix<T>::XtX(Matrix<T>& XtX) const {
5403  cerr << "Shift Matrix is used in inadequate setting" << endl;
5404 };
5405 
5406 template <typename T> void ShiftMatrix<T>::copyRow(const int ind, Vector<T>& x) const {
5407  Vector<T> sub_vec;
5408  const int mm=_inputmatrix->m();
5409  for (int i = 0; i<_shifts; ++i) {
5410  sub_vec.setData(x.rawX()+i*mm,mm);
5411  _inputmatrix->copyRow(ind+i,sub_vec);
5412  }
5413  if (_centered) x.sub(_means);
5414 };
5415 
5416 template <typename T> void ShiftMatrix<T>::copyTo(Matrix<T>& x) const {
5417  cerr << "Shift Matrix is used in inadequate setting" << endl;
5418 };
5419 
5420 
5421 template <typename T> T ShiftMatrix<T>::dot(const Matrix<T>& x) const {
5422  cerr << "Shift Matrix is used in inadequate setting" << endl;
5423  return 0;
5424 };
5425 
5426 template <typename T> void ShiftMatrix<T>::print(const string& name) const {
5427  cerr << name << endl;
5428  cerr << "Shift Matrix: " << _shifts << " shifts" << endl;
5429  _inputmatrix->print(name);
5430 };
5431 
5433 template <typename T> class DoubleRowMatrix : public AbstractMatrixB<T> {
5434  public:
5435  DoubleRowMatrix(const AbstractMatrixB<T>& inputmatrix) : _inputmatrix(&inputmatrix) {
5436  _n=inputmatrix.n();
5437  _m=2*inputmatrix.m();
5438  };
5439  int n() const { return _n; };
5440  int m() const { return _m; };
5441 
5443  void multTrans(const Vector<T>& x, Vector<T>& b,
5444  const T alpha = 1.0, const T beta = 0.0) const;
5445 
5447  virtual void mult(const SpVector<T>& x, Vector<T>& b,
5448  const T alpha = 1.0, const T beta = 0.0) const;
5449 
5450  virtual void mult(const Vector<T>& x, Vector<T>& b,
5451  const T alpha = 1.0, const T beta = 0.0) const;
5452 
5454  virtual void mult(const Matrix<T>& B, Matrix<T>& C,
5455  const bool transA = false, const bool transB = false,
5456  const T a = 1.0, const T b = 0.0) const;
5457 
5458  virtual void mult(const SpMatrix<T>& B, Matrix<T>& C,
5459  const bool transA = false, const bool transB = false,
5460  const T a = 1.0, const T b = 0.0) const;
5461 
5463  virtual void multSwitch(const Matrix<T>& B, Matrix<T>& C,
5464  const bool transA = false, const bool transB = false,
5465  const T a = 1.0, const T b = 0.0) const;
5466 
5468  virtual void XtX(Matrix<T>& XtX) const;
5469 
5470  virtual void copyRow(const int i, Vector<T>& x) const;
5471 
5472  virtual void copyTo(Matrix<T>& copy) const;
5473  virtual T dot(const Matrix<T>& x) const;
5474 
5475  virtual void print(const string& name) const;
5476 
5477  virtual ~DoubleRowMatrix() { };
5478 
5479  private:
5480  int _m;
5481  int _n;
5482  const AbstractMatrixB<T>* _inputmatrix;
5483 };
5484 
5485 
5486 template <typename T> void DoubleRowMatrix<T>::multTrans(const
5487  Vector<T>& x, Vector<T>& b, const T alpha, const T beta) const {
5488  const int mm = _inputmatrix->m();
5489  Vector<T> tmp(mm);
5490  for (int i = 0; i<mm; ++i)
5491  tmp[i]=x[2*i]+x[2*i+1];
5492  _inputmatrix->multTrans(tmp,b,alpha,beta);
5493 };
5494 
5495 
5497 template <typename T> void DoubleRowMatrix<T>::mult(const
5498  SpVector<T>& x, Vector<T>& b, const T alpha, const T beta) const {
5499  b.resize(_m);
5500  if (beta==0) {
5501  b.setZeros();
5502  } else {
5503  b.scal(beta);
5504  }
5505  const int mm = _inputmatrix->m();
5506  Vector<T> tmp(mm);
5507  _inputmatrix->mult(x,tmp,alpha);
5508  for (int i = 0; i<mm; ++i) {
5509  b[2*i]+=tmp[i];
5510  b[2*i+1]+=tmp[i];
5511  }
5512 };
5513 
5515 template <typename T> void DoubleRowMatrix<T>::mult(const
5516  Vector<T>& x, Vector<T>& b, const T alpha, const T beta) const {
5517  b.resize(_m);
5518  if (beta==0) {
5519  b.setZeros();
5520  } else {
5521  b.scal(beta);
5522  }
5523  const int mm = _inputmatrix->m();
5524  Vector<T> tmp(mm);
5525  _inputmatrix->mult(x,tmp,alpha);
5526  for (int i = 0; i<mm; ++i) {
5527  b[2*i]+=tmp[i];
5528  b[2*i+1]+=tmp[i];
5529  }
5530 };
5531 
5533 template <typename T> void DoubleRowMatrix<T>::mult(const Matrix<T>&
5534  B, Matrix<T>& C, const bool transA, const bool transB, const T a, const T
5535  b) const {
5536  FLAG(5)
5537  cerr << "Double Matrix is used in inadequate setting" << endl;
5538 }
5539 
5540 template <typename T> void DoubleRowMatrix<T>::mult(const SpMatrix<T>& B, Matrix<T>& C,
5541  const bool transA, const bool transB, const T a, const T b) const {
5542  FLAG(4)
5543  cerr << "Double Matrix is used in inadequate setting" << endl;
5544 }
5545 
5547 template <typename T> void DoubleRowMatrix<T>::multSwitch(const
5548  Matrix<T>& B, Matrix<T>& C, const bool transA, const bool transB,
5549  const T a, const T b) const {
5550  FLAG(3)
5551  cerr << "Double Matrix is used in inadequate setting" << endl;
5552 }
5553 
5554 template <typename T> void DoubleRowMatrix<T>::XtX(Matrix<T>& XtX) const {
5555  FLAG(2)
5556  cerr << "Double Matrix is used in inadequate setting" << endl;
5557 };
5558 
5559 template <typename T> void DoubleRowMatrix<T>::copyRow(const int ind, Vector<T>& x) const {
5560  const int indd2= static_cast<int>(floor(static_cast<double>(ind)/2.0));
5561  _inputmatrix->copyRow(indd2,x);
5562 };
5563 
5564 template <typename T> void DoubleRowMatrix<T>::copyTo(Matrix<T>& x) const {
5565  FLAG(1)
5566  cerr << "Double Matrix is used in inadequate setting" << endl;
5567 };
5568 
5569 
5570 template <typename T> T DoubleRowMatrix<T>::dot(const Matrix<T>& x) const {
5571  FLAG(0)
5572  cerr << "Double Matrix is used in inadequate setting" << endl;
5573  return 0;
5574 };
5575 
5576 template <typename T> void DoubleRowMatrix<T>::print(const string& name) const {
5577  cerr << name << endl;
5578  cerr << "Double Row Matrix" << endl;
5579  _inputmatrix->print(name);
5580 };
5581 
5582 }
5583 
5584 #endif
spams::Matrix::normalize
void normalize()
Normalize all columns to unit l2 norm.
Definition: linalg.h:1270
spams::ProdMatrix::add_rawCol
virtual void add_rawCol(const int i, T *DtXi, const T a) const
compute DtX(:,i)
Definition: linalg.h:5099
spams::SpMatrix::add_direct
void add_direct(const SpMatrix< T > &mat, const T a)
Definition: linalg.h:4183
spams::Matrix::norm_l1_rows
void norm_l1_rows(Vector< T > &norms) const
returns the linf norms of the columns
Definition: linalg.h:2193
spams::group
std::list< int > group
Definition: linalg.h:67
spams::Matrix::center
void center()
center the columns of the matrix
Definition: linalg.h:1299
spams::Matrix::sub
void sub(const Matrix< T > &mat)
substract the matrix mat to the current matrix
Definition: linalg.h:2080
spams::Matrix::addVecToCols
void addVecToCols(const Vector< T > &diag, const T a=1.0)
Definition: linalg.h:1526
spams::Matrix::incrDiag
void incrDiag()
add one to the diagonal
Definition: linalg.h:1502
spams::SpVector::nzmax
T nzmax() const
Definition: linalg.h:935
spams::Matrix::print
void print(const string &name) const
Print the matrix to std::cout.
Definition: linalg.h:1085
spams::SpMatrix::scal
void scal(const T a) const
scale the matrix by a
Definition: linalg.h:4240
spams::ShiftMatrix::multSwitch
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
perform C = a*B*A + b*C, possibly transposing A or B.
Definition: linalg.h:5396
spams::SpMatrix::r
int * r() const
Direct access to _r.
Definition: linalg.h:808
spams::SpMatrix::_n
int _n
number of columns
Definition: linalg.h:912
spams::Matrix::addToCols
void addToCols(const Vector< T > &diag)
Definition: linalg.h:1517
spams::SpMatrix::addVecToColsWeighted
void addVecToColsWeighted(const Vector< T > &diag, const T *weights, const T a=1.0)
Definition: linalg.h:4509
spams::Matrix::quad
T quad(const Vector< T > &vec1, const SpVector< T > &vec2) const
return vec1&#39;*A*vec2, where vec2 is sparse
Definition: linalg.h:2047
spams::Vector::normalize2
void normalize2()
normalize the vector
Definition: linalg.h:3128
spams::SpMatrix
Sparse Matrix class.
Definition: linalg.h:63
spams::Vector::nnz
int nnz() const
Definition: linalg.h:2822
spams::SpVector::asum
T asum() const
computes the sum of the magnitude of the elements
Definition: linalg.h:4891
spams::SpMatrix::toFull
void toFull(Matrix< T > &matrix) const
copy the sparse matrix into a dense matrix
Definition: linalg.h:4760
spams::SpMatrix::copyRow
void copyRow(const int i, Vector< T > &x) const
Definition: linalg.h:4479
spams::Matrix::multSwitch
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.
Definition: linalg.h:1840
spams::ShiftMatrix::multTrans
void multTrans(const Vector< T > &x, Vector< T > &b, const T alpha=1.0, const T beta=0.0) const
b <- alpha A&#39;x + beta b
Definition: linalg.h:5310
upper
static char upper
Definition: cblas_alt_template.h:23
spams::Vector::neg
void neg()
A <- -A.
Definition: linalg.h:3282
spams
Definition: dag.h:26
spams::Matrix::rank1Update
void rank1Update(const Vector< T > &vec1, const Vector< T > &vec2, const T alpha=1.0)
perform A <- A + alpha*vec1*vec2&#39;
Definition: linalg.h:2384
spams::ShiftMatrix::XtX
virtual void XtX(Matrix< T > &XtX) const
XtX = A&#39;*A.
Definition: linalg.h:5402
spams::SpVector::setL
void setL(const int L)
Definition: linalg.h:960
spams::ProdMatrix::_X
const Matrix< T > * _X
Definition: linalg.h:1038
spams::Matrix::inv_elem
void inv_elem()
inverse the elements of the matrix
Definition: linalg.h:2313
spams::Vector::sum
T sum() const
returns the sum of the vector
Definition: linalg.h:3344
spams::SpMatrix::wAAt
void wAAt(const Vector< T > &w, Matrix< T > &aat) const
aat <- sum_i w_i A(:,i)*A(:,i)&#39;
Definition: linalg.h:4619
spams::Matrix::norm_2_cols
void norm_2_cols(Vector< T > &norms) const
returns the l2 norms of the columns
Definition: linalg.h:2162
spams::Data::m
virtual int m() const =0
spams::Matrix::unwhiten
void unwhiten(Vector< T > &mean, const bool pattern=false)
whiten
Definition: linalg.h:1460
spams::Vector::Invsqrt
void Invsqrt()
A <- 1 ./ sqrt(A)
Definition: linalg.h:3093
CBLAS_TRANSPOSE
CBLAS_TRANSPOSE
Definition: utl_cblas.h:6
spams::SpVector::rawX
T * rawX() const
Definition: linalg.h:956
spams::Matrix::meanCol
void meanCol(Vector< T > &mean) const
Compute the mean of the columns.
Definition: linalg.h:2226
spams::Matrix::fmaxval
T fmaxval() const
return the 1D-index of the value of greatest magnitude
Definition: linalg.h:1193
spams::Vector::sign
void sign(Vector< T > &signs) const
puts in signs, the sign of each point in the vector
Definition: linalg.h:3351
spams::Matrix::_externAlloc
bool _externAlloc
is the data allocation external or not
Definition: linalg.h:538
misc.h
Contains miscellaneous functions.
spams::Vector::min
int min() const
returns the index of the minimum value
Definition: linalg.h:2775
spams::Matrix::sum_cols
void sum_cols(Vector< T > &sum) const
whiten
Definition: linalg.h:2215
spams::alt_log
static T alt_log(const T x)
spams::SpMatrix::pB
int pB(const int i) const
returns pB[i]
Definition: linalg.h:779
spams::Vector::Sqrt
void Sqrt(const Vector< T > &x)
A <- 1 ./ sqrt(x)
Definition: linalg.h:3088
spams::Vector::copy
void copy(const Vector< T > &x)
make a copy of x
Definition: linalg.h:2865
spams::Vector::normalize
void normalize()
normalize the vector
Definition: linalg.h:3122
cblas_alt_template.h
spams::Vector::sparseProject
void sparseProject(Vector< 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)
Definition: linalg.h:3508
allV
static char allV
Definition: cblas_alt_template.h:30
spams::Matrix::norm_2sq_rows
void norm_2sq_rows(Vector< T > &norms) const
returns the l2 norms of the columns
Definition: linalg.h:2151
spams::Matrix::add_rawCol
virtual void add_rawCol(const int i, T *DtXi, const T a) const
Copy the column i into x.
Definition: linalg.h:1121
spams::Vector::scal
void scal(const T a)
scale the vector by a
Definition: linalg.h:3277
spams::Vector::nrm2
T nrm2() const
returns ||A||_2
Definition: linalg.h:3002
spams::Vector::inv
void inv(const Vector< T > &x)
A <- 1./x.
Definition: linalg.h:3103
spams::SpMatrix::XtX
void XtX(Matrix< T > &XtX) const
XtX = A&#39;*A.
Definition: linalg.h:4569
spams::sqr
T sqr(const T x)
template version of the fabs function
spams::SpMatrix::normFsq
T normFsq() const
compute the sum of the matrix elements
Definition: linalg.h:4178
spams::Matrix::setAleat
void setAleat()
Put white Gaussian noise in the matrix.
Definition: linalg.h:1259
spams::Matrix::setn
void setn(const int n)
modify _n
Definition: linalg.h:279
spams::Data::~Data
virtual ~Data()
Definition: linalg.h:140
spams::Vector::log
void log()
replace each value by its logarithm
Definition: linalg.h:3292
spams::Vector::l1l2projectb
void l1l2projectb(Vector< T > &out, const T thrs, const T gamma, const bool pos=false, const int mode=1)
returns true if the returned vector is null
Definition: linalg.h:3567
spams::Vector::Vector
Vector()
Empty constructor.
Definition: linalg.h:2693
spams::Data::V
virtual int V() const =0
spams::Vector::fmaxval
T fmaxval() const
returns the maximum magnitude
Definition: linalg.h:2799
spams::Matrix::copyRow
void copyRow(const int i, Vector< T > &x) const
Copy the column i into x.
Definition: linalg.h:1107
CblasTrans
Definition: utl_cblas.h:6
spams::AbstractMatrixB::n
virtual int n() const =0
spams::SubMatrix::add_rawCol
void add_rawCol(const int i, T *DtXi, const T a) const
compute DtX(:,i)
Definition: linalg.h:5225
spams::SubMatrix::_indicesJ
Vector< int > _indicesJ
Definition: linalg.h:5183
spams::Vector::fmin
int fmin() const
returns the index of the value with smallest magnitude
Definition: linalg.h:2848
spams::SpMatrix::v
T * v() const
Direct access to _v.
Definition: linalg.h:810
spams::Vector::add
void add(const Vector< T > &x, const T a=1.0)
A <- A + a*x.
Definition: linalg.h:3029
spams::SpMatrix::_pB
int * _pB
indices of the beginning of columns
Definition: linalg.h:906
spams::Matrix::XtX
void XtX(Matrix< T > &XtX) const
XtX = A&#39;*A.
Definition: linalg.h:1959
spams::SpMatrix::dot_direct
T dot_direct(const SpMatrix< T > &mat) const
Definition: linalg.h:4197
spams::Vector::_X
T * _X
data
Definition: linalg.h:754
spams::Matrix::addDiag
void addDiag(const Vector< T > &diag)
Definition: linalg.h:1506
spams::Vector::KL
T KL(const Vector< T > &X)
compute the Kuhlback-Leiber divergence
Definition: linalg.h:3212
spams::Vector::applyBayerPattern
void applyBayerPattern(const int offset)
sort the vector
Definition: linalg.h:3917
utl::E
static constexpr double E
Definition: utlConstants.h:25
spams::xlogx
static T xlogx(const T x)
Definition: linalg.h:108
spams::Matrix::normalize2
void normalize2()
Normalize all columns which l2 norm is greater than one.
Definition: linalg.h:1288
spams::Matrix::merge
void merge(const Matrix< T > &B, Matrix< T > &C) const
merge two dictionaries
Definition: linalg.h:2656
reduced
static char reduced
Definition: cblas_alt_template.h:29
spams::SpMatrix::pE
int * pE() const
Direct access to _pE.
Definition: linalg.h:806
spams::Matrix::multDiagLeft
void multDiagLeft(const Vector< T > &diag)
mult by a diagonal matrix on the left
Definition: linalg.h:1917
spams::Matrix::getData
void getData(Vector< T > &data, const int i) const
Copy the column i into x.
Definition: linalg.h:1127
spams::Vector::inv
void inv()
A <- 1./A.
Definition: linalg.h:3109
spams::Vector::mult
void mult(const Vector< T > &x, const Vector< T > &y)
A <- x .* y.
Definition: linalg.h:3114
spams::Vector::sort2
void sort2(Vector< T > &out, Vector< int > &key, const bool mode) const
Definition: linalg.h:3911
spams::Matrix::setDiag
void setDiag(const Vector< T > &d)
set the diagonal
Definition: linalg.h:1994
spams::Vector::rawX
T * rawX() const
returns a modifiable reference of the data, DANGEROUS
Definition: linalg.h:593
spams::Matrix::operator[]
T operator[](const int index) const
Return the value X(i) (1D indexing)
Definition: linalg.h:232
CblasUpper
Definition: utl_cblas.h:8
spams::Vector::std
T std()
whiten
Definition: linalg.h:3265
spams::vector_groups
std::vector< group > vector_groups
Definition: linalg.h:73
spams::Matrix::transpose
void transpose(Matrix< T > &trans)
Definition: linalg.h:1489
spams::Matrix::norm_inf_rows
void norm_inf_rows(Vector< T > &norms) const
returns the linf norms of the columns
Definition: linalg.h:2184
spams::Matrix::Matrix
Matrix()
Empty constructor.
Definition: linalg.h:1066
spams::ProdMatrix::_DtX
Matrix< T > * _DtX
Depending on the mode, DtX is a matrix, or two matrices.
Definition: linalg.h:1037
spams::Matrix::subMatrixSym
void subMatrixSym(const Vector< int > &indices, Matrix< T > &subMatrix) const
extract a sub-matrix of a symmetric matrix
Definition: linalg.h:1204
spams::SpVector::L
int L() const
Definition: linalg.h:958
spams::SubMatrix::convertIndicesJ
void convertIndicesJ(Vector< int > &ind) const
Definition: linalg.h:5203
spams::Vector::mult_elementWise
void mult_elementWise(const Vector< T > &B, Vector< T > &C) const
Definition: linalg.h:670
spams::SpMatrix::_nzmax
int _nzmax
number of non-zero values
Definition: linalg.h:914
std
STL namespace.
lower
static char lower
Definition: cblas_alt_template.h:21
spams::SpMatrix::copy_direct
void copy_direct(const SpMatrix< T > &mat)
Definition: linalg.h:4190
spams::sort
static void sort(int *irOut, T *prOut, int beg, int end)
Definition: misc.h:141
spams::Matrix::copyTo
void copyTo(Matrix< T > &mat) const
make a copy of the matrix mat in the current matrix
Definition: linalg.h:260
spams::SpMatrix::wXAt
void wXAt(const Vector< T > &w, const Matrix< T > &X, Matrix< T > &XAt, const int numthreads=-1) const
XAt <- sum_i w_i X(:,i)*A(:,i)&#39;.
Definition: linalg.h:4721
spams::Vector::softmax
T softmax(const int y)
replace each value by its exponential
Definition: linalg.h:3308
spams::Vector::maxval
T maxval() const
returns the maximum value
Definition: linalg.h:2789
spams::Vector::sqr
void sqr(const Vector< T > &x)
A <- x .^ 2.
Definition: linalg.h:3077
spams::softThrs
static T softThrs(const T x, const T lambda)
Definition: linalg.h:87
spams::SpMatrix::nzmax
int nzmax() const
returns the maximum number of non-zero elements
Definition: linalg.h:785
spams::Matrix::diag
void diag(Vector< T > &d) const
extract the diagonal
Definition: linalg.h:1985
spams::SpMatrix::operator[]
T operator[](const int index) const
returns X[index]
Definition: linalg.h:4138
spams::Matrix::fillSymmetric2
void fillSymmetric2()
Definition: linalg.h:1358
spams::ProdMatrix::copyCol
void copyCol(const int i, Vector< T > &DtXi) const
compute DtX(:,i)
Definition: linalg.h:5075
spams::SpMatrix::norm_0_cols
void norm_0_cols(Vector< T > &norms) const
returns the l0 norms of the columns
Definition: linalg.h:4842
spams::Vector::setn
void setn(const int n)
Definition: linalg.h:627
spams::Vector::alltrue
bool alltrue() const
spams::Vector::copyMask
void copyMask(Vector< T > &out, Vector< bool > &mask) const
extract the rows of a matrix corresponding to a binary mask
Definition: linalg.h:4041
spams::Matrix::conjugateGradient
void conjugateGradient(const Vector< T > &b, Vector< T > &x, const T tol=1e-4, const int=4) const
compute x, such that b = Ax,
Definition: linalg.h:2482
spams::Matrix::copyCol
void copyCol(const int i, Vector< T > &x) const
Copy the column i into x.
Definition: linalg.h:1100
spams::Vector::afused
T afused() const
compute the sum of the differences
Definition: linalg.h:3336
CblasColMajor
Definition: utl_cblas.h:5
spams::Matrix::~Matrix
virtual ~Matrix()
Destructor.
Definition: linalg.h:1070
spams::Matrix::addRow
void addRow(const int i, const Vector< T > &row, const T a=1.0)
fill the matrix with the row given
Definition: linalg.h:2268
spams::SpMatrix::refCol
void refCol(int i, SpVector< T > &vec) const
reference the column i into vec
Definition: linalg.h:4116
spams::DoubleRowMatrix::XtX
virtual void XtX(Matrix< T > &XtX) const
XtX = A&#39;*A.
Definition: linalg.h:5554
spams::SubMatrix::convertIndicesI
void convertIndicesI(Vector< int > &ind) const
Definition: linalg.h:5194
spams::ProdMatrix::operator()
T operator()(const int index1, const int index2) const
returns the value of an index
Definition: linalg.h:5136
spams::Matrix::thrsmin
void thrsmin(const T nu)
perform thresholding of the matrix, with the threshold nu
Definition: linalg.h:2305
spams::ShiftMatrix::m
int m() const
Definition: linalg.h:5256
spams::Vector::unwhiten
void unwhiten(Vector< T > &mean, const bool pattern=false)
whiten
Definition: linalg.h:3233
spams::Matrix::fillSymmetric
void fillSymmetric()
Definition: linalg.h:1351
spams::Matrix::m
int m() const
Definition: linalg.h:222
spams::Matrix::_X
T * _X
pointer to the data
Definition: linalg.h:540
spams::SubMatrix::diag
void diag(Vector< T > &diag) const
compute DtX(:,i)
Definition: linalg.h:5234
spams::SpMatrix::multSwitch
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.
Definition: linalg.h:4461
utl::max
const T & max(const T &a, const T &b)
Return the maximum between a and b.
Definition: utlCore.h:263
spams::Vector::sort
void sort(Vector< T > &out, const bool mode) const
sort the vector
Definition: linalg.h:3899
spams::SpVector::SpVector
SpVector()
Empty constructor.
Definition: linalg.h:4882
spams::Matrix::svdRankOne
void svdRankOne(const Vector< T > &u0, Vector< T > &u, Vector< T > &v) const
Definition: linalg.h:1537
spams::DoubleRowMatrix::dot
virtual T dot(const Matrix< T > &x) const
Definition: linalg.h:5570
spams::SpMatrix::copyTo
void copyTo(Matrix< T > &mat) const
make a copy of the matrix mat in the current matrix
Definition: linalg.h:866
spams::Data::operator[]
virtual T operator[](const int index) const =0
spams::Matrix::hardThrshold
void hardThrshold(const T nu)
perform soft-thresholding of the matrix, with the threshold nu
Definition: linalg.h:2290
utl::Sqrt
auto Sqrt(const ExprT &expr) -> decltype(utl::F< utl::Functor::Sqrt< typename ExprT::ValueType > >(expr))
Definition: utlFunctors.h:362
spams::SubMatrix::_matrix
AbstractMatrix< T > * _matrix
Definition: linalg.h:5184
utils.h
Contains various variables and class timer.
spams::Vector::n
int n() const
returns the size of the vector
Definition: linalg.h:591
spams::ShiftMatrix::print
virtual void print(const string &name) const
Definition: linalg.h:5426
spams::Matrix::div_elementWise
void div_elementWise(const Matrix< T > &B, Matrix< T > &C) const
C = A .* B, elementwise multiplication.
Definition: linalg.h:1950
spams::Matrix::quad_mult
void quad_mult(const Vector< T > &vec1, const SpVector< T > &vec2, Vector< T > &y, const T a=1.0, const T b=0.0) const
return vec1&#39;*A*vec2, where vec2 is sparse
Definition: linalg.h:2033
utl::Matrix
NDArray< T, 2 > Matrix
Definition: utlNDArray.h:45
spams::Matrix::neg
void neg()
A <- -A.
Definition: linalg.h:1498
spams::Data
Data class, abstract class, useful in the class image.
Definition: linalg.h:130
spams::Matrix::center_rows
void center_rows()
center the columns of the matrix
Definition: linalg.h:1309
spams::ProdMatrix::~ProdMatrix
~ProdMatrix()
Constructor, D&#39;*X is represented, with optional transpositions.
Definition: linalg.h:1010
spams::Matrix::norm_1_2_col
T norm_1_2_col() const
return ||At||_{1,2} (max of l2 norm of the columns)
Definition: linalg.h:2128
spams::ShiftMatrix::copyTo
virtual void copyTo(Matrix< T > &copy) const
Definition: linalg.h:5416
spams::SpVector::_externAlloc
bool _externAlloc
external allocation
Definition: linalg.h:985
spams::ProdMatrix::n
int n() const
returns the number of columns
Definition: linalg.h:1027
MIN
#define MIN(a, b)
Definition: utils.h:47
spams::SpMatrix::_externAlloc
bool _externAlloc
if the data has been externally allocated
Definition: linalg.h:900
spams::Vector::project_sft_binary
void project_sft_binary(const Vector< T > &labels)
Definition: linalg.h:3460
spams::Matrix::refCol
void refCol(int i, Vector< T > &x) const
Reference the column i into the vector x.
Definition: linalg.h:1144
spams::ShiftMatrix
Matrix with shifts.
Definition: linalg.h:5248
CblasNoTrans
Definition: utl_cblas.h:6
spams::Vector::_n
int _n
size of the vector
Definition: linalg.h:756
spams::Matrix::copyMask
void copyMask(Matrix< T > &out, Vector< bool > &mask) const
extract the rows of a matrix corresponding to a binary mask
Definition: linalg.h:4052
spams::ShiftMatrix::mult
virtual 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
Definition: linalg.h:5332
spams::Vector::allfalse
bool allfalse() const
spams::SpMatrix::~SpMatrix
~SpMatrix()
Destructor.
Definition: linalg.h:4111
spams::AbstractMatrixB::m
virtual int m() const =0
spams::SpVector::length
T length() const
returns the length of the vector
Definition: linalg.h:937
spams::AbstractMatrixB::~AbstractMatrixB
virtual ~AbstractMatrixB()
Definition: linalg.h:184
MAX
#define MAX(a, b)
Definition: utils.h:48
spams::SpMatrix::_pE
int * _pE
indices of the end of columns
Definition: linalg.h:908
spams::Vector::setZeros
void setZeros()
Set all values to zero.
Definition: linalg.h:2871
spams::Matrix::toSparseTrans
void toSparseTrans(SpMatrix< T > &matrixTrans)
make a sparse copy of the current matrix
Definition: linalg.h:2606
utl::min
const T & min(const T &a, const T &b)
Return the minimum between a and b.
Definition: utlCore.h:257
spams::hardThrs
static T hardThrs(const T x, const T lambda)
Definition: linalg.h:98
spams::Vector::thrshold
void thrshold(const T nu)
performs soft-thresholding of the vector
Definition: linalg.h:2967
spams::Vector::asum
T asum() const
computes the sum of the magnitudes of the vector
Definition: linalg.h:3324
spams::SpMatrix::print
void print(const string &name) const
print the sparse matrix
Definition: linalg.h:4127
spams::Vector::dot
T dot(const Vector< T > &x) const
returns A&#39;x
Definition: linalg.h:3012
spams::SpMatrix::AAt
void AAt(Matrix< T > &aat) const
aat <- A*A&#39;
Definition: linalg.h:4535
FLAG
#define FLAG(a)
Definition: utils.h:57
spams::SpMatrix::v
T v(const int i) const
returns v[i]
Definition: linalg.h:783
spams::ShiftMatrix::_means
Vector< T > _means
Definition: linalg.h:5306
spams::Vector::thrsabsmin
void thrsabsmin(const T nu)
performs thresholding of the vector
Definition: linalg.h:2960
spams::SubMatrix::operator()
T operator()(const int index1, const int index2) const
Definition: linalg.h:5242
spams::DoubleRowMatrix::mult
virtual 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
Definition: linalg.h:5497
spams::Matrix::multDiagRight
void multDiagRight(const Vector< T > &diag)
mult by a diagonal matrix on the right
Definition: linalg.h:1929
spams::Vector::whiten
void whiten(Vector< T > &mean, const bool pattern=false)
whiten
Definition: linalg.h:3134
spams::DoubleRowMatrix::multTrans
void multTrans(const Vector< T > &x, Vector< T > &b, const T alpha=1.0, const T beta=0.0) const
b <- alpha A&#39;x + beta b
Definition: linalg.h:5486
spams::Matrix::toVect
void toVect(Vector< T > &vec) const
make a reference of the matrix to a vector vec
Definition: linalg.h:2647
spams::SubMatrix::_indicesI
Vector< int > _indicesI
Definition: linalg.h:5182
spams::Matrix::trace
T trace() const
return the trace of the matrix
Definition: linalg.h:2090
spams::Data::getData
virtual void getData(Vector< T > &data, const int i) const =0
spams::abs
T abs(const T x)
template version of the fabs function
spams::SpVector::refIndices
void refIndices(Vector< int > &indices) const
create a reference on the vector r
Definition: linalg.h:4920
spams::Vector::clear
void clear()
clear the vector
Definition: linalg.h:2917
spams::SubMatrix::m
int m() const
Definition: linalg.h:5171
spams::Matrix::norm_2sq_cols
void norm_2sq_cols(Vector< T > &norms) const
returns the l2 norms ^2 of the columns
Definition: linalg.h:2204
MAX_THREADS
#define MAX_THREADS
Definition: utils.h:42
spams::SpVector::print
void print(const string &name) const
print the vector to std::cerr
Definition: linalg.h:4912
spams::DoubleRowMatrix
Matrix with shifts.
Definition: linalg.h:5433
EPSILON
#define EPSILON
Definition: utils.h:60
spams::Matrix::invSym
void invSym()
inverse the matrix when it is symmetric
Definition: linalg.h:1716
spams::Matrix::fmax
int fmax() const
return the 1D-index of the value of greatest magnitude
Definition: linalg.h:1188
spams::SpVector::dot
T dot(const SpVector< T > &vec) const
dot product
Definition: linalg.h:4937
spams::fusedHomotopyAux
static T fusedHomotopyAux(const bool &sign1, const bool &sign2, const bool &sign3, const T &c1, const T &c2)
Definition: linalg.h:3654
INTT
int INTT
Definition: linalg.h:49
spams::Vector
Dense Vector class.
Definition: linalg.h:65
spams::Matrix::setData
void setData(T *X, int m, int n)
Change the data in the matrix.
Definition: linalg.h:1231
spams::Matrix::isNormalized
bool isNormalized() const
Check wether the columns of the matrix are normalized or not.
Definition: linalg.h:1158
spams::Matrix::copy
void copy(const Matrix< T > &mat)
make a copy of the matrix mat in the current matrix
Definition: linalg.h:1339
spams::Matrix::set
void set(const T a)
Set all the values to a scalar.
Definition: linalg.h:1245
spams::SpMatrix::SpMatrix
SpMatrix()
Empty constructor.
Definition: linalg.h:4096
spams::ProdMatrix
Class representing the product of two matrices.
Definition: linalg.h:998
spams::Matrix::softThrshold
void softThrshold(const T nu)
perform soft-thresholding of the matrix, with the threshold nu
Definition: linalg.h:2283
spams::Vector::logspace
void logspace(const int n, const T a, const T b)
generate logarithmically spaced values
Definition: linalg.h:2809
spams::Vector::toSparse
void toSparse(SpVector< T > &vec) const
make a sparse copy
Definition: linalg.h:4025
spams::SpMatrix::copy
void copy(const SpMatrix< T > &mat)
Definition: linalg.h:4102
spams::Vector::thrsmax
void thrsmax(const T nu)
performs soft-thresholding of the vector
Definition: linalg.h:2948
spams::Matrix::getGroup
virtual void getGroup(Matrix< T > &data, const vector_groups &groups, const int i) const
extract the group i
Definition: linalg.h:1131
spams::Matrix::exp
void exp()
each element of the matrix is replaced by its exponential
Definition: linalg.h:2010
spams::Matrix::blockThrshold
void blockThrshold(const T nu, const int sizeGroup)
performs soft-thresholding of the vector
Definition: linalg.h:2320
spams::Vector::fakeSize
void fakeSize(const int n)
change artificially the size of the vector, DANGEROUS
Definition: linalg.h:595
spams::ProdMatrix::extract_rawCol
void extract_rawCol(const int i, T *DtXi) const
compute DtX(:,i)
Definition: linalg.h:5087
spams::SpMatrix::nnz
int nnz() const
number of nonzeros elements
Definition: linalg.h:812
spams::DoubleRowMatrix::_inputmatrix
const AbstractMatrixB< T > * _inputmatrix
Definition: linalg.h:5482
spams::DoubleRowMatrix::n
int n() const
Definition: linalg.h:5439
spams::Vector::thrsmin
void thrsmin(const T nu)
performs thresholding of the vector
Definition: linalg.h:2954
spams::Matrix::fakeSize
void fakeSize(const int m, const int n)
change artificially the size of the matrix, DANGEROUS
Definition: linalg.h:307
spams::Matrix::extractRow
void extractRow(const int i, Vector< T > &row) const
fill the matrix with the row given
Definition: linalg.h:2251
spams::Vector::setPointer
void setPointer(T *X, const int n)
change the data of the vector
Definition: linalg.h:2889
spams::SpMatrix::clear
void clear()
clear the matrix
Definition: linalg.h:4204
spams::DoubleRowMatrix::print
virtual void print(const string &name) const
Definition: linalg.h:5576
spams::SpMatrix::dot
T dot(const Matrix< T > &x) const
dot product;
Definition: linalg.h:4468
spams::Matrix::Invsqrt
void Invsqrt()
Definition: linalg.h:2017
spams::Vector::exp
void exp()
replace each value by its exponential
Definition: linalg.h:3287
spams::Matrix::drop
void drop(char *fileName) const
Definition: linalg.h:2504
spams::Vector::_externAlloc
bool _externAlloc
if the data has been externally allocated
Definition: linalg.h:752
spams::Matrix::_n
int _n
number of columns
Definition: linalg.h:544
spams::Vector::set
void set(const T val)
set each value of the vector to val
Definition: linalg.h:2997
spams::Vector::thrsPos
void thrsPos()
performs soft-thresholding of the vector
Definition: linalg.h:2973
spams::SpVector::sqr
void sqr()
a <- a.^2
Definition: linalg.h:4932
spams::SpVector::_nzmax
int _nzmax
maximum number of nonzeros elements
Definition: linalg.h:993
spams::Matrix::fillRow
void fillRow(const Vector< T > &row)
fill the matrix with the row given
Definition: linalg.h:2241
spams::SpVector::fmaxval
T fmaxval() const
computes the linf norm of the vector
Definition: linalg.h:4906
spams::Vector::max
int max() const
returns the index of the largest value
Definition: linalg.h:2761
spams::Matrix::dot
T dot(const Matrix< T > &mat) const
add alpha*mat to the current matrix
Definition: linalg.h:2068
spams::SpVector::nrm2sq
T nrm2sq() const
computes the l2 norm ^2 of the vector
Definition: linalg.h:4896
spams::SpMatrix::asum
T asum() const
compute the sum of the matrix elements
Definition: linalg.h:4173
spams::Vector::operator[]
T & operator[](const int index)
returns a reference to X[index]
Definition: linalg.h:2853
spams::SpMatrix::m
int m() const
returns the number of columns
Definition: linalg.h:789
spams::Vector::fusedProjectHomotopy
void fusedProjectHomotopy(Vector< T > &out, const T lambda1, const T lambda2, const T lambda3=0, const bool penalty=true)
Definition: linalg.h:3675
spams::AbstractMatrix
Abstract matrix class.
Definition: linalg.h:188
spams::Matrix::thrsmax
void thrsmax(const T nu)
perform thresholding of the matrix, with the threshold nu
Definition: linalg.h:2298
no
static char no
Definition: cblas_alt_template.h:28
spams::ProdMatrix::m
int m() const
returns the number of rows
Definition: linalg.h:1029
spams::Data::getGroup
virtual void getGroup(Matrix< T > &data, const vector_groups &groups, const int i) const =0
spams::Matrix::mult
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
Definition: linalg.h:1783
spams::Matrix::n
int n() const
Number of columns.
Definition: linalg.h:224
spams::Matrix::eigLargestMagnSym
T eigLargestMagnSym() const
Definition: linalg.h:1695
spams::Matrix::fmin
int fmin() const
return the 1D-index of the value of lowest magnitude
Definition: linalg.h:1199
spams::Matrix::refSubMat
void refSubMat(int i, int n, Matrix< T > &mat) const
Reference the column i to i+n into the Matrix mat.
Definition: linalg.h:1153
spams::SpVector::toSpMatrix
void toSpMatrix(SpMatrix< T > &out, const int m, const int n) const
resize the vector as a sparse matrix
Definition: linalg.h:4985
spams::SpVector::_v
T * _v
data
Definition: linalg.h:987
spams::Matrix::normF
T normF() const
return ||A||_F
Definition: linalg.h:2099
spams::ProdMatrix::diag
void diag(Vector< T > &diag) const
add something to the diagonal
Definition: linalg.h:5148
spams::Matrix::singularValues
void singularValues(Vector< T > &u) const
Definition: linalg.h:1560
spams::Matrix::copyRef
void copyRef(const Matrix< T > &mat)
make a copy of the matrix mat in the current matrix
Definition: linalg.h:1345
spams::SpMatrix::norm_2sq_cols
void norm_2sq_cols(Vector< T > &norms) const
returns the l2 norms ^2 of the columns
Definition: linalg.h:4832
spams::alt_log< float >
float alt_log< float >(const float x)
Definition: linalg.h:105
incr
static char incr
Definition: cblas_alt_template.h:25
spams::SpVector::nrm2
T nrm2() const
computes the l2 norm of the vector
Definition: linalg.h:4901
spams::DoubleRowMatrix::multSwitch
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
perform C = a*B*A + b*C, possibly transposing A or B.
Definition: linalg.h:5547
spams::SpVector::_r
int * _r
indices
Definition: linalg.h:989
spams::SubMatrix::n
int n() const
Definition: linalg.h:5170
spams::SpVector
Sparse Vector class.
Definition: linalg.h:67
spams::SpVector::clear
void clear()
clears the vector
Definition: linalg.h:4958
spams::Matrix::sparseProject
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)
performs sparse projections of the columns
Definition: linalg.h:2344
spams::quick_sort
static void quick_sort(int *irOut, T *prOut, const int beg, const int end, const bool incr)
Definition: misc.h:171
spams::Matrix::operator()
T & operator()(const int i, const int j)
Return a modifiable reference to X(i,j)
Definition: linalg.h:1075
spams::Vector::nrm2sq
T nrm2sq() const
returns ||A||_2^2
Definition: linalg.h:3007
spams::SpMatrix::_r
int * _r
row indices
Definition: linalg.h:904
utl::cblas_asum
T cblas_asum(const INTT N, const T *X, const INTT incX)
spams::Matrix::setZeros
void setZeros()
Set all the values to zero.
Definition: linalg.h:1240
spams::Vector::print
void print(const char *name) const
spams::Vector::resize
void resize(const int n)
resize the vector
Definition: linalg.h:2876
spams::Data::norm_2sq_cols
virtual void norm_2sq_cols(Vector< T > &norms) const
Definition: linalg.h:139
spams::SpMatrix::V
int V() const
returns the number of columns
Definition: linalg.h:791
spams::isEqual
static bool isEqual(const T lambda1, const T lambda2)
Definition: linalg.h:81
spams::Matrix::XXt
void XXt(Matrix< T > &XXt) const
XXt = A*A&#39;.
Definition: linalg.h:1967
utl::Vector
NDArray< T, 1 > Vector
Definition: utlNDArray.h:44
spams::Matrix::clean
void clean()
clean a dictionary matrix
Definition: linalg.h:1168
spams::SpMatrix::r
int r(const int i) const
returns r[i]
Definition: linalg.h:781
spams::SpMatrix::convert2
void convert2(const Matrix< T > &v, const Vector< int > &r, const int K)
use the data from v, r for _v, _r
Definition: linalg.h:4809
spams::SpMatrix::XAt
void XAt(const Matrix< T > &X, Matrix< T > &XAt) const
XAt <- X*A&#39;.
Definition: linalg.h:4654
spams::ProdMatrix::_D
const Matrix< T > * _D
Definition: linalg.h:1039
spams::SpMatrix::resize
void resize(const int m, const int n, const int nzmax)
resize the matrix
Definition: linalg.h:4221
spams::Matrix::nrm2sq
T nrm2sq() const
return ||A||_F^2
Definition: linalg.h:446
spams::Matrix::svd
void svd(Matrix< T > &U, Vector< T > &S, Matrix< T > &V) const
Definition: linalg.h:1582
spams::Matrix::_m
int _m
number of rows
Definition: linalg.h:542
spams::Vector::softThrshold
void softThrshold(const T nu)
performs soft-thresholding of the vector
Definition: linalg.h:2925
spams::Matrix::Sqrt
void Sqrt()
each element of the matrix is replaced by its square root
Definition: linalg.h:2013
spams::SpMatrix::addVecToCols
void addVecToCols(const Vector< T > &diag, const T a=1.0)
Definition: linalg.h:4494
spams::Matrix::norm_inf_cols
void norm_inf_cols(Vector< T > &norms) const
returns the linf norms of the columns
Definition: linalg.h:2174
spams::SpVector::refVal
void refVal(Vector< T > &val) const
creates a reference on the vector val
Definition: linalg.h:4926
spams::SpMatrix::_v
T * _v
data
Definition: linalg.h:902
spams::ShiftMatrix::~ShiftMatrix
virtual ~ShiftMatrix()
Definition: linalg.h:5293
spams::Vector::hardThrshold
void hardThrshold(const T nu)
performs soft-thresholding of the vector
Definition: linalg.h:2938
spams::Matrix::meanRow
void meanRow(Vector< T > &mean) const
Compute the mean of the rows.
Definition: linalg.h:2233
spams::alt_log< double >
double alt_log< double >(const double x)
Definition: linalg.h:104
spams::Matrix::upperTriXXt
void upperTriXXt(Matrix< T > &XXt, const int L) const
XXt = A*A&#39; where A is an upper triangular matrix.
Definition: linalg.h:1975
spams::SpMatrix::getGroup
void getGroup(Matrix< T > &data, const vector_groups &groups, const int i) const
Definition: linalg.h:4159
cblas_template.h
spams::DoubleRowMatrix::~DoubleRowMatrix
virtual ~DoubleRowMatrix()
Definition: linalg.h:5477
spams::logexp
static T logexp(const T x)
Definition: linalg.h:119
spams::Matrix::eigLargestSymApprox
void eigLargestSymApprox(const Vector< T > &u0, Vector< T > &u) const
Definition: linalg.h:1629
spams::Vector::fmax
int fmax() const
returns the index of the value with largest magnitude
Definition: linalg.h:2843
spams::SpMatrix::getData
void getData(Vector< T > &data, const int index) const
Definition: linalg.h:4151
spams::Matrix::whiten
void whiten(const int V)
whiten
Definition: linalg.h:1367
spams::SpMatrix::pB
int * pB() const
Direct access to _pB.
Definition: linalg.h:804
spams::isZero
static bool isZero(const T lambda)
Definition: linalg.h:76
spams::Matrix::extract_rawCol
void extract_rawCol(const int i, T *x) const
Copy the column i into x.
Definition: linalg.h:1115
spams::SpVector::resize
void resize(const int nzmax)
resizes the vector
Definition: linalg.h:4971
spams::Vector::minval
T minval() const
returns the minimum value
Definition: linalg.h:2794
spams::init_omp
static int init_omp(const int numThreads)
Definition: misc.h:264
spams::SubMatrix::SubMatrix
SubMatrix(AbstractMatrix< T > &G, Vector< int > &indI, Vector< int > &indJ)
Definition: linalg.h:5188
spams::Vector::fminval
T fminval() const
returns the minimum magnitude
Definition: linalg.h:2804
spams::DoubleRowMatrix::DoubleRowMatrix
DoubleRowMatrix(const AbstractMatrixB< T > &inputmatrix)
Definition: linalg.h:5435
spams::Vector::fusedProject
void fusedProject(Vector< T > &out, const T lambda1, const T lambda2, const int itermax)
Definition: linalg.h:3846
spams::Matrix::norm_inf_2_col
T norm_inf_2_col() const
return ||At||_{inf,2} (max of l2 norm of the columns)
Definition: linalg.h:2115
spams::ShiftMatrix::_inputmatrix
const AbstractMatrixB< T > * _inputmatrix
Definition: linalg.h:5307
spams::ProdMatrix::setMatrices
void setMatrices(const Matrix< T > &D, const bool high_memory=true)
set_matrices
Definition: linalg.h:5059
spams::SpMatrix::n
int n() const
returns the number of rows
Definition: linalg.h:787
spams::ShiftMatrix::dot
virtual T dot(const Matrix< T > &x) const
Definition: linalg.h:5421
spams::Matrix::inv
void inv()
inverse the elements of the matrix
Definition: linalg.h:434
spams::Matrix::eye
void eye()
set the matrix to the identity;
Definition: linalg.h:1264
spams::DoubleRowMatrix::copyTo
virtual void copyTo(Matrix< T > &copy) const
Definition: linalg.h:5564
spams::Matrix::asum
T asum() const
compute the sum of the magnitude of the matrix values
Definition: linalg.h:2085
spams::Matrix::mean
T mean() const
whiten
Definition: linalg.h:2103
spams::Matrix::multTrans
void multTrans(const Vector< T > &x, Vector< T > &b, const T alpha=1.0, const T beta=0.0) const
perform b = alpha*A&#39;x + beta*b
Definition: linalg.h:1743
spams::ProdMatrix::operator[]
T operator[](const int index) const
returns the value of an index
Definition: linalg.h:5122
spams::Vector::setAleat
void setAleat()
put random values in the vector (White Gaussian Noise)
Definition: linalg.h:2912
spams::list_groups
std::list< group > list_groups
Definition: linalg.h:72
spams::SpVector::toFull
void toFull(Vector< T > &out) const
resize the vector as a sparse matrix
Definition: linalg.h:5008
spams::SpMatrix::convert
void convert(const Matrix< T > &v, const Matrix< int > &r, const int K)
use the data from v, r for _v, _r
Definition: linalg.h:4786
spams::Matrix::norm_2_rows
void norm_2_rows(Vector< T > &norms) const
returns the l2 norms of the columns
Definition: linalg.h:2139
spams::Matrix::rank1Update_mult
void rank1Update_mult(const Vector< T > &vec1, const Vector< T > &vec1b, const SpVector< T > &vec2, const T alpha=1.0)
Definition: linalg.h:2412
spams::SubMatrix::copyCol
void copyCol(const int i, Vector< T > &DtXi) const
compute DtX(:,i)
Definition: linalg.h:5220
spams::Vector::~Vector
virtual ~Vector()
Destructor.
Definition: linalg.h:2720
spams::Matrix::mult_elementWise
void mult_elementWise(const Matrix< T > &B, Matrix< T > &C) const
C = A .* B, elementwise multiplication.
Definition: linalg.h:1942
spams::Vector::div
void div(const Vector< T > &x)
A <- A ./ x.
Definition: linalg.h:3064
spams::SpMatrix::norm_1_cols
void norm_1_cols(Vector< T > &norms) const
returns the l1 norms of the columns
Definition: linalg.h:4852
spams::ShiftMatrix::n
int n() const
Definition: linalg.h:5255
spams::AbstractMatrixB
Abstract matrix class.
Definition: linalg.h:144
utl::sign
int sign(const T &x)
Return the sign of x.
Definition: utlCore.h:285
spams::Matrix::resize
void resize(int m, int n)
Resize the matrix.
Definition: linalg.h:1217
spams::Matrix::setRow
void setRow(const int i, const Vector< T > &row)
fill the matrix with the row given
Definition: linalg.h:2260
spams::Matrix::normFsq
T normFsq() const
return ||A||_F^2
Definition: linalg.h:2110
spams::SpVector::v
T v(const int i) const
access table r
Definition: linalg.h:955
spams::Vector::l1l2project
void l1l2project(Vector< T > &out, const T thrs, const T gamma, const bool pos=false) const
Definition: linalg.h:3593
spams::Matrix::thrsPos
void thrsPos()
perform soft-thresholding of the matrix, with the threshold nu
Definition: linalg.h:2376
spams::SpMatrix::sum_cols
void sum_cols(Vector< T > &sum) const
Definition: linalg.h:4524
spams::Matrix::toSparse
void toSparse(SpMatrix< T > &matrix) const
make a sparse copy of the current matrix
Definition: linalg.h:2566
spams::Matrix::X
const T * X() const
return a non-modifiable reference to the data
Definition: linalg.h:256
spams::SpMatrix::multTrans
void multTrans(const Vector< T > &x, Vector< T > &y, const T alpha=1.0, const T beta=0.0) const
y <- A&#39;*x
Definition: linalg.h:4246
spams::Vector::sub
void sub(const Vector< T > &x)
A <- A - x.
Definition: linalg.h:3052
spams::SpMatrix::toFullTrans
void toFullTrans(Matrix< T > &matrix) const
copy the sparse matrix into a dense transposed matrix
Definition: linalg.h:4772
spams::Matrix
Dense Matrix class.
Definition: linalg.h:61
spams::Matrix::NadarayaWatson
void NadarayaWatson(const Vector< int > &ind, const T sigma)
compute a Nadaraya Watson estimator
Definition: linalg.h:2519
spams::Matrix::V
int V() const
Accessor.
Definition: linalg.h:525
spams::Matrix::operator[]
T & operator[](const int index)
Return a modifiable reference to X(i) (1D indexing)
Definition: linalg.h:230
spams::SubMatrix::extract_rawCol
void extract_rawCol(const int i, T *pr) const
copy X(:,i) into Xi
Definition: linalg.h:5212
spams::Vector::logexp
void logexp()
replace each value by its exponential
Definition: linalg.h:3297
spams::ShiftMatrix::copyRow
virtual void copyRow(const int i, Vector< T > &x) const
Definition: linalg.h:5406
spams::Vector::mean
T mean()
whiten
Definition: linalg.h:3260
spams::Vector::lzero
T lzero() const
Definition: linalg.h:3328
spams::DoubleRowMatrix::m
int m() const
Definition: linalg.h:5440
spams::Vector::randperm
void randperm(int n)
put a random permutation of size n (for integral vectors)
spams::DoubleRowMatrix::copyRow
virtual void copyRow(const int i, Vector< T > &x) const
Definition: linalg.h:5559
spams::AbstractMatrix::~AbstractMatrix
virtual ~AbstractMatrix()
Definition: linalg.h:202
spams::SpVector::~SpVector
~SpVector()
Destructor.
Definition: linalg.h:4887
spams::Vector::setData
void setData(T *X, const int n)
Definition: linalg.h:607
spams::Data::n
virtual int n() const =0
spams::ProdMatrix::addDiag
void addDiag(const T diag)
add something to the diagonal
Definition: linalg.h:5112
spams::ShiftMatrix::ShiftMatrix
ShiftMatrix(const AbstractMatrixB< T > &inputmatrix, const int shifts, const bool center=false)
Definition: linalg.h:5250
spams::Matrix::setm
void setm(const int m)
modify _m
Definition: linalg.h:277
spams::Matrix::add
void add(const Matrix< T > &mat, const T alpha=1.0)
add alpha*mat to the current matrix
Definition: linalg.h:2062
spams::SpMatrix::mult
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
Definition: linalg.h:4305
spams::Vector::Sqrt
void Sqrt()
A <- 1 ./ sqrt(x)
Definition: linalg.h:3097
spams::Vector::l1project_weighted
void l1project_weighted(Vector< T > &out, const Vector< T > &weights, const T thrs, const bool residual=false) const
Definition: linalg.h:3417
INFINITY
#define INFINITY
Definition: utils.h:62
decr
static char decr
Definition: cblas_alt_template.h:26
spams::Vector::l1project
void l1project(Vector< T > &out, const T thrs, const bool simplex=false) const
Definition: linalg.h:3364
spams::SpVector::r
int r(const int i) const
access table r
Definition: linalg.h:953
spams::SpVector::_L
int _L
length
Definition: linalg.h:991
spams::Matrix::scal
void scal(const T a)
scale the matrix by the a
Definition: linalg.h:1334
spams::Vector::project_sft
void project_sft(const Vector< int > &labels, const int clas)
Definition: linalg.h:3482
spams::Matrix::rawX
T * rawX() const
reference a modifiable reference to the data, DANGEROUS
Definition: linalg.h:254
spams::SpMatrix::_m
int _m
number of rows
Definition: linalg.h:910
spams::Matrix::clear
void clear()
Clear the matrix.
Definition: linalg.h:1250
Generated by   doxygen 1.8.11