40 #ifndef _MATRIXTOOLS_H_ 41 #define _MATRIXTOOLS_H_ 43 #include "../VectorTools.h" 47 #include "../../Io/OutputStream.h" 70 template<
class MatrixA,
class MatrixO>
71 static void copy(
const MatrixA& A, MatrixO& O)
73 O.resize(A.getNumberOfRows(), A.getNumberOfColumns());
74 for (
size_t i = 0; i < A.getNumberOfRows(); i++)
76 for (
size_t j = 0; j < A.getNumberOfColumns(); j++)
89 template<
class Matrix>
93 for (
size_t i = 0; i < n; i++)
95 for (
size_t j = 0; j < n; j++) {
96 O(i, j) = (i == j) ? 1 : 0;
105 template<
class Scalar>
110 for (
size_t i = 0; i < n; i++)
112 for (
size_t j = 0; j < n; j++) { O(i, j) = (i == j) ? D[i] : 0;}
121 template<
class Scalar>
125 for (
size_t i = 0; i < n; i++)
127 for (
size_t j = 0; j < n; j++) { O(i, j) = (i == j) ? x : 0;}
136 template<
class Scalar>
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);
143 for (
size_t i = 0; i < nc; i++) { O[i] = M(i, i);}
151 template<
class Matrix,
class Scalar>
172 template<
class Matrix,
class Scalar>
179 A(i, j) = a * A(i, j) + b;
189 template<
class Scalar>
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);
198 for (
size_t i = 0; i < nrA; i++)
200 for (
size_t j = 0; j < ncB; j++)
203 for (
size_t k = 0; k < ncA; k++)
205 O(i, j) += A(i, k) * B(k, j);
223 template<
class Scalar>
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);
233 for (
size_t i = 0; i < nrA; i++)
235 for (
size_t j = 0; j < ncB; j++)
238 for (
size_t k = 0; k < ncA; k++)
240 O(i, j) += A(i, k) * B(k, j) * D[k];
262 template<
class Scalar>
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);
274 for (
size_t i = 0; i < nrA; i++)
276 for (
size_t j = 0; j < ncB; j++)
278 O(i, j) = A(i, 0) * D[0] * B(0, j);
280 O(i, j) += A(i,0) * U[0] * B(1,j);
281 for (
size_t k = 1; k < ncA-1; k++)
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));
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);
299 template<
class MatrixA,
class MatrixB>
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++)
310 for (
size_t j = 0; j < A.getNumberOfColumns(); j++)
325 template<
class MatrixA,
class MatrixB,
class Scalar>
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++)
336 for (
size_t j = 0; j < A.getNumberOfColumns(); j++)
338 A(i, j) += x*B(i, j);
354 template<
class Matrix>
357 size_t n = A.getNumberOfRows();
358 if (n != A.getNumberOfColumns())
throw DimensionException(
"MatrixTools::pow(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
393 template<
class Scalar>
396 size_t n = A.getNumberOfRows();
397 if (n != A.getNumberOfColumns())
throw DimensionException(
"MatrixTools::pow(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
400 rightEV = eigen.
getV();
401 inv(rightEV, leftEV);
415 template<
class Scalar>
418 size_t n = A.getNumberOfRows();
419 if (n != A.getNumberOfColumns())
throw DimensionException(
"MatrixTools::exp(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
422 rightEV = eigen.
getV();
423 inv(rightEV, leftEV);
437 template<
class Matrix,
class Scalar>
440 size_t n = A.getNumberOfRows();
441 if (n != A.getNumberOfColumns())
442 throw DimensionException(
"MatrixTools::pow(). nrows != ncols.", A.getNumberOfColumns(), A.getNumberOfRows());
444 getId<Matrix>(n, vO[0]);
447 for (
size_t i = 1; i < p; i++)
449 mult(vO[i], A, vO[i+1]);
457 template<
class Matrix>
462 std::vector<size_t> pos(2);
465 double currentMax = log(0.);
466 for (
size_t i = 0; i < nrows; i++)
468 for (
size_t j = 0; j < ncols; j++)
470 double currentValue = m(i, j);
472 if (currentValue > currentMax)
476 currentMax = currentValue;
489 template<
class Matrix>
494 std::vector<size_t> pos(2);
497 double currentMin = -log(0.);
498 for (
size_t i = 0; i < nrows; i++)
500 for (
size_t j = 0; j < ncols; j++)
502 double currentValue = m(i, j);
503 if (currentValue < currentMin)
507 currentMin = currentValue;
525 Real currentMax = log(0.);
526 for (
size_t i = 0; i < nrows; i++)
528 for (
size_t j = 0; j < ncols; j++)
530 Real currentValue = m(i, j);
532 if (currentValue > currentMax)
534 currentMax = currentValue;
551 Real currentMin = -log(0.);
552 for (
size_t i = 0; i < nrows; i++)
554 for (
size_t j = 0; j < ncols; j++)
556 Real currentValue = m(i, j);
557 if (currentValue < currentMin)
559 currentMin = currentValue;
572 template<
class Matrix>
576 out <<
"[" << std::endl;
582 out << m(i, j) <<
", ";
586 out <<
"]" << std::endl;
597 template<
class Matrix>
610 out << m(i, j) <<
", ";
624 template<
class Matrix>
625 static void printForR(
const Matrix& m,
const std::string& variableName =
"x", std::ostream& out = std::cout)
628 out << variableName <<
"<-matrix(c(";
638 out <<
"), nrow=" << m.
getNumberOfRows() <<
", byrow=TRUE)" << std::endl;
649 static void print(
const std::vector<Real>& v, std::ostream& out = std::cout)
651 out << v.size() << std::endl;
653 for (
size_t i = 0; i < v.size() - 1; i++)
657 if (v.size() > 0) out << v[v.size() - 1];
658 out <<
"]" << std::endl;
665 template<
class Matrix>
674 template<
class Scalar>
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);
693 template<
class Scalar>
696 if (!
isSquare(A))
throw DimensionException(
"MatrixTools::det(). Matrix A is not a square matrix.", A.getNumberOfRows(), A.getNumberOfColumns());
705 template<
class MatrixA,
class MatrixO>
708 O.resize(A.getNumberOfColumns(), A.getNumberOfRows());
709 for (
size_t i = 0; i < A.getNumberOfColumns(); i++)
711 for (
size_t j = 0; j < A.getNumberOfRows(); j++)
730 template<
class Scalar>
739 scale(O, 1. / static_cast<double>(n));
741 for (
size_t i = 0; i < r; i++)
743 for (
size_t j = 0; j < n; j++)
745 mean(i, 0) += A(i, j);
747 mean(i, 0) /=
static_cast<double>(n);
752 mult(mean, tMean, meanMat);
764 template<
class Scalar>
771 O.
resize(nrA * nrB, ncA * ncB);
772 for (
size_t ia = 0; ia < nrA; ia++)
774 for (
size_t ja = 0; ja < ncA; ja++)
776 Scalar aij = A(ia, ja);
777 for (
size_t ib = 0; ib < nrB; ib++)
779 for (
size_t jb = 0; jb < ncB; jb++)
781 O(ia * nrB + ib, ja * ncB + jb) = aij * B(ib, jb);
795 template<
class Scalar>
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);
805 for (
size_t i = 0; i < nrA; i++)
807 for (
size_t j = 0; j < ncA; j++)
809 O(i, j) = A(i, j) * B(i, j);
822 template<
class Scalar>
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);
833 for (
size_t i = 0; i < nrA; i++)
835 for (
size_t j = 0; j < ncA; j++)
837 O(i, j) = A(i, j) * B[i];
843 for (
size_t i = 0; i < nrA; i++)
845 for (
size_t j = 0; j < ncA; j++)
847 O(i, j) = A(i, j) * B[j];
860 template<
class Scalar>
867 O.
resize(nrA + nrB, ncA + ncB);
868 for (
size_t ia = 0; ia < nrA; ia++)
870 for (
size_t ja = 0; ja < ncA; ja++)
872 O(ia, ja) = A(ia, ja);
875 for (
size_t ib = 0; ib < nrB; ib++)
877 for (
size_t jb = 0; jb < nrB; jb++)
879 O(nrA + ib, ncA + jb) = B(ib, jb);
890 template<
class Scalar>
895 for (
size_t k = 0; k < vA.size(); k++)
897 nr += vA[k]->getNumberOfRows();
898 nc += vA[k]->getNumberOfColumns();
903 for (
size_t k = 0; k < vA.size(); k++)
910 O(rk + i, ck + j) = (*Ak)(i, j);
924 template<
class Scalar>
930 for (
size_t i = 0; i < n; i++)
933 for (
size_t j = 0; j < m; j++)
945 template<
class Scalar>
975 template<
class Scalar>
977 std::vector<int> &rowSol,
978 std::vector<int> &colSol,
979 std::vector<Scalar> &u,
980 std::vector<Scalar> &v)
throw (
Exception)
982 size_t dim = assignCost.getNumberOfRows();
983 if (assignCost.getNumberOfColumns() != dim)
984 throw Exception(
"MatrixTools::lap. Cost matrix should be scare.");
986 bool unassignedFound;
988 size_t numFree = 0, previousNumFree, f, k, freeRow;
990 std::vector<size_t> free(dim);
991 std::vector<size_t> pred(dim);
992 size_t j, j1, j2, endOfPath, last, low, up;
993 std::vector<size_t> colList(dim);
994 std::vector<short int> matches(dim, 0);
997 size_t uMin, uSubMin;
999 std::vector<Scalar> d(dim);
1002 for (j = dim; j > 0; j--)
1005 min = assignCost(0, j - 1);
1007 for (i = 1; i < dim; ++i) {
1008 if (assignCost(i, j - 1) <
min)
1010 min = assignCost(i, j - 1);
1016 if (++matches[iMin] == 1)
1019 rowSol[iMin] =
static_cast<int>(j - 1);
1020 colSol[j - 1] =
static_cast<int>(iMin);
1027 for (i = 0; i < dim; i++) {
1028 if (matches[i] == 0)
1029 free[numFree++] = i;
1031 if (matches[i] == 1)
1033 j1 =
static_cast<size_t>(rowSol[i]);
1035 for (j = 0; j < dim; j++)
1037 if (assignCost(i, j - 1) - v[j] <
min)
1038 min = assignCost(i, j - 1) - v[j - 1];
1039 v[j1] = v[j1] -
min;
1053 previousNumFree = numFree;
1055 while (k < previousNumFree)
1061 uMin = assignCost(i, 0) - v[0];
1063 uSubMin =
static_cast<size_t>(-log(0));
1064 for (j = 1; j < dim; j++)
1066 h = assignCost(i, j) - v[j];
1084 if (uMin < uSubMin) {
1087 v[j1] = v[j1] - (uSubMin - uMin);
1098 rowSol[i] =
static_cast<int>(j1);
1099 colSol[j1] =
static_cast<int>(i);
1102 if (uMin < uSubMin) {
1105 free[--k] =
static_cast<size_t>(i0);
1109 free[numFree++] =
static_cast<size_t>(i0);
1114 while (loopcnt < 2);
1117 for (f = 0; f < numFree; f++)
1123 for (j = 0; j < dim; j++)
1125 d[j] = assignCost(freeRow, j) - v[j];
1134 unassignedFound =
false;
1143 min = d[colList[up++]];
1144 for (k = up; k < dim; k++)
1156 colList[k] = colList[up];
1163 for (k = low; k < up; k++) {
1164 if (colSol[colList[k]] < 0)
1166 endOfPath = colList[k];
1167 unassignedFound =
true;
1173 if (!unassignedFound)
1178 i =
static_cast<size_t>(colSol[j1]);
1179 h = assignCost(i, j1) - v[j1] -
min;
1181 for (k = up; k < dim; k++)
1184 v2 = assignCost(i, j) - v[j] - h;
1193 unassignedFound =
true;
1199 colList[k] = colList[up];
1208 while (!unassignedFound);
1211 for (k = 0; k <= last; k++)
1214 v[j1] = v[j1] + d[j1] -
min;
1220 i = pred[endOfPath];
1221 colSol[endOfPath] =
static_cast<int>(i);
1223 endOfPath =
static_cast<size_t>(rowSol[i]);
1224 rowSol[i] =
static_cast<int>(j1);
1226 while (i != freeRow);
1231 for (i = 0; i < dim; i++)
1233 j =
static_cast<size_t>(rowSol[i]);
1234 u[i] = assignCost(i, j) - v[j];
1235 lapCost = lapCost + assignCost(i, j);
1310 #endif // _MATRIXTOOLS_H_ The matrix template interface.
The base class exception for zero division error.
This class allows to perform a correspondence analysis.
virtual void resize(size_t nRows, size_t nCols)=0
Resize the matrix.
virtual size_t getNumberOfColumns() const =0
const RowMatrix< Real > & getV() const
Return the eigenvector matrix.
Computes eigenvalues and eigenvectors of a real (non-complex) matrix.
const std::vector< Real > & getRealEigenValues() const
Return the real parts of the eigenvalues.
virtual size_t getNumberOfRows() const =0
Exception thrown when a dimension problem occured.