bpp-core  2.2.0
MatrixTools.h
Go to the documentation of this file.
1 //
2 // File: MatrixTools.h
3 // Created by: Julien Dutheil
4 // Created on: Mon Jan 19 16:42:25 2004
5 //
6 
7 /*
8  Copyright or © or Copr. Bio++ Development Team, (November 17, 2004)
9 
10  This software is a computer program whose purpose is to provide classes
11  for numerical calculus. This file is part of the Bio++ project.
12 
13  This software is governed by the CeCILL license under French law and
14  abiding by the rules of distribution of free software. You can use,
15  modify and/ or redistribute the software under the terms of the CeCILL
16  license as circulated by CEA, CNRS and INRIA at the following URL
17  "http://www.cecill.info".
18 
19  As a counterpart to the access to the source code and rights to copy,
20  modify and redistribute granted by the license, users are provided only
21  with a limited warranty and the software's author, the holder of the
22  economic rights, and the successive licensors have only limited
23  liability.
24 
25  In this respect, the user's attention is drawn to the risks associated
26  with loading, using, modifying and/or developing or reproducing the
27  software by the user in light of its specific status of free software,
28  that may mean that it is complicated to manipulate, and that also
29  therefore means that it is reserved for developers and experienced
30  professionals having in-depth computer knowledge. Users are therefore
31  encouraged to load and test the software's suitability as regards their
32  requirements in conditions enabling the security of their systems and/or
33  data to be ensured and, more generally, to use and operate it in the
34  same conditions as regards security.
35 
36  The fact that you are presently reading this means that you have had
37  knowledge of the CeCILL license and that you accept its terms.
38 */
39 
40 #ifndef _MATRIXTOOLS_H_
41 #define _MATRIXTOOLS_H_
42 
43 #include "../VectorTools.h"
44 #include "Matrix.h"
45 #include "LUDecomposition.h"
46 #include "EigenValue.h"
47 #include "../../Io/OutputStream.h"
48 
49 #include <cstdio>
50 #include <iostream>
51 
52 namespace bpp
53 {
58  {
59  public:
62 
63  public:
70  template<class MatrixA, class MatrixO>
71  static void copy(const MatrixA& A, MatrixO& O)
72  {
73  O.resize(A.getNumberOfRows(), A.getNumberOfColumns());
74  for (size_t i = 0; i < A.getNumberOfRows(); i++)
75  {
76  for (size_t j = 0; j < A.getNumberOfColumns(); j++)
77  {
78  O(i, j) = A(i, j);
79  }
80  }
81  }
82 
89  template<class Matrix>
90  static void getId(size_t n, Matrix& O)
91  {
92  O.resize(n, n);
93  for (size_t i = 0; i < n; i++)
94  {
95  for (size_t j = 0; j < n; j++) {
96  O(i, j) = (i == j) ? 1 : 0;
97  }
98  }
99  }
100 
105  template<class Scalar>
106  static void diag(const std::vector<Scalar>& D, Matrix<Scalar>& O)
107  {
108  size_t n = D.size();
109  O.resize(n, n);
110  for (size_t i = 0; i < n; i++)
111  {
112  for (size_t j = 0; j < n; j++) { O(i, j) = (i == j) ? D[i] : 0;}
113  }
114  }
115 
121  template<class Scalar>
122  static void diag(const Scalar x, size_t n, Matrix<Scalar>& O)
123  {
124  O.resize(n, n);
125  for (size_t i = 0; i < n; i++)
126  {
127  for (size_t j = 0; j < n; j++) { O(i, j) = (i == j) ? x : 0;}
128  }
129  }
130 
136  template<class Scalar>
137  static void diag(const Matrix<Scalar>& M, std::vector<Scalar>& O) throw (DimensionException)
138  {
139  size_t nc = M.getNumberOfColumns();
140  size_t nr = M.getNumberOfRows();
141  if (nc != nr) throw DimensionException("MatrixTools::diag(). M must be a square matrix.", nr, nc);
142  O.resize(nc);
143  for (size_t i = 0; i < nc; i++) { O[i] = M(i, i);}
144  }
145 
151  template<class Matrix, class Scalar>
152  static void fill(Matrix& M, Scalar x)
153  {
154  for (size_t i = 0; i < M.getNumberOfRows(); i++)
155  {
156  for (size_t j = 0; j < M.getNumberOfColumns(); j++)
157  {
158  M(i, j) = x;
159  }
160  }
161  }
162 
172  template<class Matrix, class Scalar>
173  static void scale(Matrix& A, Scalar a, Scalar b = 0)
174  {
175  for (size_t i = 0; i < A.getNumberOfRows(); i++)
176  {
177  for (size_t j = 0; j < A.getNumberOfColumns(); j++)
178  {
179  A(i, j) = a * A(i, j) + b;
180  }
181  }
182  }
183 
189  template<class Scalar>
190  static void mult(const Matrix<Scalar>& A, const Matrix<Scalar>& B, Matrix<Scalar>& O) throw (DimensionException)
191  {
192  size_t ncA = A.getNumberOfColumns();
193  size_t nrA = A.getNumberOfRows();
194  size_t nrB = B.getNumberOfRows();
195  size_t ncB = B.getNumberOfColumns();
196  if (ncA != nrB) throw DimensionException("MatrixTools::mult(). nrows B != ncols A.", nrB, ncA);
197  O.resize(nrA, ncB);
198  for (size_t i = 0; i < nrA; i++)
199  {
200  for (size_t j = 0; j < ncB; j++)
201  {
202  O(i, j) = 0;
203  for (size_t k = 0; k < ncA; k++)
204  {
205  O(i, j) += A(i, k) * B(k, j);
206  }
207  }
208  }
209  }
210 
223  template<class Scalar>
224  static void mult(const Matrix<Scalar>& A, const std::vector<Scalar>& D, const Matrix<Scalar>& B, Matrix<Scalar>& O) throw (DimensionException)
225  {
226  size_t ncA = A.getNumberOfColumns();
227  size_t nrA = A.getNumberOfRows();
228  size_t nrB = B.getNumberOfRows();
229  size_t ncB = B.getNumberOfColumns();
230  if (ncA != nrB) throw DimensionException("MatrixTools::mult(). nrows B != ncols A.", nrB, ncA);
231  if (ncA != D.size()) throw DimensionException("MatrixTools::mult(). Vector size is not equal to matrix size.", D.size(), ncA);
232  O.resize(nrA, ncB);
233  for (size_t i = 0; i < nrA; i++)
234  {
235  for (size_t j = 0; j < ncB; j++)
236  {
237  O(i, j) = 0;
238  for (size_t k = 0; k < ncA; k++)
239  {
240  O(i, j) += A(i, k) * B(k, j) * D[k];
241  }
242  }
243  }
244  }
245 
262  template<class Scalar>
263  static void mult(const Matrix<Scalar>& A, const std::vector<Scalar>& D, const std::vector<Scalar>& U, const std::vector<Scalar>& L, const Matrix<Scalar>& B, Matrix<Scalar>& O) throw (DimensionException)
264  {
265  size_t ncA = A.getNumberOfColumns();
266  size_t nrA = A.getNumberOfRows();
267  size_t nrB = B.getNumberOfRows();
268  size_t ncB = B.getNumberOfColumns();
269  if (ncA != nrB) throw DimensionException("MatrixTools::mult(). nrows B != ncols A.", nrB, ncA);
270  if (ncA != D.size()) throw DimensionException("MatrixTools::mult(). Vector size is not equal to matrix size.", D.size(), ncA);
271  if (ncA != U.size()+1) throw DimensionException("MatrixTools::mult(). Vector size is not equal to matrix size-1.", U.size(), ncA);
272  if (ncA != L.size()+1) throw DimensionException("MatrixTools::mult(). Vector size is not equal to matrix size-1.", L.size(), ncA);
273  O.resize(nrA, ncB);
274  for (size_t i = 0; i < nrA; i++)
275  {
276  for (size_t j = 0; j < ncB; j++)
277  {
278  O(i, j) = A(i, 0) * D[0] * B(0, j);
279  if (nrB>1)
280  O(i, j) += A(i,0) * U[0] * B(1,j);
281  for (size_t k = 1; k < ncA-1; k++)
282  {
283  O(i, j) += A(i, k) * (L[k-1] * B(k-1, j) + D[k] * B(k, j) + U[k] * B(k+1,j));
284  }
285  if (ncA>=2)
286  O(i, j) += A(i, ncA-1) * L[ncA-2] * B(ncA-2, j);
287  O(i,j) += A(i, ncA-1) * D[ncA-1] * B(ncA-1, j);
288  }
289  }
290  }
291 
299  template<class MatrixA, class MatrixB>
300  static void add(MatrixA& A, const MatrixB& B) throw (DimensionException)
301  {
302  size_t ncA = A.getNumberOfColumns();
303  size_t nrA = A.getNumberOfRows();
304  size_t nrB = B.getNumberOfRows();
305  size_t ncB = B.getNumberOfColumns();
306  if (ncA != ncB) throw DimensionException("MatrixTools::operator+(). A and B must have the same number of colums.", ncB, ncA);
307  if (nrA != nrB) throw DimensionException("MatrixTools::operator+(). A and B must have the same number of rows.", nrB, nrA);
308  for (size_t i = 0; i < A.getNumberOfRows(); i++)
309  {
310  for (size_t j = 0; j < A.getNumberOfColumns(); j++)
311  {
312  A(i, j) += B(i, j);
313  }
314  }
315  }
316 
325  template<class MatrixA, class MatrixB, class Scalar>
326  static void add(MatrixA& A, Scalar& x, const MatrixB& B) throw (DimensionException)
327  {
328  size_t ncA = A.getNumberOfColumns();
329  size_t nrA = A.getNumberOfRows();
330  size_t nrB = B.getNumberOfRows();
331  size_t ncB = B.getNumberOfColumns();
332  if (ncA != ncB) throw DimensionException("MatrixTools::operator+(). A and B must have the same number of colums.", ncB, ncA);
333  if (nrA != nrB) throw DimensionException("MatrixTools::operator+(). A and B must have the same number of rows.", nrB, nrA);
334  for (size_t i = 0; i < A.getNumberOfRows(); i++)
335  {
336  for (size_t j = 0; j < A.getNumberOfColumns(); j++)
337  {
338  A(i, j) += x*B(i, j);
339  }
340  }
341  }
342 
354  template<class Matrix>
355  static void pow(const Matrix& A, size_t p, Matrix& O) throw (DimensionException)
356  {
357  size_t n = A.getNumberOfRows();
358  if (n != A.getNumberOfColumns()) throw DimensionException("MatrixTools::pow(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
359  switch(p){
360  case 0:
361  getId<Matrix>(n, O);
362  break;
363  case 1:
364  copy(A,O);
365  break;
366  case 2:
367  mult(A,A,O);
368  break;
369  default:
370  Matrix tmp;
371  if (p%2){
372  pow(A,p/2,tmp);
373  pow(tmp,2,O);
374  }
375  else{
376  pow(A,(p-1)/2,tmp);
377  pow(tmp,2,O);
378  mult(A,O,tmp);
379  copy(tmp,O);
380  }
381  }
382  }
383 
393  template<class Scalar>
394  static void pow(const Matrix<Scalar>& A, double p, Matrix<Scalar>& O) throw (DimensionException)
395  {
396  size_t n = A.getNumberOfRows();
397  if (n != A.getNumberOfColumns()) throw DimensionException("MatrixTools::pow(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
398  EigenValue<Scalar> eigen(A);
399  RowMatrix<Scalar> rightEV, leftEV;
400  rightEV = eigen.getV();
401  inv(rightEV, leftEV);
402  mult(rightEV, VectorTools::pow(eigen.getRealEigenValues(), p), leftEV, O);
403  }
404 
415  template<class Scalar>
416  static void exp(const Matrix<Scalar>& A, Matrix<Scalar>& O) throw (DimensionException)
417  {
418  size_t n = A.getNumberOfRows();
419  if (n != A.getNumberOfColumns()) throw DimensionException("MatrixTools::exp(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
420  EigenValue<Scalar> eigen(A);
421  RowMatrix<Scalar> rightEV, leftEV;
422  rightEV = eigen.getV();
423  inv(rightEV, leftEV);
424  mult(rightEV, VectorTools::exp(eigen.getRealEigenValues()), leftEV, O);
425  }
426 
437  template<class Matrix, class Scalar>
438  static void Taylor(const Matrix& A, size_t p, std::vector< RowMatrix<Scalar> > & vO) throw (DimensionException)
439  {
440  size_t n = A.getNumberOfRows();
441  if (n != A.getNumberOfColumns())
442  throw DimensionException("MatrixTools::pow(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
443  vO.resize(p+1);
444  getId<Matrix>(n, vO[0]);
445  copy(A,vO[1]);
446 
447  for (size_t i = 1; i < p; i++)
448  {
449  mult(vO[i], A, vO[i+1]);
450  }
451  }
452 
457  template<class Matrix>
458  static std::vector<size_t> whichMax(const Matrix& m)
459  {
460  size_t nrows = m.getNumberOfRows();
461  size_t ncols = m.getNumberOfColumns();
462  std::vector<size_t> pos(2);
463  size_t imax = 0;
464  size_t jmax = 0;
465  double currentMax = log(0.);
466  for (size_t i = 0; i < nrows; i++)
467  {
468  for (size_t j = 0; j < ncols; j++)
469  {
470  double currentValue = m(i, j);
471  // cout << currentValue << "\t" << (currentValue > currentMax) << endl;
472  if (currentValue > currentMax)
473  {
474  imax = i;
475  jmax = j;
476  currentMax = currentValue;
477  }
478  }
479  }
480  pos[0] = imax;
481  pos[1] = jmax;
482  return pos;
483  }
484 
489  template<class Matrix>
490  static std::vector<size_t> whichMin(const Matrix& m)
491  {
492  size_t nrows = m.getNumberOfRows();
493  size_t ncols = m.getNumberOfColumns();
494  std::vector<size_t> pos(2);
495  size_t imin = 0;
496  size_t jmin = 0;
497  double currentMin = -log(0.);
498  for (size_t i = 0; i < nrows; i++)
499  {
500  for (size_t j = 0; j < ncols; j++)
501  {
502  double currentValue = m(i, j);
503  if (currentValue < currentMin)
504  {
505  imin = i;
506  jmin = j;
507  currentMin = currentValue;
508  }
509  }
510  }
511  pos[0] = imin;
512  pos[1] = jmin;
513  return pos;
514  }
515 
520  template<class Real>
521  static Real max(const Matrix<Real>& m)
522  {
523  size_t nrows = m.getNumberOfRows();
524  size_t ncols = m.getNumberOfColumns();
525  Real currentMax = log(0.);
526  for (size_t i = 0; i < nrows; i++)
527  {
528  for (size_t j = 0; j < ncols; j++)
529  {
530  Real currentValue = m(i, j);
531  // cout << currentValue << "\t" << (currentValue > currentMax) << endl;
532  if (currentValue > currentMax)
533  {
534  currentMax = currentValue;
535  }
536  }
537  }
538  return currentMax;
539  }
540 
541 
546  template<class Real>
547  static Real min(const Matrix<Real>& m)
548  {
549  size_t nrows = m.getNumberOfRows();
550  size_t ncols = m.getNumberOfColumns();
551  Real currentMin = -log(0.);
552  for (size_t i = 0; i < nrows; i++)
553  {
554  for (size_t j = 0; j < ncols; j++)
555  {
556  Real currentValue = m(i, j);
557  if (currentValue < currentMin)
558  {
559  currentMin = currentValue;
560  }
561  }
562  }
563  return currentMin;
564  }
565 
572  template<class Matrix>
573  static void print(const Matrix& m, std::ostream& out = std::cout)
574  {
575  out << m.getNumberOfRows() << "x" << m.getNumberOfColumns() << std::endl;
576  out << "[" << std::endl;
577  for (size_t i = 0; i < m.getNumberOfRows(); i++)
578  {
579  out << "[";
580  for (size_t j = 0; j < m.getNumberOfColumns() - 1; j++)
581  {
582  out << m(i, j) << ", ";
583  }
584  if (m.getNumberOfColumns() > 0) out << m(i, m.getNumberOfColumns() - 1) << "]" << std::endl;
585  }
586  out << "]" << std::endl;
587  }
588 
597  template<class Matrix>
598  static void print(const Matrix& m, bpp::OutputStream& out, char pIn = '(', char pOut = ')')
599  {
600  out << pIn;
601 
602  for (size_t i = 0; i < m.getNumberOfRows(); i++)
603  {
604  if (i!=0)
605  out << ",";
606 
607  out << pIn;
608  for (size_t j = 0; j < m.getNumberOfColumns() - 1; j++)
609  {
610  out << m(i, j) << ", ";
611  }
612  if (m.getNumberOfColumns() > 0) out << m(i, m.getNumberOfColumns() - 1) << pOut;
613  }
614  out << pOut;
615  }
616 
624  template<class Matrix>
625  static void printForR(const Matrix& m, const std::string& variableName = "x", std::ostream& out = std::cout)
626  {
627  out.precision(12);
628  out << variableName << "<-matrix(c(";
629  for (size_t i = 0; i < m.getNumberOfRows(); i++)
630  {
631  for (size_t j = 0; j < m.getNumberOfColumns(); j++)
632  {
633  if (i > 0 || j > 0)
634  out << ", ";
635  out << m(i, j);
636  }
637  }
638  out << "), nrow=" << m.getNumberOfRows() << ", byrow=TRUE)" << std::endl;
639  }
640 
641 
648  template<class Real>
649  static void print(const std::vector<Real>& v, std::ostream& out = std::cout)
650  {
651  out << v.size() << std::endl;
652  out << "[";
653  for (size_t i = 0; i < v.size() - 1; i++)
654  {
655  out << v[i] << ", ";
656  }
657  if (v.size() > 0) out << v[v.size() - 1];
658  out << "]" << std::endl;
659  }
660 
665  template<class Matrix>
666  static bool isSquare(const Matrix& A) { return A.getNumberOfRows() == A.getNumberOfColumns(); }
667 
674  template<class Scalar>
676  {
677  if (!isSquare(A)) throw DimensionException("MatrixTools::inv(). Matrix A is not a square matrix.", A.getNumberOfRows(), A.getNumberOfColumns());
680  getId(A.getNumberOfRows(), I);
681  return lu.solve(I, O);
682  }
683 
693  template<class Scalar>
694  static double det(const Matrix<Scalar>& A) throw (DimensionException)
695  {
696  if (!isSquare(A)) throw DimensionException("MatrixTools::det(). Matrix A is not a square matrix.", A.getNumberOfRows(), A.getNumberOfColumns());
698  return lu.det();
699  }
700 
705  template<class MatrixA, class MatrixO>
706  static void transpose(const MatrixA& A, MatrixO& O)
707  {
708  O.resize(A.getNumberOfColumns(), A.getNumberOfRows());
709  for (size_t i = 0; i < A.getNumberOfColumns(); i++)
710  {
711  for (size_t j = 0; j < A.getNumberOfRows(); j++)
712  {
713  O(i, j) = A(j, i);
714  }
715  }
716  }
717 
730  template<class Scalar>
731  static void covar(const Matrix<Scalar>& A, Matrix<Scalar>& O)
732  {
733  size_t r = A.getNumberOfRows();
734  size_t n = A.getNumberOfColumns();
735  O.resize(r, r);
737  transpose(A, tA);
738  mult(A, tA, O);
739  scale(O, 1. / static_cast<double>(n));
740  RowMatrix<Scalar> mean(r, 1);
741  for (size_t i = 0; i < r; i++)
742  {
743  for (size_t j = 0; j < n; j++)
744  {
745  mean(i, 0) += A(i, j);
746  }
747  mean(i, 0) /= static_cast<double>(n);
748  }
749  RowMatrix<Scalar> tMean;
750  transpose(mean, tMean);
751  RowMatrix<Scalar> meanMat;
752  mult(mean, tMean, meanMat);
753  scale(meanMat, -1.);
754  add(O, meanMat);
755  }
756 
764  template<class Scalar>
765  static void kroneckerMult(const Matrix<Scalar>& A, const Matrix<Scalar>& B, Matrix<Scalar>& O)
766  {
767  size_t ncA = A.getNumberOfColumns();
768  size_t nrA = A.getNumberOfRows();
769  size_t nrB = B.getNumberOfRows();
770  size_t ncB = B.getNumberOfColumns();
771  O.resize(nrA * nrB, ncA * ncB);
772  for (size_t ia = 0; ia < nrA; ia++)
773  {
774  for (size_t ja = 0; ja < ncA; ja++)
775  {
776  Scalar aij = A(ia, ja);
777  for (size_t ib = 0; ib < nrB; ib++)
778  {
779  for (size_t jb = 0; jb < ncB; jb++)
780  {
781  O(ia * nrB + ib, ja * ncB + jb) = aij * B(ib, jb);
782  }
783  }
784  }
785  }
786  }
787 
795  template<class Scalar>
796  static void hadamardMult(const Matrix<Scalar>& A, const Matrix<Scalar>& B, Matrix<Scalar>& O)
797  {
798  size_t ncA = A.getNumberOfColumns();
799  size_t nrA = A.getNumberOfRows();
800  size_t nrB = B.getNumberOfRows();
801  size_t ncB = B.getNumberOfColumns();
802  if (nrA != nrB) throw DimensionException("MatrixTools::hadamardMult(). nrows A != nrows B.", nrA, nrB);
803  if (ncA != ncB) throw DimensionException("MatrixTools::hadamardMult(). ncols A != ncols B.", ncA, ncB);
804  O.resize(nrA, ncA);
805  for (size_t i = 0; i < nrA; i++)
806  {
807  for (size_t j = 0; j < ncA; j++)
808  {
809  O(i, j) = A(i, j) * B(i, j);
810  }
811  }
812  }
813 
822  template<class Scalar>
823  static void hadamardMult(const Matrix<Scalar>& A, const std::vector<Scalar>& B, Matrix<Scalar>& O, bool row = true)
824  {
825  size_t ncA = A.getNumberOfColumns();
826  size_t nrA = A.getNumberOfRows();
827  size_t sB = B.size();
828  if (row == true && nrA != sB) throw DimensionException("MatrixTools::hadamardMult(). nrows A != size of B.", nrA, sB);
829  if (row == false && ncA != sB) throw DimensionException("MatrixTools::hadamardMult(). ncols A != size of B.", ncA, sB);
830  O.resize(nrA, ncA);
831  if (row)
832  {
833  for (size_t i = 0; i < nrA; i++)
834  {
835  for (size_t j = 0; j < ncA; j++)
836  {
837  O(i, j) = A(i, j) * B[i];
838  }
839  }
840  }
841  else
842  {
843  for (size_t i = 0; i < nrA; i++)
844  {
845  for (size_t j = 0; j < ncA; j++)
846  {
847  O(i, j) = A(i, j) * B[j];
848  }
849  }
850  }
851  }
852 
860  template<class Scalar>
861  static void kroneckerSum(const Matrix<Scalar>& A, const Matrix<Scalar>& B, Matrix<Scalar>& O)
862  {
863  size_t ncA = A.getNumberOfColumns();
864  size_t nrA = A.getNumberOfRows();
865  size_t nrB = B.getNumberOfRows();
866  size_t ncB = B.getNumberOfColumns();
867  O.resize(nrA + nrB, ncA + ncB);
868  for (size_t ia = 0; ia < nrA; ia++)
869  {
870  for (size_t ja = 0; ja < ncA; ja++)
871  {
872  O(ia, ja) = A(ia, ja);
873  }
874  }
875  for (size_t ib = 0; ib < nrB; ib++)
876  {
877  for (size_t jb = 0; jb < nrB; jb++)
878  {
879  O(nrA + ib, ncA + jb) = B(ib, jb);
880  }
881  }
882  }
883 
890  template<class Scalar>
891  static void kroneckerSum(const std::vector< Matrix<Scalar>*>& vA, Matrix<Scalar>& O)
892  {
893  size_t nr = 0;
894  size_t nc = 0;
895  for (size_t k = 0; k < vA.size(); k++)
896  {
897  nr += vA[k]->getNumberOfRows();
898  nc += vA[k]->getNumberOfColumns();
899  }
900  O.resize(nr, nc);
901  size_t rk = 0; // Row counter
902  size_t ck = 0; // Col counter
903  for (size_t k = 0; k < vA.size(); k++)
904  {
905  const Matrix<Scalar>* Ak = vA[k];
906  for (size_t i = 0; i < Ak->getNumberOfRows(); i++)
907  {
908  for (size_t j = 0; j < Ak->getNumberOfColumns(); j++)
909  {
910  O(rk + i, ck + j) = (*Ak)(i, j);
911  }
912  }
913  rk += Ak->getNumberOfRows();
914  ck += Ak->getNumberOfColumns();
915  }
916  }
917 
924  template<class Scalar>
925  static void toVVdouble(const Matrix<Scalar>& M, std::vector< std::vector<Scalar> >& vO)
926  {
927  size_t n = M.getNumberOfRows();
928  size_t m = M.getNumberOfColumns();
929  vO.resize(n);
930  for (size_t i = 0; i < n; i++)
931  {
932  vO[i].resize(m);
933  for (size_t j = 0; j < m; j++)
934  {
935  vO[i][j] = M(i, j);
936  }
937  }
938  }
939 
945  template<class Scalar>
946  static Scalar sumElements(const Matrix<Scalar>& M)
947  {
948  Scalar sum = 0;
949  for (size_t i = 0; i < M.getNumberOfRows(); i++)
950  {
951  for (size_t j = 0; j < M.getNumberOfColumns(); j++)
952  {
953  sum += M(i, j);
954  }
955  }
956  return sum;
957  }
958 
959 
960 
975  template<class Scalar>
976  static Scalar lap(Matrix<Scalar>& assignCost,
977  std::vector<int> &rowSol,
978  std::vector<int> &colSol,
979  std::vector<Scalar> &u,
980  std::vector<Scalar> &v) throw (Exception)
981  {
982  size_t dim = assignCost.getNumberOfRows();
983  if (assignCost.getNumberOfColumns() != dim)
984  throw Exception("MatrixTools::lap. Cost matrix should be scare.");
985 
986  bool unassignedFound;
987  size_t i, iMin;
988  size_t numFree = 0, previousNumFree, f, k, freeRow;
989  int i0;
990  std::vector<size_t> free(dim); // list of unassigned rows.
991  std::vector<size_t> pred(dim); // row-predecessor of column in augmenting/alternating path.
992  size_t j, j1, j2, endOfPath, last, low, up;
993  std::vector<size_t> colList(dim); // list of columns to be scanned in various ways.
994  std::vector<short int> matches(dim, 0); // counts how many times a row could be assigned.
995  Scalar min;
996  Scalar h;
997  size_t uMin, uSubMin;
998  Scalar v2;
999  std::vector<Scalar> d(dim); // 'cost-distance' in augmenting path calculation.
1000 
1001  // Column reduction
1002  for (j = dim; j > 0; j--) // reverse order gives better results.
1003  {
1004  // find minimum cost over rows.
1005  min = assignCost(0, j - 1);
1006  iMin = 0;
1007  for (i = 1; i < dim; ++i) {
1008  if (assignCost(i, j - 1) < min)
1009  {
1010  min = assignCost(i, j - 1);
1011  iMin = i;
1012  }
1013  }
1014  v[j - 1] = min;
1015 
1016  if (++matches[iMin] == 1)
1017  {
1018  // init assignment if minimum row assigned for first time.
1019  rowSol[iMin] = static_cast<int>(j - 1);
1020  colSol[j - 1] = static_cast<int>(iMin);
1021  }
1022  else
1023  colSol[j - 1] = -1; // row already assigned, column not assigned.
1024  }
1025 
1026  // Reduction tranfer
1027  for (i = 0; i < dim; i++) {
1028  if (matches[i] == 0) // fill list of unassigned 'free' rows.
1029  free[numFree++] = i;
1030  else {
1031  if (matches[i] == 1) // transfer reduction from rows that are assigned once.
1032  {
1033  j1 = static_cast<size_t>(rowSol[i]); //rowSol[i] is >= 0 here
1034  min = -log(0);
1035  for (j = 0; j < dim; j++)
1036  if (j != j1)
1037  if (assignCost(i, j - 1) - v[j] < min)
1038  min = assignCost(i, j - 1) - v[j - 1];
1039  v[j1] = v[j1] - min;
1040  }
1041  }
1042  }
1043 
1044  // Augmenting row reduction
1045  short loopcnt = 0; // do-loop to be done twice.
1046  do
1047  {
1048  loopcnt++;
1049 
1050  // scan all free rows.
1051  // in some cases, a free row may be replaced with another one to be scanned next.
1052  k = 0;
1053  previousNumFree = numFree;
1054  numFree = 0; // start list of rows still free after augmenting row reduction.
1055  while (k < previousNumFree)
1056  {
1057  i = free[k];
1058  k++;
1059 
1060  // find minimum and second minimum reduced cost over columns.
1061  uMin = assignCost(i, 0) - v[0];
1062  j1 = 0;
1063  uSubMin = static_cast<size_t>(-log(0));
1064  for (j = 1; j < dim; j++)
1065  {
1066  h = assignCost(i, j) - v[j];
1067  if (h < uSubMin) {
1068  if (h >= uMin)
1069  {
1070  uSubMin = h;
1071  j2 = j;
1072  }
1073  else
1074  {
1075  uSubMin = uMin;
1076  uMin = h;
1077  j2 = j1;
1078  j1 = j;
1079  }
1080  }
1081  }
1082 
1083  i0 = colSol[j1];
1084  if (uMin < uSubMin) {
1085  // change the reduction of the minimum column to increase the minimum
1086  // reduced cost in the row to the subminimum.
1087  v[j1] = v[j1] - (uSubMin - uMin);
1088  } else { // minimum and subminimum equal.
1089  if (i0 >= 0) // minimum column j1 is assigned.
1090  {
1091  // swap columns j1 and j2, as j2 may be unassigned.
1092  j1 = j2;
1093  i0 = colSol[j2];
1094  }
1095  }
1096 
1097  // (re-)assign i to j1, possibly de-assigning an i0.
1098  rowSol[i] = static_cast<int>(j1);
1099  colSol[j1] = static_cast<int>(i);
1100 
1101  if (i0 >= 0) { // minimum column j1 assigned earlier.
1102  if (uMin < uSubMin) {
1103  // put in current k, and go back to that k.
1104  // continue augmenting path i - j1 with i0.
1105  free[--k] = static_cast<size_t>(i0);
1106  } else {
1107  // no further augmenting reduction possible.
1108  // store i0 in list of free rows for next phase.
1109  free[numFree++] = static_cast<size_t>(i0);
1110  }
1111  }
1112  }
1113  }
1114  while (loopcnt < 2); // repeat once.
1115 
1116  // Augment solution for each free row.
1117  for (f = 0; f < numFree; f++)
1118  {
1119  freeRow = free[f]; // start row of augmenting path.
1120 
1121  // Dijkstra shortest path algorithm.
1122  // runs until unassigned column added to shortest path tree.
1123  for (j = 0; j < dim; j++)
1124  {
1125  d[j] = assignCost(freeRow, j) - v[j];
1126  pred[j] = freeRow;
1127  colList[j] = j; // init column list.
1128  }
1129 
1130  low = 0; // columns in 0..low-1 are ready, now none.
1131  up = 0; // columns in low..up-1 are to be scanned for current minimum, now none.
1132  // columns in up..dim-1 are to be considered later to find new minimum,
1133  // at this stage the list simply contains all columns
1134  unassignedFound = false;
1135  do
1136  {
1137  if (up == low) // no more columns to be scanned for current minimum.
1138  {
1139  last = low - 1;
1140 
1141  // scan columns for up..dim-1 to find all indices for which new minimum occurs.
1142  // store these indices between low..up-1 (increasing up).
1143  min = d[colList[up++]];
1144  for (k = up; k < dim; k++)
1145  {
1146  j = colList[k];
1147  h = d[j];
1148  if (h <= min)
1149  {
1150  if (h < min) // new minimum.
1151  {
1152  up = low; // restart list at index low.
1153  min = h;
1154  }
1155  // new index with same minimum, put on undex up, and extend list.
1156  colList[k] = colList[up];
1157  colList[up++] = j;
1158  }
1159  }
1160 
1161  // check if any of the minimum columns happens to be unassigned.
1162  // if so, we have an augmenting path right away.
1163  for (k = low; k < up; k++) {
1164  if (colSol[colList[k]] < 0)
1165  {
1166  endOfPath = colList[k];
1167  unassignedFound = true;
1168  break;
1169  }
1170  }
1171  }
1172 
1173  if (!unassignedFound)
1174  {
1175  // update 'distances' between freerow and all unscanned columns, via next scanned column.
1176  j1 = colList[low];
1177  low++;
1178  i = static_cast<size_t>(colSol[j1]);
1179  h = assignCost(i, j1) - v[j1] - min;
1180 
1181  for (k = up; k < dim; k++)
1182  {
1183  j = colList[k];
1184  v2 = assignCost(i, j) - v[j] - h;
1185  if (v2 < d[j])
1186  {
1187  pred[j] = i;
1188  if (v2 == min) { // new column found at same minimum value
1189  if (colSol[j] < 0)
1190  {
1191  // if unassigned, shortest augmenting path is complete.
1192  endOfPath = j;
1193  unassignedFound = true;
1194  break;
1195  }
1196  // else add to list to be scanned right away.
1197  else
1198  {
1199  colList[k] = colList[up];
1200  colList[up++] = j;
1201  }
1202  }
1203  d[j] = v2;
1204  }
1205  }
1206  }
1207  }
1208  while (!unassignedFound);
1209 
1210  // update column prices.
1211  for (k = 0; k <= last; k++)
1212  {
1213  j1 = colList[k];
1214  v[j1] = v[j1] + d[j1] - min;
1215  }
1216 
1217  // reset row and column assignments along the alternating path.
1218  do
1219  {
1220  i = pred[endOfPath];
1221  colSol[endOfPath] = static_cast<int>(i);
1222  j1 = endOfPath;
1223  endOfPath = static_cast<size_t>(rowSol[i]);
1224  rowSol[i] = static_cast<int>(j1);
1225  }
1226  while (i != freeRow);
1227  }
1228 
1229  // calculate optimal cost.
1230  Scalar lapCost = 0;
1231  for (i = 0; i < dim; i++)
1232  {
1233  j = static_cast<size_t>(rowSol[i]);
1234  u[i] = assignCost(i, j) - v[j];
1235  lapCost = lapCost + assignCost(i, j);
1236  }
1237 
1238  return lapCost;
1239  }
1240 
1241  };
1242 
1243 /* DEPRECATED
1244  namespace MatrixOperators {
1245 
1246  MatrixB operator*(const MatrixA & A, const MatrixB & B) throw (DimensionException)
1247  {
1248  return MatrixTools::mult<MatrixA, MatrixB>(A, B);
1249  }
1250 
1251  template<class MatrixA, class MatrixB>
1252  MatrixA operator+(const MatrixA & A, const MatrixB & B) throw (DimensionException)
1253  {
1254  MatrixA C = A;
1255  MatrixTools::add<MatrixA, MatrixB>(C, B);
1256  return C;
1257  }
1258 
1259  template<class MatrixA, class MatrixB>
1260  MatrixA operator+=(MatrixA & A, const MatrixB & B) throw (DimensionException)
1261  {
1262  MatrixTools::add<MatrixA, MatrixB>(A, B);
1263  return A;
1264  }
1265 
1266  template<class Matrix>
1267  Matrix operator-(const Matrix A)
1268  {
1269  Matrix B(A.getNumberOfRows(), A.getNumberOfColumns());
1270  for(size_t i = 0; i < B.getNumberOfRows(); i++) {
1271  for(size_t j = 0; j < B.getNumberOfColumns(); j++) {
1272  B(i, j) = -A(i, j);
1273  }
1274  }
1275  return B;
1276  }
1277 
1278  template<class MatrixA, class MatrixB>
1279  MatrixA operator-(const MatrixA & A, const MatrixB & B) throw (DimensionException)
1280  {
1281  // size_t ncA = A.getNumberOfColumns();
1282  // size_t nrA = A.getNumberOfRows();
1283  // size_t nrB = B.getNumberOfRows();
1284  // size_t ncB = B.getNumberOfColumns();
1285  // if(ncA != ncB) throw DimensionException("MatrixTools::operator-(). A and B must have the same number of colums.", ncB, ncA);
1286  // if(nrA != nrB) throw DimensionException("MatrixTools::operator-(). A and B must have the same number of rows.", nrB, nrA);
1287  // MatrixB C(A.getNumberOfRows(), A.getNumberOfColumns());
1288  // for(size_t i = 0; i < A.getNumberOfRows(); i++) {
1289  // for(size_t j = 0; j < A.getNumberOfColumns(); j++) {
1290  // C(i, j) = A(i, j) - B(i, j);
1291  // }
1292  // }
1293  // return C;
1294  MatrixA C = A;
1295  MatrixTools::add<MatrixA, MatrixB>(C, -B);
1296  return C;
1297  }
1298 
1299  template<class MatrixA, class MatrixB>
1300  MatrixA operator-=(MatrixA & A, const MatrixB & B) throw (DimensionException)
1301  {
1302  MatrixTools::add<MatrixA, MatrixB>(A, -B);
1303  return A;
1304  }
1305 
1306  };
1307 */
1308 } // end of namespace bpp.
1309 
1310 #endif // _MATRIXTOOLS_H_
The matrix template interface.
Definition: Matrix.h:58
static void add(MatrixA &A, const MatrixB &B)
Add matrix B to matrix A.
Definition: MatrixTools.h:300
static void printForR(const Matrix &m, const std::string &variableName="x", std::ostream &out=std::cout)
Print a matrix to a stream, so that it is read by R.
Definition: MatrixTools.h:625
The base class exception for zero division error.
Definition: Exceptions.h:149
static std::vector< T > pow(const std::vector< T > &v1, T &b)
Definition: VectorTools.h:824
static void kroneckerSum(const std::vector< Matrix< Scalar > *> &vA, Matrix< Scalar > &O)
Compute the Kronecker sum of n row matrices.
Definition: MatrixTools.h:891
static bool isSquare(const Matrix &A)
Definition: MatrixTools.h:666
static void diag(const std::vector< Scalar > &D, Matrix< Scalar > &O)
Definition: MatrixTools.h:106
This class allows to perform a correspondence analysis.
virtual void resize(size_t nRows, size_t nCols)=0
Resize the matrix.
static void Taylor(const Matrix &A, size_t p, std::vector< RowMatrix< Scalar > > &vO)
Compute a vector of the first powers of a given matrix.
Definition: MatrixTools.h:438
static void kroneckerSum(const Matrix< Scalar > &A, const Matrix< Scalar > &B, Matrix< Scalar > &O)
Compute the Kronecker sum of two row matrices.
Definition: MatrixTools.h:861
static Real min(const Matrix< Real > &m)
Definition: MatrixTools.h:547
static void print(const Matrix &m, std::ostream &out=std::cout)
Print a matrix to a stream.
Definition: MatrixTools.h:573
static std::vector< size_t > whichMin(const Matrix &m)
Definition: MatrixTools.h:490
static void hadamardMult(const Matrix< Scalar > &A, const std::vector< Scalar > &B, Matrix< Scalar > &O, bool row=true)
Compute the "Hadamard" product of a row matrix and a vector containing weights, according to rows or ...
Definition: MatrixTools.h:823
static void copy(const MatrixA &A, MatrixO &O)
Copy operation. This function supplies the lack of inheritence of the assigment operator :D ...
Definition: MatrixTools.h:71
static void mult(const Matrix< Scalar > &A, const std::vector< Scalar > &D, const std::vector< Scalar > &U, const std::vector< Scalar > &L, const Matrix< Scalar > &B, Matrix< Scalar > &O)
Compute A . (U+D+L) . B where D is a diagonal matrix, U (resp. L) is a matrix in which the only non-z...
Definition: MatrixTools.h:263
static void exp(const Matrix< Scalar > &A, Matrix< Scalar > &O)
Perform matrix exponentiation using diagonalization.
Definition: MatrixTools.h:416
Matrix storage by row.
Definition: Matrix.h:129
static void print(const Matrix &m, bpp::OutputStream &out, char pIn='(', char pOut=')')
Print a matrix to a stream.
Definition: MatrixTools.h:598
static void diag(const Matrix< Scalar > &M, std::vector< Scalar > &O)
Definition: MatrixTools.h:137
virtual size_t getNumberOfColumns() const =0
const RowMatrix< Real > & getV() const
Return the eigenvector matrix.
Definition: EigenValue.h:1174
static Scalar inv(const Matrix< Scalar > &A, Matrix< Scalar > &O)
Definition: MatrixTools.h:675
static void getId(size_t n, Matrix &O)
Get a identity matrix of a given size.
Definition: MatrixTools.h:90
static Scalar sumElements(const Matrix< Scalar > &M)
Sum all elements in M.
Definition: MatrixTools.h:946
static Real max(const Matrix< Real > &m)
Definition: MatrixTools.h:521
static std::vector< size_t > whichMax(const Matrix &m)
Definition: MatrixTools.h:458
static void mult(const Matrix< Scalar > &A, const Matrix< Scalar > &B, Matrix< Scalar > &O)
Definition: MatrixTools.h:190
static void pow(const Matrix &A, size_t p, Matrix &O)
Compute the power of a given matrix.
Definition: MatrixTools.h:355
Computes eigenvalues and eigenvectors of a real (non-complex) matrix.
Definition: EigenValue.h:105
OutputStream interface.
Definition: OutputStream.h:64
LU Decomposition.
static void pow(const Matrix< Scalar > &A, double p, Matrix< Scalar > &O)
Compute the power of a given matrix, using eigen value decomposition.
Definition: MatrixTools.h:394
static std::vector< double > exp(const std::vector< T > &v1)
Definition: VectorTools.h:776
static void transpose(const MatrixA &A, MatrixO &O)
Definition: MatrixTools.h:706
Exception base class.
Definition: Exceptions.h:57
static void scale(Matrix &A, Scalar a, Scalar b=0)
Multiply all elements of a matrix by a given value, and add a constant.
Definition: MatrixTools.h:173
static void diag(const Scalar x, size_t n, Matrix< Scalar > &O)
Definition: MatrixTools.h:122
static void hadamardMult(const Matrix< Scalar > &A, const Matrix< Scalar > &B, Matrix< Scalar > &O)
Compute the Hadamard product of two row matrices with same dimensions.
Definition: MatrixTools.h:796
const std::vector< Real > & getRealEigenValues() const
Return the real parts of the eigenvalues.
Definition: EigenValue.h:1181
static void print(const std::vector< Real > &v, std::ostream &out=std::cout)
Print a vector to a stream.
Definition: MatrixTools.h:649
static void mult(const Matrix< Scalar > &A, const std::vector< Scalar > &D, const Matrix< Scalar > &B, Matrix< Scalar > &O)
Compute A . D . B where D is a diagonal matrix in O(n^3).
Definition: MatrixTools.h:224
static void toVVdouble(const Matrix< Scalar > &M, std::vector< std::vector< Scalar > > &vO)
Convert to a vector of vector.
Definition: MatrixTools.h:925
virtual size_t getNumberOfRows() const =0
static void kroneckerMult(const Matrix< Scalar > &A, const Matrix< Scalar > &B, Matrix< Scalar > &O)
Compute the Kronecker product of two row matrices.
Definition: MatrixTools.h:765
Exception thrown when a dimension problem occured.
static void fill(Matrix &M, Scalar x)
Set all elements in M to value x.
Definition: MatrixTools.h:152
static Scalar lap(Matrix< Scalar > &assignCost, std::vector< int > &rowSol, std::vector< int > &colSol, std::vector< Scalar > &u, std::vector< Scalar > &v)
Linear Assignment Problem.
Definition: MatrixTools.h:976
static void add(MatrixA &A, Scalar &x, const MatrixB &B)
Add matrix x.B to matrix A.
Definition: MatrixTools.h:326
static void covar(const Matrix< Scalar > &A, Matrix< Scalar > &O)
Compute the variance-covariance matrix of an input matrix.
Definition: MatrixTools.h:731
Functions dealing with matrices.
Definition: MatrixTools.h:57
static double det(const Matrix< Scalar > &A)
Get determinant of a square matrix.
Definition: MatrixTools.h:694