Index: trunk/tests/test_Matrix.cc =================================================================== --- trunk/tests/test_Matrix.cc (revision 690) +++ trunk/tests/test_Matrix.cc (revision 691) @@ -1,1391 +1,1439 @@ #include #include #include "UnitTest++.h" #include "test_utils.hh" #include "npstat/nm/Matrix.hh" #include "npstat/nm/SimpleFunctors.hh" using namespace npstat; using namespace std; namespace { struct InverterD { inline double operator()(const double& d) const {return 1.0/d;} }; void fillRandomInts(Matrix& m, const int maxint) { const unsigned nRows = m.nRows(); const unsigned nCols = m.nColumns(); for (unsigned row=0; row& m, const double location=0.0, + const double scale=1.0) + { + const unsigned nRows = m.nRows(); + const unsigned nCols = m.nColumns(); + for (unsigned row=0; row& m, const int maxint) { const unsigned nRows = m.nRows(); std::vector diag(nRows); for (unsigned i=0; i m1; Matrix m2; Matrix m3; Matrix m4(5, 5); Matrix m6(10, 10, 0); CHECK(m6.isDiagonal()); Matrix m7(10, 10, 1); CHECK(m7.isDiagonal()); Matrix > m8(10, 10, 0); Matrix > m9(10, 10, 1); Matrix > m10 = m9; CHECK_EQUAL(m9, m10); Matrix > m11(m9); CHECK_EQUAL(m9, m11); Matrix m12; CHECK(!m12.isDiagonal()); double data[] = {1, 2, 3}; m12.diagFill(data, sizeof(data)/sizeof(data[0])); CHECK(m12.isDiagonal()); } TEST(Matrix_equality) { Matrix m1(2,2,1); Matrix m2(2,2,0); CHECK_EQUAL(m1, m1); CHECK_EQUAL(false, m1 == m2); CHECK(m1 != m2); m2[0][0] = 1; m2[1][1] = 1; CHECK_EQUAL(m1, m2); } TEST(Matrix_Matrix_multiplication) { int a1[2][3] = {{1, 4, 2}, {2, 8, 13}}; int a2[3][4] = {{2, 5, -3, 8}, {8, 3, 16, 7}, {7, 4, 2, -3}}; int a3[2][4] = {{48, 25, 65, 30}, {159, 86, 148, 33}}; int a5[3][4] = {{6, 15, -9, 24}, {56, 21, 112, 49}, {77, 44, 22, -33}}; int a6[3][4] = {{6, 35, -33, 104}, {24, 21, 176, 91}, {21, 28, 22, -39}}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(3, 4, &a2[0][0]); Matrix m3(2, 4, &a3[0][0]); Matrix m5(3, 4, &a5[0][0]); Matrix m6(3, 4, &a6[0][0]); CHECK(!m3.isDiagonal()); CHECK_EQUAL(-39, m6.minValue()); CHECK_EQUAL(176, m6.maxValue()); CHECK_EQUAL(176, m6.maxAbsValue()); CHECK(std::make_pair(2U,3U) == m6.argmin()); CHECK(std::make_pair(1U,2U) == m6.argmax()); CHECK(std::make_pair(1U,2U) == m6.argmaxAbs()); CHECK_EQUAL(m1*m1.T(), m1.timesT()); CHECK_EQUAL(m2*m2.T(), m2.timesT()); CHECK_EQUAL(m3*m3.T(), m3.timesT()); CHECK_EQUAL(m5*m5.T(), m5.timesT()); CHECK_EQUAL(m1.T()*m1, m1.TtimesThis()); CHECK_EQUAL(m2.T()*m2, m2.TtimesThis()); CHECK_EQUAL(m3.T()*m3, m3.TtimesThis()); CHECK_EQUAL(m5.T()*m5, m5.TtimesThis()); Matrix diag33(3, 3); Matrix diag44(4, 4); int idiag[4] = {3, 7, 11, 13}; diag33.diagFill(idiag, 3); diag44.diagFill(idiag, 4); CHECK(diag33.isDiagonal()); CHECK_EQUAL(diag33*diag33.T(), diag33.timesT()); CHECK_EQUAL(diag44*diag44.T(), diag44.timesT()); CHECK_EQUAL(diag33.T()*diag33, diag33.TtimesThis()); CHECK_EQUAL(diag44.T()*diag44, diag44.TtimesThis()); const Matrix& m6t = m6.T(); const Matrix& trtest = m5*m6t; CHECK_EQUAL(trtest.tr(), m5.productTr(m6t)); CHECK_EQUAL(trtest.tr(), m6t.productTr(m5)); CHECK_EQUAL(m3, m1*m2); CHECK_EQUAL(m5, diag33*m2); CHECK_EQUAL(m6, m2*diag44); Matrix temp; m1.times(m2, &temp); CHECK_EQUAL(m3, temp); diag33.times(m2, &temp); CHECK_EQUAL(m5, temp); m2.times(diag44, &temp); CHECK_EQUAL(m6, temp); int a4[3][3] = {{1, 4, 2}, {2, 8, 13}, {27, 33, 45}}; Matrix m4(3, 3, &a4[0][0]); CHECK(!m4.isDiagonal()); Matrix unit(3, 3, 1); CHECK_EQUAL(m4, m4*unit); CHECK_EQUAL(m4, unit*m4); } TEST(Matrix_Ttimes_1) { const unsigned nCycles = 1000; for (unsigned ic=0; ic m1(n, n); fillRandomDiag(m1, 100); Matrix m2(n, k); fillRandomInts(m2, 100); CHECK_EQUAL(m2.T()*m1, m2.Ttimes(m1)); CHECK_EQUAL(m1.T()*m2, m1.Ttimes(m2)); } } TEST(Matrix_Ttimes_2) { const unsigned nCycles = 1000; for (unsigned ic=0; ic m1(n, m); fillRandomInts(m1, 100); Matrix m2(n, k); fillRandomInts(m2, 100); CHECK_EQUAL(m1.T()*m2, m1.Ttimes(m2)); } } TEST(Matrix_timesT_1) { const unsigned nCycles = 1000; for (unsigned ic=0; ic m1(n, n); fillRandomDiag(m1, 100); Matrix m2(k, n); fillRandomInts(m2, 100); CHECK_EQUAL(m1*m2.T(), m1.timesT(m2)); CHECK_EQUAL(m2*m1.T(), m2.timesT(m1)); } } TEST(Matrix_timesT_2) { const unsigned nCycles = 1000; for (unsigned ic=0; ic m1(m, n); fillRandomInts(m1, 100); Matrix m2(k, n); fillRandomInts(m2, 100); CHECK_EQUAL(m1*m2.T(), m1.timesT(m2)); } } TEST(Matrix_directSum) { const int a1[2][3] = {{1, 4, 2}, {2, 8, 13}}; const Matrix m1(2, 3, &a1[0][0]); const int a2[4][1] = {{7}, {25}, {36}, {49}}; const Matrix m2(4, 1, &a2[0][0]); const int a3[6][4] = {{1, 4, 2, 0}, {2, 8, 13, 0}, {0, 0, 0, 7}, {0, 0, 0, 25}, {0, 0, 0, 36}, {0, 0, 0, 49}}; const Matrix m3(6, 4, &a3[0][0]); CHECK_EQUAL(m3, m1.directSum(m2)); const int a4[6][4] = {{ 7, 0, 0, 0}, {25, 0, 0, 0}, {36, 0, 0, 0}, {49, 0, 0, 0}, { 0, 1, 4, 2}, { 0, 2, 8, 13}}; const Matrix m4(6, 4, &a4[0][0]); CHECK_EQUAL(m4, m2.directSum(m1)); } TEST(Matrix_Vector_multiplication) { int a1[2][3] = {{1, 4, 2}, {2, 8, 13}}; unsigned a2[3] = {6, 3, 15}; int a3[2] = {48, 231}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(3, 1, &a2[0]); Matrix m3(2, 1, &a3[0]); CHECK(m3 == m1*m2); CHECK_EQUAL(m3, m1.timesVector(a2, 3)); for (unsigned i=0; i<2; ++i) CHECK_EQUAL(m3[i][0], m1.timesVector(i, a2, 3)); int buf[2]; m1.timesVector(a2, 3, buf, 2); for (unsigned i=0; i<2; ++i) CHECK_EQUAL(m3[i][0], buf[i]); } TEST(Matrix_row_multiplication) { int a1[2][3] = {{1, 4, 2}, {2, 8, 13}}; unsigned a2[2] = {7, 11}; int a3[3] = {29, 116, 157}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(1, 2, &a2[0]); Matrix m3(1, 3, &a3[0]); CHECK(m3 == m2*m1); CHECK_EQUAL(m3, m1.rowMultiply(a2, 2)); for (unsigned i=0; i<3; ++i) CHECK_EQUAL(m3[0][i], m1.rowMultiply(i, a2, 2)); int buf[3]; m1.rowMultiply(a2, 2, buf, 3); for (unsigned i=0; i<3; ++i) CHECK_EQUAL(m3[0][i], buf[i]); } TEST(Matrix_unary_ops) { int a1[2][3] = {{1, 2, 3}, {4, 5, 6}}; Matrix m1(2, 3, &a1[0][0]); CHECK_EQUAL(-m1, m1*(-1)); CHECK_EQUAL(+m1, m1); } TEST(Matrix_Const_ops) { int a1[2][3] = {{1, 2, 3}, {4, 5, 6}}; int a3[2][3] = {{5, 10, 15}, {20, 25, 30}}; Matrix m1(2, 3, &a1[0][0]); Matrix m3(2, 3, &a3[0][0]); Matrix m2 = m1*5; Matrix m4 = m1; m1 *= 5; CHECK_EQUAL(m1, m2); CHECK_EQUAL(m1, m3); m1 /= 5; CHECK_EQUAL(m1, m4); } TEST(Vector_Matrix_multiplication) { int a1[2][3] = {{1, 4, 2}, {2, 8, 13}}; int a2[2] = {3, 4}; int a3[3] = {11, 44, 58}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(1, 2, &a2[0]); Matrix m3(1, 3, &a3[0]); CHECK_EQUAL(m3, m2*m1); } TEST(Matrix_transpose) { int a1[2][3] = {{1, 2, 3}, {4, 5, 6}}; int a2[3][2] = {{1, 4}, {2, 5}, {3, 6}}; int a3[3][3] = {{1, 2, 3}, {4, 5, 6}, {9, 8, 7}}; int a4[3][3] = {{1, 4, 9}, {2, 5, 8}, {3, 6, 7}}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(3, 2, &a2[0][0]); Matrix m3(3, 3, &a3[0][0]); Matrix m4(3, 3, &a4[0][0]); CHECK_EQUAL(m2, m1.T()); CHECK_EQUAL(m4, m3.T()); m4.Tthis(); CHECK_EQUAL(m4, m3); } TEST(Matrix_transpose_2) { int a1[2][3] = {{1, 2, 3}, {4, 5, 6}}; int a2[3][2] = {{1, 4}, {2, 5}, {3, 6}}; int a3[3][3] = {{1, 2, 3}, {4, 5, 6}, {9, 8, 7}}; int a4[3][3] = {{1, 4, 9}, {2, 5, 8}, {3, 6, 7}}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(3, 2, &a2[0][0]); Matrix m3(3, 3, &a3[0][0]); Matrix m4(3, 3, &a4[0][0]); CHECK_EQUAL(m2, m1.T()); CHECK_EQUAL(m4, m3.T()); m4.Tthis(); CHECK_EQUAL(m4, m3); } TEST(Matrix_Tthis) { const unsigned nCycles = 100; for (unsigned ic=0; ic m2(n, k); fillRandomInts(m2, 100); const Matrix mt = m2.T(); m2.Tthis(); CHECK_EQUAL(m2, mt); } } TEST(Matrix_pow) { Matrix m(3, 3); for (unsigned ntry=0; ntry<10; ++ntry) { for (unsigned i=0; i<3; ++i) for (unsigned j=0; j<3; ++j) m[i][j] = static_cast(test_rng()*100 - 50); CHECK_EQUAL(m.pow(0), Matrix(3, 3, 1)); CHECK_EQUAL(m.pow(1), m); CHECK_EQUAL(m.pow(2), m*m); CHECK_EQUAL(m.pow(3), m*m*m); CHECK_EQUAL(m.pow(4), m*m*m*m); } } TEST(Matrix_subtraction) { int a1[2][3] = {{1, 4, 2}, {2, 8, 13}}; int a2[2][3] = {{7, 5, 4}, {8, 3, -5}}; int a3[2][3] = {{-6, -1, -2,}, {-6, 5, 18}}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(2, 3, &a2[0][0]); Matrix m3(2, 3, &a3[0][0]); CHECK_EQUAL(m3, m1 - m2); m1 -= m2; CHECK_EQUAL(m3, m1); } TEST(Matrix_addition) { int a1[2][3] = {{1, 4, 2}, {2, 8, 13}}; int a2[2][3] = {{7, 5, 4}, {8, 3, -5}}; int a3[2][3] = {{-6, -1, -2,}, {-6, 5, 18}}; Matrix m1(2, 3, &a1[0][0]); Matrix m2(2, 3, &a2[0][0]); Matrix m3(2, 3, &a3[0][0]); CHECK_EQUAL(m1, m3 + m2); CHECK_EQUAL(m1, m2 + m3); m3 += m2; CHECK_EQUAL(m1, m3); } TEST(Matrix_setRow) { double a1[3][3] = {{1, 4, 2}, {2, 8, 13}, {27, 33, 45}}; double a2[3][3] = {{1, 4, 2}, {2, 8, 13}, {-1, 6, 5}}; int arr[] = {-1, 6, 5}; Matrix m1(3, 3, &a1[0][0]); m1.setRow(2, &arr[0], 3); Matrix m2(3, 3, &a2[0][0]); CHECK_EQUAL(m1, m2); } TEST(Matrix_setColumn) { double a1[3][3] = {{1, 4, 2}, {2, 8, 13}, {27, 33, 45}}; double a2[3][3] = {{1, -1, 2}, {2, 6, 13}, {27, 5, 45}}; int arr[] = {-1, 6, 5}; Matrix m1(3, 3, &a1[0][0]); m1.setColumn(1, &arr[0], 3); Matrix m2(3, 3, &a2[0][0]); CHECK_EQUAL(m1, m2); } TEST(Matrix_symmetrize) { double a1[3][3] = {{1, 4, 2}, {2, 8, 13}, {27, 33, 45}}; double a2[3][3] = {{1, 3, 14.5}, {3, 8, 23}, {14.5, 23, 45}}; Matrix m1(3, 3, &a1[0][0]); Matrix m2(3, 3, &a2[0][0]); CHECK(!m1.isSymmetric()); CHECK(!m1.isAntiSymmetric()); CHECK(m1.symmetrize().isSymmetric()); CHECK(m1.antiSymmetrize().isAntiSymmetric()); CHECK_EQUAL(m1.symmetrize(), m2); CHECK_EQUAL((m1 + m1.T())/2.0, m1.symmetrize()); CHECK_EQUAL((m1 - m1.T())/2.0, m1.antiSymmetrize()); CHECK_EQUAL(m1, m1.symmetrize()+m1.antiSymmetrize()); } TEST(Matrix_bilinear_1) { double a1[3][3] = {{1, 4, 2}, {2, 8, 13}, {27, 33, 45}}; Matrix m1(3, 3, &a1[0][0]); double vec1[3] = {7, 11, 17}; double vec2[3] = {23, 29, 31}; CHECK_EQUAL(26537.0, m1.bilinear(vec1, 3)); CHECK_EQUAL(116535.0, m1.bilinear(vec2, 3)); } TEST(Matrix_bilinear_2) { const unsigned nCycles = 1000; for (unsigned ic=0; ic m1(n, n); fillRandomInts(m1, 100); Matrix m2(n, k); fillRandomInts(m2, 100); Matrix m3(k, k); fillRandomInts(m1, 100); CHECK_EQUAL(m2.T()*m1*m2, m1.bilinear(m2)); CHECK_EQUAL(m2*m3*m2.T(), m3.bilinearT(m2)); } } TEST(Matrix_bilinear_3) { const unsigned nCycles = 1000; for (unsigned ic=0; ic m1(n, n); fillRandomDiag(m1, 100); Matrix m2(n, k); fillRandomInts(m2, 100); Matrix m3(k, k); fillRandomDiag(m3, 100); CHECK_EQUAL(m2.T()*m1*m2, m1.bilinear(m2)); CHECK_EQUAL(m2*m3*m2.T(), m3.bilinearT(m2)); } } TEST(Matrix_bilinear_4) { const unsigned nCycles = 1000; for (unsigned ic=0; ic m1(n, n); fillRandomInts(m1, 100); Matrix m2(n, n); fillRandomDiag(m2, 100); CHECK_EQUAL(m2.T()*m1*m2, m1.bilinear(m2)); CHECK_EQUAL(m2*m1*m2.T(), m1.bilinearT(m2)); } } TEST(Matrix_bilinear_5) { const unsigned nCycles = 100; for (unsigned ic=0; ic m1(n, n); fillRandomDiag(m1, 100); Matrix m2(n, n); fillRandomDiag(m2, 100); CHECK_EQUAL(m2.T()*m1*m2, m1.bilinear(m2)); CHECK_EQUAL(m2*m1*m2.T(), m1.bilinearT(m2)); } } TEST(Matrix_trace_det) { double a1[3][3] = {{1, 4, 2}, {2, 8, 13}, {27, 33, 45}}; Matrix m1(3, 3, &a1[0][0]); CHECK_EQUAL(54.0, m1.sp()); CHECK_EQUAL(54.0, m1.tr()); CHECK_EQUAL(675.0, m1.det()); } TEST(Matrix_inverse_double) { double a1[3][3] = {{4, 2.4, 9.0}, {2.4, 9, 10.5}, {9.0, 10.5, 25}}; double a2[3][3] = {{2.89772727273, 0.871212121212, -1.40909090909}, {0.871212121212, 0.479797979798, -0.515151515152}, {-1.40909090909, -0.515151515152, 0.763636363636}}; Matrix m1(3, 3, &a1[0][0]); Matrix m2(3, 3, &a2[0][0]); Matrix unit(3, 3, 1); CHECK_CLOSE(0.0, (m1*m2 - unit).maxAbsValue(), 1.0e-10); CHECK_CLOSE(0.0, (m1.symPDInv() - m2).maxAbsValue(), 1.0e-10); CHECK_CLOSE(0.0, (m1.symPDEigenInv(1.e-12) - m2).maxAbsValue(), 1.0e-10); CHECK_CLOSE(0.0, (m1.symInv() - m2).maxAbsValue(), 1.0e-10); CHECK_CLOSE(0.0, (m1.inv() - m2).maxAbsValue(), 1.0e-10); CHECK_CLOSE(0.0, (m1.symFcn(InverterD()) - m2).maxAbsValue(), 1.0e-10); } TEST(Matrix_inverse_float) { float a1[3][3] = {{4, 2.4, 9.0}, {2.4, 9, 10.5}, {9.0, 10.5, 25}}; float a2[3][3] = {{2.89772727273, 0.871212121212, -1.40909090909}, {0.871212121212, 0.479797979798, -0.515151515152}, {-1.40909090909, -0.515151515152, 0.763636363636}}; Matrix m1(3, 3, &a1[0][0]); Matrix m2(3, 3, &a2[0][0]); Matrix unit(3, 3, 1); CHECK_CLOSE(0.0, (m1*m2 - unit).maxAbsValue(), 2.0e-6); CHECK_CLOSE(0.0, (m1.symPDInv() - m2).maxAbsValue(), 2.0e-6); CHECK_CLOSE(0.0, (m1.symPDEigenInv(1.e-7) - m2).maxAbsValue(), 2.0e-5); CHECK_CLOSE(0.0, (m1.symInv() - m2).maxAbsValue(), 4.0e-6); CHECK_CLOSE(0.0, (m1.inv() - m2).maxAbsValue(), 4.0e-6); } + TEST(Matrix_leftInv) + { + const double eps = 1.0e-9; + + for (unsigned ncols=2; ncols<6; ++ncols) + { + const Matrix unit(ncols, ncols, 1); + + for (unsigned nrows=ncols; nrows<6; ++nrows) + { + Matrix m(nrows, ncols); + fillRandomDoubles(m); + const Matrix& inv = m.leftInv(); + const Matrix& prod = inv*m; + CHECK_CLOSE(0.0, (prod - unit).maxAbsValue(), eps); + } + } + } + + TEST(Matrix_rightInv) + { + const double eps = 1.0e-9; + + for (unsigned nrows=2; nrows<6; ++nrows) + { + const Matrix unit(nrows, nrows, 1); + + for (unsigned ncols=nrows; ncols<6; ++ncols) + { + Matrix m(nrows, ncols); + fillRandomDoubles(m); + const Matrix& inv = m.rightInv(); + const Matrix& prod = m*inv; + CHECK_CLOSE(0.0, (prod - unit).maxAbsValue(), eps); + } + } + } + TEST(Matrix_covarToCorr) { double a1[2][2] = {{9.0, 10.5}, {10.5, 25.0}}; double a2[2][2] = {{1.0, 0.7}, {0.7, 1.0}}; Matrix m1(2, 2, &a1[0][0]); Matrix m2(2, 2, &a2[0][0]); Matrix m3(m1.covarToCorr()); CHECK_CLOSE(0.0, (m3 - m2).maxAbsValue(), 1.0e-15); double covs[2] = {9.0, 25.0}; Matrix m4(m3.corrToCovar(covs, 2)); CHECK_CLOSE(0.0, (m4 - m1).maxAbsValue(), 1.0e-15); } TEST(Matrix_coarseSum) { int a1[2][4] = {{1, 4, 2, 5}, {2, -10, 13, 11}}; int a2[2][2] = {{5, 7}, {-8, 24}}; int a3[1][4] = {{3, -6, 15, 16}}; int a4[2][1] = {{12}, {16}}; int a5[1][2] = {{-3, 31}}; Matrix m1(2, 4, &a1[0][0]); Matrix m2(2, 2, &a2[0][0]); Matrix m3(1, 4, &a3[0][0]); Matrix m4(2, 1, &a4[0][0]); Matrix m5(1, 2, &a5[0][0]); Matrix temp; m1.coarseSum(1, 1, &temp); CHECK(m1 == temp); m1.coarseSum(1, 2, &temp); CHECK(m2 == temp); m1.coarseSum(2, 1, &temp); CHECK_EQUAL(m3, temp); m1.coarseSum(1, 4, &temp); CHECK(m4 == temp); m1.coarseSum(2, 2, &temp); CHECK(m5 == temp); } TEST(Matrix_sqrt) { double a1[3][3] = {{4, 2.4, 9.0}, {2.4, 9, 10.5}, {9.0, 10.5, 25}}; Matrix m1(3, 3, &a1[0][0]); Matrix msqrt(m1.symFcn(FcnFunctor1(sqrt))); CHECK_CLOSE(0.0, (msqrt*msqrt - m1).maxAbsValue(), 1.0e-12); } TEST(Matrix_gen_eigensystem_double) { const unsigned ncycles = 5; const double eps = 1.0e-12; for (unsigned ndim = 2; ndim < 10; ++ndim) { vector > evalues(ndim); for (unsigned icycle = 0; icycle < ncycles; ++icycle) { Matrix m(ndim, ndim); for (unsigned row=0; row left(ndim, ndim); Matrix right(ndim, ndim); m.genEigen(&evalues[0], evalues.size(), &right, &left); bool firstComplex = true; for (unsigned i=0; i& prod = m.timesVector(right.data() + i*ndim, ndim); for (unsigned row=0; row& lprod = m.rowMultiply(left.data() + i*ndim, ndim); for (unsigned col=0; col > m2(m); vector > evec(ndim); for (unsigned j=0; j(*(right.data() + i*ndim + j), *(right.data() + (i+1)*ndim + j)); else evec[j] = complex(*(right.data() + (i-1)*ndim + j), -*(right.data() + i*ndim + j)); const Matrix >& prod = m2.timesVector(&evec[0], ndim); for (unsigned row=0; row(*(left.data() + i*ndim + j), -*(left.data() + (i+1)*ndim + j)); else evec[j] = complex(*(left.data() + (i-1)*ndim + j), *(left.data() + i*ndim + j)); const Matrix >& lprod = m2.rowMultiply(&evec[0], ndim); for (unsigned col=0; col > evalues(ndim); for (unsigned icycle = 0; icycle < ncycles; ++icycle) { Matrix m(ndim, ndim); for (unsigned row=0; row left(ndim, ndim); Matrix right(ndim, ndim); m.genEigen(&evalues[0], evalues.size(), &right, &left); bool firstComplex = true; for (unsigned i=0; i& prod = m.timesVector(right.data() + i*ndim, ndim); for (unsigned row=0; row& lprod = m.rowMultiply(left.data() + i*ndim, ndim); for (unsigned col=0; col > m2(m); vector > evec(ndim); for (unsigned j=0; j(*(right.data() + i*ndim + j), *(right.data() + (i+1)*ndim + j)); else evec[j] = complex(*(right.data() + (i-1)*ndim + j), -*(right.data() + i*ndim + j)); const Matrix >& prod = m2.timesVector(&evec[0], ndim); for (unsigned row=0; row(*(left.data() + i*ndim + j), -*(left.data() + (i+1)*ndim + j)); else evec[j] = complex(*(left.data() + (i-1)*ndim + j), *(left.data() + i*ndim + j)); const Matrix >& lprod = m2.rowMultiply(&evec[0], ndim); for (unsigned col=0; col m1(4, 4, &a1[0][0]); Matrix m2(4, 4); double eigenvalues[4]; m1.tdSymEigen(eigenvalues, 4, &m2); for (unsigned i=0; i<4; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-14); for (unsigned i=0; i<4; ++i) { double sum = 0.0; for (unsigned j=0; j<4; ++j) sum += a2[i][j]*a2[i][j]; sum = sqrt(sum); for (unsigned j=0; j<4; ++j) a2[i][j] /= sum; } for (unsigned i=0; i<4; ++i) { double ratio = 1.0; if (m2[i][0]*a2[i][0] < 0.0) ratio = -1.0; for (unsigned j=0; j<4; ++j) CHECK_CLOSE(ratio*a2[i][j], m2[i][j], 1.0e-13); } } TEST(Matrix_sym_eigensystem_double) { double a1[3][3] = {{4, 2.4, 9.0}, {2.4, 9, 10.5}, {9.0, 10.5, 25}}; double evalues[3] = {0.25663680988831909567, 4.6647571372875151856, 33.078606052824165719}; double a2[3][3] = {{-1.9833553555164647014, 1.5988627419919654596, 0.34836285960951042746}, {-0.6564919038536665325, -3.307143575299064186, 0.47079431584176874967}, { 1.0, 1.0, 1.0}}; Matrix m1(3, 3, &a1[0][0]); Matrix m2(3, 3); double eigenvalues[3]; m1.symEigen(eigenvalues, 3, NULL, EIGEN_RRR); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-14); m1.symEigen(eigenvalues, 3, NULL, EIGEN_D_AND_C); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-14); m1.symEigen(eigenvalues, 3, NULL); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-14); for (unsigned i=0; i<3; ++i) { double sum = 0.0; for (unsigned j=0; j<3; ++j) sum += a2[j][i]*a2[j][i]; sum = sqrt(sum); for (unsigned j=0; j<3; ++j) a2[j][i] /= sum; } m1.symEigen(eigenvalues, 3, &m2, EIGEN_RRR); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-13); for (unsigned i=0; i<3; ++i) { double ratio = 1.0; if (m2[i][0]*a2[0][i] < 0.0) ratio = -1.0; for (unsigned j=0; j<3; ++j) CHECK_CLOSE(ratio*a2[j][i], m2[i][j], 1.0e-13); } m1.symEigen(eigenvalues, 3, &m2, EIGEN_D_AND_C); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-14); for (unsigned i=0; i<3; ++i) { double ratio = 1.0; if (m2[i][0]*a2[0][i] < 0.0) ratio = -1.0; for (unsigned j=0; j<3; ++j) CHECK_CLOSE(ratio*a2[j][i], m2[i][j], 1.0e-14); } m1.symEigen(eigenvalues, 3, &m2); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-14); for (unsigned i=0; i<3; ++i) { double ratio = 1.0; if (m2[i][0]*a2[0][i] < 0.0) ratio = -1.0; for (unsigned j=0; j<3; ++j) CHECK_CLOSE(ratio*a2[j][i], m2[i][j], 1.0e-14); } Matrix m3(3, 3); m3.diagFill(eigenvalues, 3); CHECK((m2.T()*m3*m2 - m1).maxAbsValue() < 1.0e-13); Matrix m3U(3, 3, 1); CHECK((m2.T()*m2 - m3U).maxAbsValue() < 1.0e-13); CHECK((m2*m2.T() - m3U).maxAbsValue() < 1.0e-13); std::ostringstream os; m1.classId().write(os); m1.write(os); std::istringstream is(os.str()); gs::ClassId id(is, 1); Matrix readback; Matrix::restore(id, is, &readback); CHECK(readback == m1); } TEST(Matrix_nonzeros) { double ad[3][3] = {{4.0, 0.0, 9.0}, {-5.0, 9, 10.5}, {11.0, 3.5, 0.0}}; Matrix md(3, 3, &ad[0][0]); CHECK_EQUAL(7U, md.nonZeros()); } TEST(Matrix_least_squares_double) { double ad[3][2] = {{1., 7.}, {-2., 5.}, {3., 11.}}; Matrix md(3, 2, &ad[0][0]); double b[3] = {1., 2., 3.}; double sol[2]; CHECK(md.linearLeastSquares(b, 3, sol, 2)); CHECK_CLOSE(-11.0/61, sol[0], 1.0e-14); CHECK_CLOSE(52.0/183, sol[1], 1.0e-14); } TEST(Matrix_least_squares_float) { float ad[3][2] = {{1., 7.}, {-2., 5.}, {3., 11.}}; Matrix md(3, 2, &ad[0][0]); float b[3] = {1., 2., 3.}; float sol[2]; CHECK(md.linearLeastSquares(b, 3, sol, 2)); CHECK_CLOSE(-11.0/61, sol[0], 1.0e-6); CHECK_CLOSE(52.0/183, sol[1], 1.0e-6); } TEST(Matrix_weighted_least_squares_double) { const double eps = 1.0e-12; double ad[3][2] = {{1., 7.}, {-2., 5.}, {3., 11.}}; Matrix md(3, 2, &ad[0][0]); double b[3] = {1., 2., 3.}; double sol[2]; double chisq = 0.0; Matrix cov; Matrix vinv(3, 3, 0); vinv[0][0] = 1; vinv[1][1] = 1.0/4.0; vinv[2][2] = 1.0/25.0; double ecov[2][2] = {{2003/3598.0, -97/1799.0}, {-97/1799.0, 118/5397.0}}; Matrix expectedCov(2, 2, &ecov[0][0]); CHECK(md.weightedLeastSquares(b, 3, vinv, sol, 2, &chisq, &cov)); CHECK_CLOSE(-689/1799.0, sol[0], eps); CHECK_CLOSE(1172/5397.0, sol[1], eps); CHECK_CLOSE(800/5397.0, chisq, eps); CHECK((cov - expectedCov).maxAbsValue() < eps); } TEST(Matrix_underdeterminedLinearSystem_double_1) { double ad[3][3] = {{1., 1./5, 1./7}, {1./5, 2., 1./11}, {1./7, 1./11, 3.}}; Matrix V(3, 3, &ad[0][0]); Matrix constr(1, 3); constr[0][0] = 13; constr[0][1] = 17; constr[0][2] = -19; double cval[1] = {1.0 + 1.0/23.0}; double sol[3]; double chisq; Matrix A; const double eps = 1.0e-14; CHECK(constr.underdeterminedLinearSystem(cval, 1, V, sol, 3, &chisq, &A)); CHECK_CLOSE(31614.0/3960623, sol[0], eps); CHECK_CLOSE(80556.0/3960623, sol[1], eps); CHECK_CLOSE(-123810.0/3960623, sol[2], eps); CHECK_CLOSE(55440.0/91094329, chisq, eps); const Matrix& expectedA = V*constr.T()*(constr*V*constr.T()).symPDInv(); CHECK((A - expectedA).maxAbsValue() < eps); } TEST(Matrix_underdeterminedLinearSystem_double_2) { double ad[3][3] = {{1., 1./5, 1./7}, {1./5, 2., 1./11}, {1./7, 1./11, 3.}}; Matrix V(3, 3, &ad[0][0]); Matrix constr(2, 3); constr[0][0] = 13; constr[0][1] = 17; constr[0][2] = -19; constr[1][0] = 1; constr[1][1] = 1; constr[1][2] = 1; double cval[2] = {1.0 + 1.0/23.0, 1.5}; double sol[3]; double chisq; const double eps = 1.0e-14; CHECK(constr.underdeterminedLinearSystem(cval, 2, V, sol, 3, &chisq)); CHECK_CLOSE(1083077325.0/3483415504, sol[0], eps); CHECK_CLOSE(947968553.0/1741707752, sol[1], eps); CHECK_CLOSE(2246108825.0/3483415504, sol[2], eps); CHECK_CLOSE(105821683275.0/320474226368, chisq, eps); } TEST(Matrix_constrained_least_squares_double) { double ad[4][2] = {{1., 7.}, {-2., 5.}, {3., 11.}, {-1., 1.}}; Matrix md(4, 2, &ad[0][0]); double b[4] = {1., 2., 3., 4.}; double sol[2]; Matrix constr(1, 2); constr[0][0] = 1; constr[0][1] = 1; double d[1] = {-422.0/2099}; double chisq; // The following call will invoke the Lapack algorithm CHECK(md.constrainedLeastSquares(b, 4, constr, d, 1, sol, 2, &chisq)); CHECK_CLOSE(-1174.0/2099, sol[0], 1.0e-14); CHECK_CLOSE(752.0/2099, sol[1], 1.0e-14); CHECK_CLOSE(24710.0/2099, chisq, 1.0e-13); // The following call will invoke the "textbook" algorithm double sol0[2]; Matrix cov0, cov, proj, A; CHECK(md.constrainedLeastSquares(b, 4, constr, d, 1, sol, 2, &chisq, &cov, sol0, &cov0, &proj, &A)); CHECK_CLOSE(-1174.0/2099, sol[0], 1.0e-14); CHECK_CLOSE(752.0/2099, sol[1], 1.0e-14); CHECK_CLOSE(24710.0/2099, chisq, 1.0e-13); const Matrix& expSol = proj*Matrix(2, 1, sol0) + A*Matrix(1, 1, d); CHECK_CLOSE(expSol[0][0], sol[0], 1.0e-14); CHECK_CLOSE(expSol[1][0], sol[1], 1.0e-14); const Matrix& expCov = proj * cov0 * proj.T(); CHECK((cov - expCov).maxAbsValue() < 1.0e-13); } TEST(Matrix_linear_system_double) { static const double eps = 1.0e-14; double a1[2][2] = {{2, 3}, {6, 9}}; Matrix m1(2, 2, &a1[0][0]); double rhs1[2] = {1, 2}; double solution[3]; CHECK_EQUAL(false, m1.solveLinearSystem(rhs1, 2, solution)); double a2[3][3] = {{4.0, 7.0, 9.0}, {-5.0, 9, 10.5}, {11.0, 3.5, 25}}; Matrix m2(3, 3, &a2[0][0]); for (unsigned i=0; i<3; ++i) for (unsigned j=0; j<3; ++j) CHECK_EQUAL(m2[i][j], a2[i][j]); double x[3] = {1.0, 2.0, 3.0}; const Matrix& rhs2 = m2.timesVector(x, 3); CHECK_EQUAL(true, m2.solveLinearSystem(rhs2.data(), 3, solution)); for (unsigned j=0; j<3; ++j) CHECK_CLOSE(x[j], solution[j], eps); double y[3] = {7.0, 1.0, 3.0}; const Matrix& rhs3 = m2.timesVector(y, 3); CHECK_EQUAL(true, m2.solveLinearSystem(rhs3.data(), 3, solution)); for (unsigned j=0; j<3; ++j) CHECK_CLOSE(y[j], solution[j], eps); Matrix rhs4(3, 2); for (unsigned row=0; row<3; ++row) rhs4[row][0] = rhs2.data()[row]; for (unsigned row=0; row<3; ++row) rhs4[row][1] = rhs3.data()[row]; Matrix sol; m2.solveLinearSystems(rhs4, &sol); for (unsigned row=0; row<3; ++row) CHECK_CLOSE(x[row], sol[row][0], eps); for (unsigned row=0; row<3; ++row) CHECK_CLOSE(y[row], sol[row][1], eps); } TEST(Matrix_linear_system_float) { float a1[2][2] = {{2, 3}, {6, 9}}; Matrix m1(2, 2, &a1[0][0]); float rhs1[2] = {1, 2}; float solution[3]; CHECK_EQUAL(false, m1.solveLinearSystem(rhs1, 2, solution)); float a2[3][3] = {{4.0, 7.0, 9.0}, {-5.0, 9, 10.5}, {11.0, 3.5, 25}}; float x[3] = {1.0, 2.0, 3.0}; Matrix m2(3, 3, &a2[0][0]); Matrix rhs2(m2.timesVector(x, 3)); CHECK_EQUAL(true, m2.solveLinearSystem(rhs2.data(), 3, solution)); for (unsigned j=0; j<3; ++j) CHECK_CLOSE(x[j], solution[j], 1.0e-6); } TEST(Matrix_svd_double) { const double eps = 1.0e-12; double a1[3][4] = {{2., 5., -3., 8.}, {8., 3., 16., 7.}, {7., 4., 2., -11.}}; Matrix m1(3, 4, &a1[0][0]); double svalues[3]; Matrix U, V; m1.svd(svalues, 3, &U, &V, SVD_SIMPLE); const Matrix& UT = U.T(); const Matrix& m1svd = UT*diag(svalues, 3, 4)*V; const Matrix& delta = m1svd - m1; CHECK(delta.maxAbsValue() < eps); Matrix u3(3, 3, 1); Matrix u4(4, 4, 1); CHECK((UT*U - u3).maxAbsValue() < eps); CHECK((U*UT - u3).maxAbsValue() < eps); CHECK((V.T()*V - u4).maxAbsValue() < eps); CHECK((V*V.T() - u4).maxAbsValue() < eps); // Check singular vectors const Matrix& mt = m1.T(); for (unsigned i=0; i<3; ++i) { const Matrix& mv = m1.timesVector(V[i], V.nColumns()); Matrix u(3, 1, U[i]); CHECK((u*svalues[i] - mv).maxAbsValue() < eps); const Matrix& mtu = mt.timesVector(U[i], U.nColumns()); Matrix v(4, 1, V[i]); CHECK((v*svalues[i] - mtu).maxAbsValue() < eps); } CHECK(m1.timesVector(V[3], V.nColumns()).maxAbsValue() < eps); } TEST(Matrix_svd_double_2) { const double eps = 1.0e-12; double a1[3][4] = {{2., 5., -3., 8.}, {8., 3., 16., 7.}, {7., 4., 2., -11.}}; Matrix m1(3, 4, &a1[0][0]); double svalues[3]; Matrix U, V; m1.svd(svalues, 3, &U, &V, SVD_D_AND_C); const Matrix& UT = U.T(); const Matrix& m1svd = UT*diag(svalues, 3, 4)*V; const Matrix& delta = m1svd - m1; CHECK(delta.maxAbsValue() < 1.0e-12); Matrix u3(3, 3, 1); Matrix u4(4, 4, 1); CHECK((UT*U - u3).maxAbsValue() < eps); CHECK((U*UT - u3).maxAbsValue() < eps); CHECK((V.T()*V - u4).maxAbsValue() < eps); CHECK((V*V.T() - u4).maxAbsValue() < eps); // Check singular vectors const Matrix& mt = m1.T(); for (unsigned i=0; i<3; ++i) { const Matrix& mv = m1.timesVector(V[i], V.nColumns()); Matrix u(3, 1, U[i]); CHECK((u*svalues[i] - mv).maxAbsValue() < eps); const Matrix& mtu = mt.timesVector(U[i], U.nColumns()); Matrix v(4, 1, V[i]); CHECK((v*svalues[i] - mtu).maxAbsValue() < eps); } CHECK(m1.timesVector(V[3], V.nColumns()).maxAbsValue() < eps); } TEST(Matrix_svd_float) { float a1[3][4] = {{2., 5., -3., 8.}, {8., 3., 16., 7.}, {7., 4., 2., -3.}}; Matrix m1(3, 4, &a1[0][0]); float svalues[3]; Matrix U, V; m1.svd(svalues, 3, &U, &V, SVD_SIMPLE); const Matrix& m1svd = U.T()*diag(svalues, 3, 4)*V; const Matrix& delta = m1svd - m1; CHECK(delta.maxAbsValue() < 1.0e-5); } TEST(Matrix_svd_float_2) { float a1[3][4] = {{2., 5., -3., 8.}, {8., 3., 16., 7.}, {7., 4., 2., -3.}}; Matrix m1(3, 4, &a1[0][0]); float svalues[3]; Matrix U, V; m1.svd(svalues, 3, &U, &V, SVD_D_AND_C); const Matrix& m1svd = U.T()*diag(svalues, 3, 4)*V; const Matrix& delta = m1svd - m1; CHECK(delta.maxAbsValue() < 1.0e-5); } TEST(Matrix_eigensystem_float) { float a1[3][3] = {{4, 2.4, 9.0}, {2.4, 9, 10.5}, {9.0, 10.5, 25}}; float evalues[3] = {0.25663680988831909567, 4.6647571372875151856, 33.078606052824165719}; float a2[3][3] = {{-1.9833553555164647014, 1.5988627419919654596, 0.34836285960951042746}, {-0.6564919038536665325, -3.307143575299064186, 0.47079431584176874967}, { 1.0, 1.0, 1.0}}; Matrix m1(3, 3, &a1[0][0]); Matrix m2(3, 3); float eigenvalues[3]; m1.symEigen(eigenvalues, 3, NULL, EIGEN_RRR); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-5); m1.symEigen(eigenvalues, 3, NULL, EIGEN_D_AND_C); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-5); m1.symEigen(eigenvalues, 3, NULL); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-5); for (unsigned i=0; i<3; ++i) { float sum = 0.0; for (unsigned j=0; j<3; ++j) sum += a2[j][i]*a2[j][i]; sum = sqrt(sum); for (unsigned j=0; j<3; ++j) a2[j][i] /= sum; } m1.symEigen(eigenvalues, 3, &m2, EIGEN_RRR); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-4); for (unsigned i=0; i<3; ++i) { float ratio = 1.0; if (m2[i][0]*a2[0][i] < 0.0) ratio = -1.0; for (unsigned j=0; j<3; ++j) CHECK_CLOSE(ratio*a2[j][i], m2[i][j], 1.0e-4); } m1.symEigen(eigenvalues, 3, &m2, EIGEN_D_AND_C); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-5); for (unsigned i=0; i<3; ++i) { float ratio = 1.0; if (m2[i][0]*a2[0][i] < 0.0) ratio = -1.0; for (unsigned j=0; j<3; ++j) CHECK_CLOSE(ratio*a2[j][i], m2[i][j], 1.0e-5); } m1.symEigen(eigenvalues, 3, &m2); for (unsigned i=0; i<3; ++i) CHECK_CLOSE(evalues[i], eigenvalues[i], 1.0e-5); for (unsigned i=0; i<3; ++i) { float ratio = 1.0; if (m2[i][0]*a2[0][i] < 0.0) ratio = -1.0; for (unsigned j=0; j<3; ++j) CHECK_CLOSE(ratio*a2[j][i], m2[i][j], 1.0e-5); } } TEST(Matrix_symPSDefEffectiveRank) { for (unsigned dim=1; dim<20; ++dim) { const double ddim = dim; const Matrix m(dim, dim, 1); const double effR = m.symPSDefEffectiveRank(); CHECK_CLOSE(ddim, effR, 1.0e-12); } } TEST(Matrix_removeRowAndColumn) { const unsigned nrows = 4; const unsigned ncols = 3; Matrix m(nrows, ncols); for (unsigned row=0; row& rm = m.removeRowAndColumn( removedRow, removedCol); if (removedRow < nrows) CHECK_EQUAL(nrows-1, rm.nRows()); else CHECK_EQUAL(nrows, rm.nRows()); if (removedCol < ncols) CHECK_EQUAL(ncols-1, rm.nColumns()); else CHECK_EQUAL(ncols, rm.nColumns()); for (unsigned row=0; row removedRow) --torow; if (col > removedCol) --tocol; CHECK_EQUAL(m[row][col], rm.at(torow, tocol)); } } } } Index: trunk/npstat/nm/Matrix.icc =================================================================== --- trunk/npstat/nm/Matrix.icc (revision 690) +++ trunk/npstat/nm/Matrix.icc (revision 691) @@ -1,3026 +1,3046 @@ #include #include #include #include #include #include #include "geners/GenericIO.hh" #include "geners/IOException.hh" #include "npstat/nm/allocators.hh" #include "npstat/nm/lapack_interface.hh" #include "npstat/nm/PreciseType.hh" #include "npstat/nm/ComplexComparesAbs.hh" namespace npstat { template inline Matrix::Matrix() : data_(0), nrows_(0), ncols_(0), len_(0), isDiagonal_(false), diagonalityKnown_(false) { } template inline Matrix::Matrix(const Matrix& r) : data_(0), nrows_(r.nrows_), ncols_(r.ncols_), len_(r.len_), isDiagonal_(r.isDiagonal_), diagonalityKnown_(r.diagonalityKnown_) { if (len_) { if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; for (unsigned i=0; i inline Matrix::Matrix(Matrix&& r) : data_(0), nrows_(r.nrows_), ncols_(r.ncols_), len_(r.len_), isDiagonal_(r.isDiagonal_), diagonalityKnown_(r.diagonalityKnown_) { if (len_) { if (len_ <= Len) { data_ = local_; for (unsigned i=0; i template inline Matrix::Matrix(const Matrix& r) : data_(0), nrows_(r.nrows_), ncols_(r.ncols_), len_(r.len_), isDiagonal_(false), diagonalityKnown_(false) { if (len_) { if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; copyBuffer(data_, r.data_, len_); } else { nrows_ = 0U; ncols_ = 0U; } } template template inline Matrix& Matrix::setData( const Num2* data, const unsigned dataLength) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::setData: uninitialized matrix"); if (dataLength != len_) throw std::invalid_argument( "In npstat::Matrix::setData: incompatible input length"); assert(data); copyBuffer(data_, data, len_); diagonalityKnown_ = false; isDiagonal_ = false; return *this; } template template inline Matrix& Matrix::setRow( const unsigned row, const Num2* data, const unsigned dataLength) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::setRow: uninitialized matrix"); if (row >= nrows_) throw std::invalid_argument( "In npstat::Matrix::setRow: row number out of range"); if (dataLength != ncols_) throw std::invalid_argument( "In npstat::Matrix::setRow: incompatible input length"); assert(data); copyBuffer(data_+row*ncols_, data, len_); diagonalityKnown_ = false; isDiagonal_ = false; return *this; } template template inline Matrix& Matrix::setColumn( const unsigned col, const Num2* data, const unsigned dataLength) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::setColumn: uninitialized matrix"); if (col >= ncols_) throw std::invalid_argument( "In npstat::Matrix::setColumn: column number out of range"); if (dataLength != nrows_) throw std::invalid_argument( "In npstat::Matrix::setColumn: incompatible input length"); assert(data); Numeric* d = data_ + col; for (unsigned i=0; i template inline Matrix& Matrix::setFromTriplets( Iterator first, Iterator const last) { zeroOut(); for (; first != last; ++first) { const unsigned r = static_cast(first->row()); if (r >= nrows_) throw std::invalid_argument( "In npstat::Matrix::setFromTriplets: row number out of range"); const unsigned c = static_cast(first->col()); if (c >= ncols_) throw std::invalid_argument( "In npstat::Matrix::setFromTriplets: column number out of range"); data_[r*ncols_ + c] += first->value(); } return *this; } template template inline Matrix::Matrix( const Matrix& r, const unsigned rowMin, const unsigned rowMax, const unsigned columnMin, const unsigned columnMax) : data_(0), nrows_(rowMax - rowMin), ncols_(columnMax - columnMin), len_(nrows_*ncols_), isDiagonal_(false), diagonalityKnown_(false) { if (rowMin > rowMax || columnMin > columnMax) throw std::invalid_argument( "In npstat::Matrix subrange constructor: " "invalid subrange specification"); if (len_) { if (rowMax > r.nrows_) throw std::invalid_argument( "In npstat::Matrix subrange constructor: " "maximum row number out of range"); if (columnMax > r.ncols_) throw std::invalid_argument( "In npstat::Matrix subrange constructor: " "maximum column number out of range"); if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; for (unsigned row=0; row inline void Matrix::initialize(const unsigned nrows, const unsigned ncols) { const unsigned long sz = static_cast(nrows)*ncols; if (sz > static_cast(UINT_MAX)) throw std::invalid_argument( "In npstat::Matrix::validateRowsCols: requested matrix size is too large"); if (!sz) throw std::invalid_argument( "In npstat::Matrix::validateRowsCols: both number of " "rows and number of columns must be positive"); len_ = sz; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } template inline Matrix::Matrix(const unsigned nrows, const unsigned ncols) : data_(0), nrows_(nrows), ncols_(ncols), len_(0), isDiagonal_(false), diagonalityKnown_(false) { initialize(nrows, ncols); } template inline Matrix::Matrix(const unsigned nrows, const unsigned ncols, const Numeric* data) : data_(0), nrows_(nrows), ncols_(ncols), len_(0), isDiagonal_(false), diagonalityKnown_(false) { initialize(nrows, ncols); assert(data); for (unsigned i=0; i inline Matrix& Matrix::zeroOut() { const Numeric zero = Numeric(); for (unsigned i=0; i inline Matrix& Matrix::constFill(const Numeric c) { for (unsigned i=0; i inline Matrix::Matrix(const unsigned nrows, const unsigned ncols, const int initCode) : data_(0), nrows_(nrows), ncols_(ncols), len_(0), isDiagonal_(false), diagonalityKnown_(false) { initialize(nrows, ncols); zeroOut(); switch (initCode) { case 0: break; case 1: { if (nrows_ != ncols_) { if (len_ > Len) delete [] data_; throw std::invalid_argument( "In npstat::Matrix constructor: " "unit matrix must be square"); } const Numeric one(static_cast(1)); for (unsigned i=0; i Len) delete [] data_; throw std::invalid_argument( "In npstat::Matrix constructor: " "invalid initialization code"); } } } template inline Matrix::~Matrix() { if (data_ != local_) delete [] data_; } template inline Matrix& Matrix::uninitialize() { if (data_ != local_) delete [] data_; data_ = 0; nrows_ = 0; ncols_ = 0; len_ = 0; diagonalityKnown_ = false; isDiagonal_ = false; return *this; } template inline void Matrix::invertDiagonal(Matrix* result) const { assert(result != this); Numeric* rdata = result->data_; const unsigned len = nrows_*ncols_; const Numeric zero = Numeric(); const Numeric one = static_cast(1); for (unsigned i=0; idiagonalityKnown_ = true; result->isDiagonal_ = true; } template inline Matrix Matrix::symPDInv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symPDInv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symPDInv: matrix is not square"); Matrix result(nrows_, nrows_); if (isDiagonal()) invertDiagonal(&result); else invert_posdef_sym_matrix(data_, nrows_, result.data_); return result; } template inline Matrix Matrix::symInv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symInv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symInv: matrix is not square"); Matrix result(nrows_, nrows_); if (isDiagonal()) invertDiagonal(&result); else invert_sym_matrix(data_, nrows_, result.data_); return result; } template inline Matrix Matrix::inv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::inv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::inv: matrix is not square"); Matrix result(nrows_, nrows_); if (isDiagonal()) invertDiagonal(&result); else invert_general_matrix(data_, nrows_, result.data_); return result; } template + inline Matrix Matrix::leftInv() const + { + if (!nrows_) throw std::runtime_error( + "In npstat::Matrix::leftInv: uninitialized matrix"); + if (nrows_ < ncols_) throw std::domain_error( + "In npstat::Matrix::leftInv: more columns than rows"); + return TtimesThis().symInv().timesT(*this); + } + + template + inline Matrix Matrix::rightInv() const + { + if (!nrows_) throw std::runtime_error( + "In npstat::Matrix::rightInv: uninitialized matrix"); + if (nrows_ > ncols_) throw std::domain_error( + "In npstat::Matrix::rightInv: more rows than columns"); + return Ttimes(timesT().symInv()); + } + + template inline void Matrix::genEigen( typename GeneralizedComplex::type *eigvals, const unsigned lenEigvals, Matrix* rightEigenvectors, Matrix* leftEigenvectors) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::genEigen: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::genEigen: matrix is not square"); if (lenEigvals < ncols_) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "insufficient length of the eigenvalue buffer"); assert(eigvals); Numeric* right = 0; if (rightEigenvectors) { if (rightEigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "can not replace matrix with its own right eigenvectors"); if (rightEigenvectors->data_ == 0) *rightEigenvectors = *this; if (!(nrows_ == rightEigenvectors->nrows_ && nrows_ == rightEigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "incompatible matrix for storing right eigenvectors"); right = rightEigenvectors->data_; rightEigenvectors->isDiagonal_ = false; rightEigenvectors->diagonalityKnown_ = false; } Numeric* left = 0; if (leftEigenvectors) { if (leftEigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "can not replace matrix with its own left eigenvectors"); if (leftEigenvectors->data_ == 0) *leftEigenvectors = *this; if (!(nrows_ == leftEigenvectors->nrows_ && nrows_ == leftEigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "incompatible matrix for storing left eigenvectors"); left = leftEigenvectors->data_; leftEigenvectors->isDiagonal_ = false; leftEigenvectors->diagonalityKnown_ = false; } if (left || right) gen_matrix_eigensystem(data_, ncols_, eigvals, right, left); else gen_matrix_eigenvalues(data_, ncols_, eigvals); } template inline void Matrix::svd(Numeric* singularValues, const unsigned lenSingularValues, Matrix* U, Matrix* V, const SvdMethod m) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::svd: uninitialized matrix"); if (lenSingularValues < std::min(nrows_, ncols_)) throw std::invalid_argument( "In npstat::Matrix::svd: " "insufficient length of the singular values buffer"); assert(singularValues); assert(U); assert(V); if (U == this || V == this) throw std::invalid_argument( "In npstat::Matrix::svd: " "can not replace matrix with its own decomposition"); U->resize(nrows_, nrows_); V->resize(ncols_, ncols_); bool status = false; switch (m) { case SVD_SIMPLE: status = gen_matrix_svd(data_, nrows_, ncols_, U->data_, singularValues, V->data_); break; case SVD_D_AND_C: status = gen_matrix_svd_dc(data_, nrows_, ncols_, U->data_, singularValues, V->data_); break; default: assert(!"In npstat::Matrix::svd: incomplete switch " "statement. This is a bug. Please report."); } if (!status) throw std::runtime_error( "In npstat::Matrix::svd: algorithm did not converge"); } template inline void Matrix::symEigen(Numeric* eigenvalues, const unsigned lenEigenvalues, Matrix* eigenvectors, const EigenMethod m) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symEigen: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symEigen: matrix is not square"); if (lenEigenvalues < ncols_) throw std::invalid_argument( "In npstat::Matrix::symEigen: " "insufficient length of the eigenvalue buffer"); assert(eigenvalues); if (eigenvectors) { if (eigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::symEigen: " "can not replace matrix with its own eigenvectors"); if (eigenvectors->data_ == 0) *eigenvectors = *this; if (!(nrows_ == eigenvectors->nrows_ && nrows_ == eigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::symEigen: " "incompatible matrix for storing eigenvectors"); eigenvectors->isDiagonal_ = false; eigenvectors->diagonalityKnown_ = false; switch (m) { case EIGEN_SIMPLE: sym_matrix_eigensystem(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_D_AND_C: sym_matrix_eigensystem_dc(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_RRR: sym_matrix_eigensystem_rrr(data_, ncols_, eigenvalues, eigenvectors->data_); break; default: assert(!"In npstat::Matrix::symEigen: incomplete switch " "statement 1. This is a bug. Please report."); } } else switch (m) { case EIGEN_SIMPLE: sym_matrix_eigenvalues(data_, ncols_, eigenvalues); break; case EIGEN_D_AND_C: sym_matrix_eigenvalues_dc(data_, ncols_, eigenvalues); break; case EIGEN_RRR: sym_matrix_eigenvalues_rrr(data_, ncols_, eigenvalues); break; default: assert(!"In npstat::Matrix::symEigen: incomplete switch " "statement 2. This is a bug. Please report."); } } template inline void Matrix::tdSymEigen(Numeric* eigenvalues, const unsigned lenEigenvalues, Matrix* eigenvectors, const EigenMethod m) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::tdSymEigen: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::tdSymEigen: matrix is not square"); if (lenEigenvalues < ncols_) throw std::invalid_argument( "In npstat::Matrix::tdSymEigen: " "insufficient length of the eigenvalue buffer"); assert(eigenvalues); if (eigenvectors) { if (eigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::tdSymEigen: " "can not replace matrix with its own eigenvectors"); if (eigenvectors->data_ == 0) *eigenvectors = *this; if (!(nrows_ == eigenvectors->nrows_ && nrows_ == eigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::tdSymEigen: " "incompatible matrix for storing eigenvectors"); eigenvectors->isDiagonal_ = false; eigenvectors->diagonalityKnown_ = false; switch (m) { case EIGEN_SIMPLE: td_sym_matrix_eigensystem(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_D_AND_C: td_sym_matrix_eigensystem_dc(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_RRR: td_sym_matrix_eigensystem_rrr(data_, ncols_, eigenvalues, eigenvectors->data_); break; default: assert(!"In npstat::Matrix::tdSymEigen: incomplete switch " "statement 1. This is a bug. Please report."); } } else switch (m) { case EIGEN_SIMPLE: td_sym_matrix_eigenvalues(data_, ncols_, eigenvalues); break; case EIGEN_D_AND_C: td_sym_matrix_eigenvalues_dc(data_, ncols_, eigenvalues); break; case EIGEN_RRR: td_sym_matrix_eigenvalues_rrr(data_, ncols_, eigenvalues); break; default: assert(!"In npstat::Matrix::tdSymEigen: incomplete switch " "statement 2. This is a bug. Please report."); } } template inline Matrix Matrix::symPDEigenInv( const double tol, const bool extend, const EigenMethod eigenm) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symPDEigenInv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symPDEigenInv: matrix is not square"); if (!(tol > 0.0 && tol < 1.0)) throw std::domain_error( "In npstat::Matrix::symPDEigenInv: invalid tolerance parameter"); Matrix evec; std::vector evBuf(nrows_); Numeric* eigenvalues = &evBuf[0]; this->symEigen(eigenvalues, nrows_, &evec, eigenm); const unsigned nm1 = nrows_ - 1U; const double emax = eigenvalues[nm1]; if (emax <= 0.0) throw std::domain_error( "In npstat::Matrix::symPDEigenInv: " "largest eigenvalue is not positive"); const Numeric ecut = emax*tol; const Numeric replacementValue = 1.0/(emax*tol); const Numeric one = 1.0; eigenvalues[nm1] = one/eigenvalues[nm1]; for (unsigned i=0; i ecut) eigenvalues[i] = one/eigenvalues[i]; else if (extend) eigenvalues[i] = replacementValue; else eigenvalues[i] = 0; } Matrix dmat(nrows_, nrows_); dmat.diagFill(eigenvalues, nrows_); return dmat.bilinear(evec); } template inline double Matrix::symPSDefEffectiveRank( const double tol, const EigenMethod m, double* eigenSum) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symPSDefEffectiveRank: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symPSDefEffectiveRank: matrix is not square"); if (tol < 0.0) throw std::domain_error( "In npstat::Matrix::symPSDefEffectiveRank: " "tolerance can not be negative"); std::vector eigenvaluesVec(nrows_); this->symEigen(&eigenvaluesVec[0], nrows_, 0, m); const Numeric* eigenvalues = &eigenvaluesVec[0]; if (eigenvalues[nrows_ - 1] <= static_cast(0)) { if (eigenSum) *eigenSum = 0.0; return 0.0; } const Numeric cutoff = tol*eigenvalues[nrows_ - 1]; unsigned ifirst = nrows_ - 1U; for (; ifirst; --ifirst) if (eigenvalues[ifirst - 1] <= cutoff) break; long double sum = 0.0L, entropy = 0.0L; Numeric largest = Numeric(); for (unsigned i=ifirst; i largest) largest = eigenvalues[i]; entropy -= eigenvalues[i]*std::log(eigenvalues[i]); } if (eigenSum) *eigenSum = sum/static_cast(largest); return sum*std::exp(entropy/sum); } template template inline Matrix Matrix::symFcn(const Functor& fcn) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symFcn: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symFcn: matrix is not square"); Matrix evectors(nrows_, nrows_); Matrix evalues(nrows_, nrows_); symEigen(evalues.data_, nrows_, &evectors); const Numeric zero = Numeric(); for (unsigned i=1; i inline Matrix Matrix::pow(const unsigned degree) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::pow: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::pow: matrix is not square"); if (degree == 0U) return Matrix(nrows_, ncols_, 1); if (degree == 1U) return *this; Matrix result(*this); for (unsigned i=1; i inline bool Matrix::isSymmetric() const { if (nrows_ == 0) return false; if (nrows_ != ncols_) return false; for (unsigned i=1; i inline void Matrix::calcDiagonal() { isDiagonal_ = false; if (len_) { const Numeric null = Numeric(); isDiagonal_ = (nrows_ == ncols_); for (unsigned i=0; i inline bool Matrix::isDiagonal() const { if (!diagonalityKnown_) (const_cast(this))->calcDiagonal(); return isDiagonal_; } template inline bool Matrix::isZero() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::isZero: uninitialized matrix"); const Numeric zero = Numeric(); for (unsigned i=0; i inline bool Matrix::isAntiSymmetric() const { if (nrows_ == 0) return false; if (nrows_ != ncols_) return false; for (unsigned i=1; i inline Matrix Matrix::directSum(const Matrix& r) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::directSum: uninitialized matrix"); if (!r.nrows_) throw std::runtime_error( "In npstat::Matrix::directSum: uninitialized argument"); const unsigned newCols = ncols_+r.ncols_; Matrix result(nrows_+r.nrows_, newCols, 0); for (unsigned row=0; row inline Matrix Matrix::symmetrize() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symmetrize: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symmetrize: matrix is not square"); Matrix result(ncols_, nrows_); const Numeric two(static_cast(2)); for (unsigned i=0; i inline Matrix Matrix::antiSymmetrize() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::antiSymmetrize: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::antiSymmetrize: matrix is not square"); Matrix result(ncols_, nrows_); const Numeric two(static_cast(2)); const Numeric zero = Numeric(); for (unsigned i=0; i inline Matrix& Matrix::resize( const unsigned nrows, const unsigned ncols) { if (nrows != nrows_ || ncols != ncols_) { uninitialize(); len_ = nrows*ncols; if (len_) { nrows_ = nrows; ncols_ = ncols; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } } else { isDiagonal_ = false; diagonalityKnown_ = false; } return *this; } template inline Matrix& Matrix::tagAsDiagonal( const bool b) { if (!data_) throw std::runtime_error( "In npstat::Matrix::tagAsDiagonal: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::tagAsDiagonal: matrix is not square"); diagonalityKnown_ = true; isDiagonal_ = b; return *this; } template inline Matrix& Matrix::clearMainDiagonal() { const Numeric null = Numeric(); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::clearMainDiagonal: matrix is not square"); for (unsigned i=0; i inline Matrix& Matrix::Tthis() { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::Tthis: uninitialized matrix"); if (!(diagonalityKnown_ && isDiagonal_)) { if (nrows_ == ncols_) { for (unsigned i=0; i inline Matrix& Matrix::makeDiagonal() { const Numeric null = Numeric(); tagAsDiagonal(); for (unsigned i=0; i template inline Matrix& Matrix::diagFill(const Num2* data, const unsigned dataLen) { if (dataLen) assert(data); if (data_) { if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::diagFill: matrix is not square"); if (nrows_ != dataLen) throw std::invalid_argument( "In npstat::Matrix::diagFill: " "incompatible size of the input array"); } else if (dataLen) { nrows_ = dataLen; ncols_ = dataLen; len_ = dataLen*dataLen; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } if (dataLen) { if (!(diagonalityKnown_ && isDiagonal_)) zeroOut(); for (unsigned i=0; i(data[i]); diagonalityKnown_ = true; isDiagonal_ = true; } return *this; } template inline bool Matrix::isCompatible(const Matrix& r) const { if (!nrows_) return false; else return nrows_ == r.nrows_ && ncols_ == r.ncols_; } template inline Matrix& Matrix::operator=(const Matrix& r) { if (this == &r) return *this; if (data_) { if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator=: " "incompatible argument dimensions"); } else { nrows_ = r.nrows_; ncols_ = r.ncols_; len_ = r.len_; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } for (unsigned i=0; i inline Matrix& Matrix::operator=(Matrix&& r) { if (this == &r) return *this; if (data_) { if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator=: " "incompatible argument dimensions"); if (data_ != local_) { delete [] data_; data_ = 0; } } nrows_ = r.nrows_; ncols_ = r.ncols_; len_ = r.len_; if (len_ <= Len) { data_ = local_; for (unsigned i=0; i template inline Matrix& Matrix::operator=(const Matrix& r) { if (reinterpret_cast(this) == reinterpret_cast(&r)) return *this; if (data_) { if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator=: " "incompatible argument dimensions"); } else { nrows_ = r.nrows_; ncols_ = r.ncols_; len_ = r.len_; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } copyBuffer(data_, r.data_, len_); if (r.diagonalityKnown_ && r.isDiagonal_) { diagonalityKnown_ = true; isDiagonal_ = true; } else { diagonalityKnown_ = false; isDiagonal_ = false; } return *this; } template inline bool Matrix::operator==(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator==: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator==: " "incompatible argument dimensions"); for (unsigned i=0; i inline bool Matrix::operator!=(const Matrix& r) const { return !(*this == r); } template inline const Numeric* Matrix::operator[](const unsigned i) const { return data_ + i*ncols_; } template inline Numeric* Matrix::operator[](const unsigned i) { diagonalityKnown_ = false; return data_ + i*ncols_; } template inline Numeric Matrix::operator()( const unsigned row, const unsigned column) const { return data_[row*ncols_ + column]; } template inline Numeric Matrix::at( const unsigned row, const unsigned column) const { if (row >= nrows_) throw std::out_of_range("npstat::Matrix::at: " "row index out of range"); if (column >= ncols_) throw std::out_of_range("npstat::Matrix::at: " "column index out of range"); return data_[row*ncols_ + column]; } template inline Numeric Matrix::rowSum(const unsigned row) const { typedef typename PreciseType::type Precise; if (row >= nrows_) throw std::out_of_range("npstat::Matrix::rowSum: " "row index out of range"); Precise sum = Precise(); const Numeric* rData = data_ + row*ncols_; for (unsigned col=0; col(sum); } template inline Numeric Matrix::columnSum(const unsigned column) const { typedef typename PreciseType::type Precise; if (column >= ncols_) throw std::out_of_range("npstat::Matrix::columnSum: " "column index out of range"); Precise sum = Precise(); const Numeric* columnData = data_ + column; for (unsigned row=0; row(sum); } template inline Matrix& Matrix::set( const unsigned row, const unsigned column, const Numeric value) { if (row >= nrows_) throw std::out_of_range("npstat::Matrix::set: " "row index out of range"); if (column >= ncols_) throw std::out_of_range("npstat::Matrix::set: " "column index out of range"); data_[row*ncols_ + column] = value; diagonalityKnown_ = false; return *this; } template inline Matrix Matrix::operator+(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator+: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator+: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline Matrix Matrix::hadamardProduct(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::hadamardProduct: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::hadamardProduct: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline Matrix Matrix::hadamardRatio(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::hadamardRatio: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::hadamardRatio: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i(0)) throw std::domain_error( "In npstat::Matrix::hadamardRatio: division by zero"); result.data_[i] = data_[i]/r.data_[i]; } return result; } template inline void Matrix::plus(const Matrix& r, Matrix* result) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::plus: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::plus: " "incompatible argument dimensions"); assert(result); if (result != this) result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; i inline Matrix Matrix::operator-(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator-: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator-: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline void Matrix::minus(const Matrix& r, Matrix* result) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::minus: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::minus: " "incompatible argument dimensions"); assert(result); result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; i inline Matrix Matrix::T() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::T: uninitialized matrix"); Matrix result(ncols_, nrows_); if (diagonalityKnown_ && isDiagonal_) copyBuffer(result.data_, data_, nrows_*ncols_); else transposeBuffer(result.data_, data_, nrows_, ncols_); result.isDiagonal_ = isDiagonal_; result.diagonalityKnown_ = diagonalityKnown_; return result; } template inline Matrix Matrix::TtimesThis() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::TtimesThis: uninitialized matrix"); Matrix result(ncols_, ncols_); if (isDiagonal()) { result.zeroOut(); for (unsigned row=0; row inline Matrix Matrix::bilinearT(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::bilinearT: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::bilinearT: matrix is not square"); if (ncols_ != r.ncols_) throw std::domain_error( "In npstat::Matrix::bilinearT: incompatible matrix dimensions"); Matrix result(r.nrows_, r.nrows_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row Len) buf = new Numeric[nrows_]; for (unsigned col=0; col inline Matrix Matrix::bilinear(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::bilinear: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::bilinear: matrix is not square"); if (ncols_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::bilinear: incompatible matrix dimensions"); Matrix result(r.ncols_, r.ncols_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row Len) buf = new Numeric[nrows_]; for (unsigned col=0; col inline Matrix Matrix::Ttimes(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::Ttimes: uninitialized matrix"); if (nrows_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::Ttimes: incompatible matrix dimensions"); Matrix result(ncols_, r.ncols_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row inline Matrix Matrix::timesT(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::timesT: uninitialized matrix"); if (ncols_ != r.ncols_) throw std::domain_error( "In npstat::Matrix::timesT: incompatible matrix dimensions"); Matrix result(nrows_, r.nrows_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row inline Matrix Matrix::timesT() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesT: uninitialized matrix"); Matrix result(nrows_, nrows_); if (isDiagonal()) { result.zeroOut(); for (unsigned row=0; row inline double Matrix::frobeniusNorm() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::frobeniusNorm: uninitialized matrix"); long double sum = 0.0L; for (unsigned i=0; i inline Matrix Matrix::covarToCorr() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::covarToCorr: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::covarToCorr: matrix is not square"); Matrix result(ncols_, nrows_); Numeric* rdata = result.data_; for (unsigned row=0; row(1); return result; } template inline Matrix Matrix::corrToCovar( const double* variances, const unsigned lenVariances) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::corrToCovar: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::corrToCovar: matrix is not square"); if (lenVariances != nrows_) throw std::invalid_argument( "In npstat::Matrix::corrToCovar: " "incompatible length of covariance array"); assert(variances); for (unsigned row=0; row inline void Matrix::makeCoarse( const unsigned n, const unsigned m, Matrix* result) const { assert(result); if (result == this) throw std::invalid_argument( "In npstat::Matrix::makeCoarse: " "can not replace matrix with its own coarse version"); const unsigned newRows = nrows_ / n; const unsigned newCols = ncols_ / m; result->resize(newRows, newCols); if (m == 1U && n == 1U) { copyBuffer(result->data_, data_, len_); result->isDiagonal_ = isDiagonal_; result->diagonalityKnown_ = diagonalityKnown_; } else if (m == 1U) { // Sum the rows for (unsigned row=0; rowdata_ + row*ncols_; for (unsigned col=0; coldata_ + row*newCols; const Numeric* off = data_ + row*ncols_; for (unsigned col=0; coldata_ + row*newCols; const unsigned rstart = row*n; const unsigned rstop = rstart + n; for (unsigned col=0; col inline void Matrix::coarseSum( const unsigned n, const unsigned m, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::coarseSum: uninitialized matrix"); if (!(n && m)) throw std::invalid_argument( "In npstat::Matrix::coarseSum: " "number of divisions must be positive"); if (nrows_ % n) throw std::invalid_argument( "In npstat::Matrix::coarseSum: " "n is not a divisor of the number or rows"); if (ncols_ % m) throw std::invalid_argument( "In npstat::Matrix::coarseSum: " "m is not a divisor of the number or columns"); makeCoarse(n, m, result); } template inline void Matrix::coarseAverage( const unsigned n, const unsigned m, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::coarseAverage: uninitialized matrix"); const unsigned nm = n*m; if (!nm) throw std::invalid_argument( "In npstat::Matrix::coarseAverage: " "number of divisions must be positive"); if (nrows_ % n) throw std::invalid_argument( "In npstat::Matrix::coarseAverage: " "n is not a divisor of the number or rows"); if (ncols_ % m) throw std::invalid_argument( "In npstat::Matrix::coarseAverage: " "m is not a divisor of the number or columns"); makeCoarse(n, m, result); if (nm > 1U) *result /= static_cast(nm); } template inline Matrix Matrix::removeRowAndColumn( const unsigned rRow, const unsigned rColumn) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::removeRowAndColumn: uninitialized matrix"); unsigned newRows = nrows_; unsigned newCols = ncols_; if (rRow < nrows_) --newRows; if (rColumn < ncols_) --newCols; if (!(newRows && newCols)) throw std::invalid_argument( "In npstat::Matrix::removeRowAndColumn: can not remove the only " "row and/or column"); Matrix result(newRows, newCols); unsigned fromrow = 0U; for (unsigned torow=0U; torow inline Matrix Matrix::outer(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::outer: uninitialized matrix"); const unsigned nrows = nrows_*r.nrows_; const unsigned ncols = ncols_*r.ncols_; Matrix result(nrows, ncols); for (unsigned tr=0; tr inline Matrix Matrix::operator*(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator*: uninitialized matrix"); if (ncols_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::operator*: incompatible matrix dimensions"); Matrix result(nrows_, r.ncols_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row inline void Matrix::times(const Matrix& r, Matrix* result) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::times: uninitialized matrix"); if (ncols_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::times: incompatible matrix dimensions"); assert(result); if (result == this || result == &r) throw std::invalid_argument( "In npstat::Matrix::times: " "can not have result set to one of the arguments"); result->resize(nrows_, r.ncols_); Numeric* const rdata = result->data_; const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result->zeroOut(); for (unsigned row=0; rowdiagonalityKnown_ = true; result->isDiagonal_ = true; } else if (thisDiag) { // scale each row of the matrix r for (unsigned row=0; row template inline Matrix Matrix::timesVector( const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesVector: uninitialized matrix"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::timesVector: incompatible vector length"); assert(data); Matrix result(nrows_, 1U); Numeric* res = result[0]; if (isDiagonal()) { for (unsigned row=0; row template inline Numeric Matrix::timesVector( const unsigned row, const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesVector: uninitialized matrix"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::timesVector: incompatible vector length"); if (row >= nrows_) throw std::invalid_argument( "In npstat::Matrix::timesVector: row number out of range"); assert(data); if (isDiagonal()) return data_[row*ncols_ + row]*data[row]; else { const Numeric* prow = data_ + row*ncols_; Numeric sum = Numeric(); for (unsigned k=0; k template inline void Matrix::timesVector( const Num2* data, const unsigned dataLen, Num3* res, const unsigned resultLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesVector: uninitialized matrix"); if (nrows_ > resultLen) throw std::invalid_argument( "In npstat::Matrix::timesVector: " "insufficient length of the output buffer"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::timesVector: incompatible vector length"); assert(data); assert(res); if (isDiagonal()) { for (unsigned row=0; row template inline Matrix Matrix::rowMultiply( const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::rowMultiply: uninitialized matrix"); if (nrows_ != dataLen) throw std::domain_error( "In npstat::Matrix::rowMultiply: incompatible row length"); assert(data); Matrix result(1U, ncols_); Numeric* res = result[0]; if (isDiagonal()) { for (unsigned col=0; col template inline Numeric Matrix::rowMultiply( const unsigned col, const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::rowMultiply: uninitialized matrix"); if (nrows_ != dataLen) throw std::domain_error( "In npstat::Matrix::rowMultiply: incompatible row length"); if (col >= ncols_) throw std::invalid_argument( "In npstat::Matrix::rowMultiply: column number out of range"); assert(data); if (isDiagonal()) return data[col]*data_[col*ncols_ + col]; else { const Numeric* pcol = data_ + col; Numeric sum = Numeric(); for (unsigned k=0; k template inline void Matrix::rowMultiply( const Num2* data, const unsigned dataLen, Num3* res, const unsigned resultLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::rowMultiply: uninitialized matrix"); if (ncols_ > resultLen) throw std::invalid_argument( "In npstat::Matrix::rowMultiply: " "insufficient length of the output buffer"); if (nrows_ != dataLen) throw std::domain_error( "In npstat::Matrix::rowMultiply: incompatible row length"); assert(data); assert(res); if (isDiagonal()) { for (unsigned col=0; col template inline Numeric Matrix::bilinear( const Num2* data, const unsigned dataLen) const { typedef typename PreciseType::type Precise; if (!nrows_) throw std::runtime_error( "In npstat::Matrix::bilinear: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::bilinear: matrix is not square"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::bilinear: incompatible vector length"); assert(data); Precise sum = Precise(); if (isDiagonal()) { for (unsigned row=0; row(sum); } template inline bool Matrix::solveLinearSystem( const Numeric* rhs, const unsigned lenRhs, Numeric* solution) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::solveLinearSystem: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::solveLinearSystem: matrix is not square"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::solveLinearSystem: " "incompatible dimensionality of the right hand side vector"); assert(rhs); assert(solution); return solve_linear_system(data_, lenRhs, rhs, solution); } template inline bool Matrix::underdeterminedLinearSystem( const Numeric* rhs, const unsigned lenRhs, const Matrix& V, Numeric* solution, const unsigned lenSolution, Numeric* resultNormSq, Matrix* pa) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::underdeterminedLinearSystem: " "uninitialized matrix"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "incompatible dimensionality of the solution vector"); if (nrows_ >= ncols_) throw std::domain_error( "In npstat::Matrix::underdeterminedLinearSystem: " "the system is not underdetermined"); if (!V.isSquare()) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "covariance matrix is not square"); if (V.nRows() != ncols_) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "incompatible covariance matrix dimensionality"); assert(rhs); assert(solution); const Matrix& vct = V.timesT(*this); const Matrix& cvctinv = (*this * vct).symPDInv(); const Matrix& A = vct*cvctinv; A.timesVector(rhs, lenRhs, solution, lenSolution); if (pa) { if (!A.isCompatible(*pa)) pa->uninitialize(); *pa = A; } if (resultNormSq) { const Matrix& B = this->Ttimes(cvctinv); *resultNormSq = Numeric(); for (unsigned i=0; i inline bool Matrix::linearLeastSquares( const Numeric* rhs, const unsigned lenRhs, Numeric* solution, const unsigned lenSolution) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::linearLeastSquares: uninitialized matrix"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::linearLeastSquares: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::linearLeastSquares: " "incompatible dimensionality of the solution vector"); if (nrows_ <= ncols_) throw std::domain_error( "In npstat::Matrix::linearLeastSquares: " "the system is not overdetermined"); assert(rhs); assert(solution); return linear_least_squares(data_, nrows_, ncols_, rhs, solution); } template inline bool Matrix::weightedLeastSquares( const Numeric* rhsVec, const unsigned lenRhs, const Matrix& vinv, Numeric* solution, const unsigned lenSolution, Numeric* resultChiSquare, Matrix* resultCovarianceMatrix) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::weightedLeastSquares: uninitialized matrix"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "incompatible dimensionality of the solution vector"); if (vinv.nrows_ != nrows_) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "incompatible inverse covariance matrix dimensionality"); if (!vinv.isSquare()) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "inverse covariance matrix is not square"); if (nrows_ <= ncols_) throw std::domain_error( "In npstat::Matrix::weightedLeastSquares: " "the system is not overdetermined"); assert(rhsVec); assert(solution); Matrix rhs(lenRhs, 1, rhsVec); const Matrix& rcov = vinv.bilinear(*this).symPDInv(); const Matrix& sol = rcov * Ttimes(vinv * rhs); assert(sol.nrows_ == lenSolution); copyBuffer(solution, sol.data_, lenSolution); if (resultChiSquare) { rhs -= *this * sol; *resultChiSquare = vinv.bilinear(rhs.data_, lenRhs); } if (resultCovarianceMatrix) { if (!rcov.isCompatible(*resultCovarianceMatrix)) resultCovarianceMatrix->uninitialize(); *resultCovarianceMatrix = rcov; } return true; } template inline bool Matrix::constrainedLeastSquares( const Numeric* rhs1, const unsigned lenRhs1, const Matrix& B, const Numeric* rhs2, const unsigned lenRhs2, Numeric* solution, const unsigned lenSolution, Numeric* resultChiSquare, Matrix* resultCovarianceMatrix, Numeric* unconstrainedSolution, Matrix* unconstrainedCovmat, Matrix* projectionMatrix, Matrix* pA) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::constrainedLeastSquares: uninitialized matrix"); if (nrows_ != lenRhs1) throw std::invalid_argument( "In npstat::Matrix::constrainedLeastSquares: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::constrainedLeastSquares: " "incompatible dimensionality of the solution vector"); if (!B.nrows_) throw std::runtime_error( "In npstat::Matrix::constrainedLeastSquares: " "uninitialized constraint matrix"); if (B.ncols_ != lenSolution) throw std::runtime_error( "In npstat::Matrix::constrainedLeastSquares: " "incompatible constraint matrix"); if (B.nrows_ != lenRhs2) throw std::invalid_argument( "In npstat::Matrix::constrainedLeastSquares: " "incompatible dimensionality of the constraint vector"); if (nrows_ + lenRhs2 <= ncols_) throw std::domain_error( "In npstat::Matrix::constrainedLeastSquares: " "the system is not overdetermined"); assert(rhs1); assert(rhs2); assert(solution); bool status = false; if (resultCovarianceMatrix || unconstrainedSolution || unconstrainedCovmat || projectionMatrix || pA) { if (nrows_ < ncols_) throw std::domain_error( "In npstat::Matrix::constrainedLeastSquares: " "unconstrained system is underdetermined"); // Solution of the unconstrained problem const Matrix& cov0 = TtimesThis().symPDInv(); const Matrix& sol0 = cov0 * T().timesVector(rhs1, lenRhs1); assert(sol0.nrows_ == lenSolution); if (unconstrainedSolution) copyBuffer(unconstrainedSolution, sol0.data_, lenSolution); if (unconstrainedCovmat) { if (!unconstrainedCovmat->isCompatible(cov0)) unconstrainedCovmat->uninitialize(); *unconstrainedCovmat = cov0; } // Solve the constrained problem const Matrix& tmp = cov0.timesT(B); const Matrix& A = tmp * (B * tmp).symPDInv(); const Matrix& CP = Matrix(lenSolution, lenSolution, 1) - A * B; const Matrix& sol = CP * sol0 + A.timesVector(rhs2, lenRhs2); copyBuffer(solution, sol.data_, lenSolution); if (projectionMatrix) { if (!projectionMatrix->isCompatible(CP)) projectionMatrix->uninitialize(); *projectionMatrix = CP; } if (pA) { if (!pA->isCompatible(A)) pA->uninitialize(); *pA = A; } if (resultCovarianceMatrix) { // Evaluate the result covariance matrix const Matrix& rcov = cov0.bilinearT(CP); if (!rcov.isCompatible(*resultCovarianceMatrix)) resultCovarianceMatrix->uninitialize(); *resultCovarianceMatrix = rcov; } status = true; } else status = constrained_least_squares(data_, nrows_, ncols_, rhs1, B.data_, B.nrows_, B.ncols_, rhs2, solution); if (resultChiSquare) { if (status) { long double chisq = 0.0L; for (unsigned i=0; i(chisq); } else *resultChiSquare = static_cast(-1); } return status; } template inline bool Matrix::solveLinearSystems( const Matrix& rhs, Matrix* solution) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::solveLinearSystems: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::solveLinearSystems: matrix is not square"); if (nrows_ != rhs.nrows_) throw std::invalid_argument( "In npstat::Matrix::solveLinearSystems: " "incompatible dimensionality of the right hand side matrix"); assert(solution); if (solution == this || solution == &rhs) throw std::invalid_argument( "In npstat::Matrix::solveLinearSystems: " "can not have result set to one of the arguments"); solution->resize(rhs.nrows_, rhs.ncols_); const bool ok = solve_linear_systems(data_, rhs.nrows_, rhs.ncols_, rhs.data_, solution->data_); if (!ok) { const Numeric zero = Numeric(); const unsigned len = solution->length(); Numeric* sdata = solution->data_; for (unsigned i=0; i inline Numeric Matrix::minValue() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::minValue: uninitialized matrix"); Numeric minval = data_[0]; for (unsigned i=1; i::less(data_[i], minval)) minval = data_[i]; return minval; } template inline std::pair Matrix::argmin() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::minValue: uninitialized matrix"); Numeric minval = data_[0]; std::pair location(0U, 0U); for (unsigned i=1; i::less(data_[i], minval)) { minval = data_[i]; location = std::pair(i/ncols_, i%ncols_); } return location; } template inline Numeric Matrix::maxValue() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxValue: uninitialized matrix"); Numeric maxval = data_[0]; for (unsigned i=1; i::less(maxval, data_[i])) maxval = data_[i]; return maxval; } template inline std::pair Matrix::argmax() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxValue: uninitialized matrix"); Numeric maxval = data_[0]; std::pair location(0U, 0U); for (unsigned i=1; i::less(maxval, data_[i])) { maxval = data_[i]; location = std::pair(i/ncols_, i%ncols_); } return location; } template inline Numeric Matrix::maxAbsValue() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxAbsValue: uninitialized matrix"); Numeric maxval = std::abs(data_[0]); for (unsigned i=1; i maxval) maxval = std::abs(data_[i]); return maxval; } template inline std::pair Matrix::argmaxAbs() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxAbsValue: uninitialized matrix"); Numeric maxval = std::abs(data_[0]); std::pair location(0U, 0U); for (unsigned i=1; i maxval) { maxval = std::abs(data_[i]); location = std::pair(i/ncols_, i%ncols_); } return location; } template inline Numeric Matrix::tr() const { typedef typename PreciseType::type Precise; if (!nrows_) throw std::runtime_error( "In npstat::Matrix::tr: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::tr: matrix is not square"); Precise sum = Precise(); for (unsigned row=0; row(sum); } template inline Numeric Matrix::productTr(const Matrix& r) const { typedef typename PreciseType::type Precise; if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::productTr: uninitialized matrix"); if (ncols_ != r.nrows_ || nrows_ != r.ncols_) throw std::domain_error( "In npstat::Matrix::productTr: incompatible matrix dimensions"); Precise dsum = Precise(); for (unsigned row=0; row(dsum); } template inline Numeric Matrix::det() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::det: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::det: matrix is not square"); Numeric dvalue = Numeric(); if (nrows_ == 1) dvalue = data_[0]; else if (nrows_ == 2) dvalue = data_[0]*data_[3] - data_[1]*data_[2]; else if (nrows_ == 3) dvalue = -(data_[2]*data_[4]*data_[6]) + data_[1]*data_[5]*data_[6] + data_[2]*data_[3]*data_[7] - data_[0]*data_[5]*data_[7] - data_[1]*data_[3]*data_[8] + data_[0]*data_[4]*data_[8]; else if (nrows_ == 4) dvalue = data_[3]*data_[6]*data_[9]*data_[12] - data_[2]*data_[7]*data_[9]*data_[12] - data_[3]*data_[5]*data_[10]*data_[12] + data_[1]*data_[7]*data_[10]*data_[12] + data_[2]*data_[5]*data_[11]*data_[12] - data_[1]*data_[6]*data_[11]*data_[12] - data_[3]*data_[6]*data_[8]*data_[13] + data_[2]*data_[7]*data_[8]*data_[13] + data_[3]*data_[4]*data_[10]*data_[13] - data_[0]*data_[7]*data_[10]*data_[13] - data_[2]*data_[4]*data_[11]*data_[13] + data_[0]*data_[6]*data_[11]*data_[13] + data_[3]*data_[5]*data_[8]*data_[14] - data_[1]*data_[7]*data_[8]*data_[14] - data_[3]*data_[4]*data_[9]*data_[14] + data_[0]*data_[7]*data_[9]*data_[14] + data_[1]*data_[4]*data_[11]*data_[14] - data_[0]*data_[5]*data_[11]*data_[14] - data_[2]*data_[5]*data_[8]*data_[15] + data_[1]*data_[6]*data_[8]*data_[15] + data_[2]*data_[4]*data_[9]*data_[15] - data_[0]*data_[6]*data_[9]*data_[15] - data_[1]*data_[4]*data_[10]*data_[15] + data_[0]*data_[5]*data_[10]*data_[15]; else dvalue = lu_decomposition_matrix_det(data_, nrows_); return dvalue; } template inline Matrix Matrix::operator*(const Numeric r) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator*: uninitialized matrix"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline void Matrix::times(const Numeric r, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::times: uninitialized matrix"); assert(result); if (result != this) result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; idiagonalityKnown_ = true; result->isDiagonal_ = true; } } template inline Matrix Matrix::operator/(const Numeric r) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator/: uninitialized matrix"); if (r == static_cast(0)) throw std::domain_error( "In npstat::Matrix::operator/: division by zero"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline void Matrix::over(const Numeric r, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::over: uninitialized matrix"); if (r == static_cast(0)) throw std::domain_error( "In npstat::Matrix::over: division by zero"); assert(result); if (result != this) result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; idiagonalityKnown_ = true; result->isDiagonal_ = true; } } template inline Matrix Matrix::operator+() const { return *this; } template inline Matrix Matrix::operator-() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator-(unary): uninitialized matrix"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline Matrix& Matrix::operator*=(const Numeric r) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator*=: uninitialized matrix"); for (unsigned i=0; i 0 if (!(diagonalityKnown_ && isDiagonal_)) diagonalityKnown_ = false; return *this; } template inline Matrix& Matrix::operator/=(const Numeric r) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator/=: uninitialized matrix"); if (r == static_cast(0)) throw std::domain_error( "In npstat::Matrix::operator/=: division by zero"); for (unsigned i=0; i 0 if (!(diagonalityKnown_ && isDiagonal_)) diagonalityKnown_ = false; return *this; } template inline Matrix& Matrix::operator+=(const Matrix& r) { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator+=: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator+=: " "incompatible argument dimensions"); for (unsigned i=0; i inline Matrix& Matrix::operator-=(const Matrix& r) { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator-=: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator-=: " "incompatible argument dimensions"); for (unsigned i=0; i inline const char* Matrix::classname() { static const std::string name( gs::template_class_name("npstat::Matrix")); return name.c_str(); } template inline bool Matrix::write(std::ostream& os) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::write: uninitialized matrix"); assert(data_); gs::write_pod(os, nrows_); gs::write_pod(os, ncols_); gs::write_array(os, data_, len_); return !os.fail(); } template inline void Matrix::restore(const gs::ClassId& id, std::istream& in, Matrix* m) { static const gs::ClassId current( gs::ClassId::makeId >()); current.ensureSameId(id); assert(m); m->uninitialize(); gs::read_pod(in, &m->nrows_); gs::read_pod(in, &m->ncols_); m->len_ = m->nrows_ * m->ncols_; if (m->len_ <= Len) m->data_ = m->local_; else m->data_ = new Numeric[m->len_]; gs::read_array(in, m->data_, m->len_); if (in.fail()) { m->uninitialize(); throw gs::IOReadFailure("In npstat::Matrix::restore: " "input stream failure"); } } template inline unsigned Matrix::nonZeros() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::nonZeros: uninitialized matrix"); const Numeric null = Numeric(); unsigned cnt = 0U; for (unsigned i=0; i inline Matrix diag(const Numeric* data, const unsigned nrows) { Matrix m(nrows, nrows); m.diagFill(data, nrows); return m; } template inline Matrix diag(const Numeric* data, const unsigned nrows, const unsigned ncols) { Matrix m(nrows, ncols, 0); assert(data); const unsigned imax = std::min(nrows, ncols); for (unsigned i=0; i inline std::ostream& operator<<( std::ostream& os, const npstat::Matrix& m) { // Figure out the largest string length needed to print an element const unsigned nrows = m.nRows(); const unsigned ncols = m.nColumns(); const unsigned len = nrows*ncols; if (!len) throw std::runtime_error( "In ::operator<<: can not write uninitialized npstat::Matrix object"); const N* data = m.data(); unsigned maxlen = 0; for (unsigned i=0; i maxlen) maxlen = siz; } ++maxlen; for (unsigned i=0; i #include #include #include #include "geners/ClassId.hh" #include "npstat/nm/EigenMethod.hh" #include "npstat/nm/SvdMethod.hh" #include "npstat/nm/GeneralizedComplex.hh" namespace npstat { /** // A simple helper class for matrix manipulations. Depending on how much // space is provided with the "Len" parameter, the data will be placed // either on the stack or on the heap. Do not set "Len" to 0, the code // assumes that the stack space is available for at least one element. // // C storage convention is used for internal data. In the element access // operations, array bounds are not checked. // // A number of operations (solving linear systems, calculating eigenvalues // and eigenvectors, singular value decomposition, etc) are performed // by calling appropriate routines from LAPACK. This usually limits // the Numeric template parameter types for which these operations are // available to float and double. If the operation is unavailable for // the template parameter used, std::invalid_argument exception is raised. // // Note that for simple matrix operations this class is likely to be // slower than matrix classes based on expression templates (such as // those in boost uBLAS or in Blitz++). If speed is really important // for your calculations, consider using a dedicated matrix library. */ template class Matrix { template friend class Matrix; public: typedef Numeric value_type; /** // Default constructor creates an unitialized matrix // which can be assigned from other matrices */ Matrix(); /** // This constructor creates an unitialized matrix // which can be assigned element-by-element or from // another matrix with the same dimensions */ Matrix(unsigned nrows, unsigned ncols); /** // This constructor initializes the matrix as follows: // // initCode = 0: all elements are initialized to 0. // // initCode = 1: matrix must be square; diagonal elements are // initialized to 1, all other elements to 0. */ Matrix(unsigned nrows, unsigned ncols, int initCode); /** // This constructor initializes the matrix from the // given 1-d array using C storage conventions */ Matrix(unsigned nrows, unsigned ncols, const Numeric* data); /** Copy constructor */ Matrix(const Matrix&); /** Move constructor */ Matrix(Matrix&&); /** Converting copy constructor */ template Matrix(const Matrix&); /** // Constructor from a subrange of another matrix. The minimum row // and column numbers are included, maximum numbers are excluded. */ template Matrix(const Matrix&, unsigned rowMin, unsigned rowMax, unsigned columnMin, unsigned columnMax); ~Matrix(); /** // Assignment operator. The matrix on the left must either be // in an uninitialized state or have compatible dimensions with // the matrix on the right. */ Matrix& operator=(const Matrix&); /** Move assignment operator */ Matrix& operator=(Matrix&&); /** Converting assignment operator */ template Matrix& operator=(const Matrix&); /** Set from triplets (for compatibility with Eigen) */ template Matrix& setFromTriplets(Iterator first, Iterator last); /** Set all data from an external array */ template Matrix& setData(const Num2* data, unsigned dataLength); /** Set a row from an external array */ template Matrix& setRow( unsigned row, const Num2* data, unsigned dataLength); /** Set a column from an external array */ template Matrix& setColumn(unsigned col, const Num2* data, unsigned dataLength); /** // Tag this matrix as diagonal. This may improve the speed // of certain operations. Of course, this should be done only // if you know for sure that this matrix is, indeed, diagonal. // Works for square matrices only. // // This method can also be used to tag matrices as non-diagonal. */ Matrix& tagAsDiagonal(bool b=true); /** // Fill the main diagonal and set all other elements to zero. // This operation tags the matrix as diagonal. */ template Matrix& diagFill(const Num2* data, unsigned dataLen); //@{ /** A simple inspection of matrix properties */ inline unsigned nRows() const {return nrows_;} inline unsigned nColumns() const {return ncols_;} inline unsigned long length() const {return static_cast(nrows_)*ncols_;} inline Numeric* data() const {return data_;} inline bool isSquare() const {return data_ && ncols_ == nrows_;} bool isSymmetric() const; bool isAntiSymmetric() const; bool isDiagonal() const; bool isZero() const; //@} /** // This method resets the object to an unintialized state // Can be applied in order to force the assignment operators to work. */ Matrix& uninitialize(); /** // Check if this matrix has the same number of rows and columns // as the other one. "false" will be returned in case any one of // the two matrices (or both) is not initialized. */ bool isCompatible(const Matrix& other) const; /** // This method changes the object dimensions. // All data is lost in the process. */ Matrix& resize(unsigned nrows, unsigned ncols); /** This method sets all elements to 0 */ Matrix& zeroOut(); /** // This method sets all diagonal elements to 0. // It can only be used with square matrices. */ Matrix& clearMainDiagonal(); /** // This method sets all non-diagonal elements to 0. // This operation is only valid for square matrices. // It tags the matrix as diagonal. */ Matrix& makeDiagonal(); /** This method sets all elements to the given value */ Matrix& constFill(Numeric c); /** // Method for transposing this matrix. Uses std::swap for square // matrices but might allocate memory from the heap if the matrix // is not square. */ Matrix& Tthis(); //@{ /** Compare two matrices for equality */ bool operator==(const Matrix&) const; bool operator!=(const Matrix&) const; //@} /** // Non-const access to the data (works like 2-d array in C). // No bounds checking. */ Numeric* operator[](unsigned); /** // Const access to the data (works like 2-d array in C). // No bounds checking. */ const Numeric* operator[](unsigned) const; /** Data modification method with bounds checking */ Matrix& set(unsigned row, unsigned column, Numeric value); /** Access by value without bounds checking */ Numeric operator()(unsigned row, unsigned column) const; /** Access by value with bounds checking */ Numeric at(unsigned row, unsigned column) const; /** Sum of the values in the given row */ Numeric rowSum(unsigned row) const; /** Sum of the values in the given column */ Numeric columnSum(unsigned column) const; /** Remove row and/or column. Indices out of range are ignored */ Matrix removeRowAndColumn(unsigned row, unsigned column) const; /** Number of non-zero elements */ unsigned nonZeros() const; /** // Replace every n x m square block in the matrix by its sum. // The number of rows must be divisible by n. The number of // columns must be divisible by m. */ void coarseSum(unsigned n, unsigned m, Matrix* result) const; /** // Replace every n x m square block in the matrix by its average // The number of rows must be divisible by n. The number of // columns must be divisible by m. */ void coarseAverage(unsigned n, unsigned m, Matrix* result) const; /** Unary plus */ Matrix operator+() const; /** Unary minus */ Matrix operator-() const; //@{ /** Binary algebraic operation with matrix or scalar */ Matrix operator*(const Matrix& r) const; Matrix operator*(Numeric r) const; Matrix operator/(Numeric r) const; Matrix operator+(const Matrix& r) const; Matrix operator-(const Matrix& r) const; //@} /** Hadamard product (a.k.a. Schur product) */ Matrix hadamardProduct(const Matrix& r) const; /** Hadamard ratio. All elements of the denominator must be non-zero */ Matrix hadamardRatio(const Matrix& denominator) const; //@{ /** // Binary algebraic operation which writes its result into // an existing matrix potentially avoiding memory allocation */ void times(const Matrix& r, Matrix* result) const; void times(Numeric r, Matrix* result) const; void over(Numeric r, Matrix* result) const; void plus(const Matrix& r, Matrix* result) const; void minus(const Matrix& r, Matrix* result) const; //@} //@{ /** In-place algebraic operations with matrices and scalars */ Matrix& operator*=(Numeric r); Matrix& operator/=(Numeric r); Matrix& operator+=(const Matrix& r); Matrix& operator-=(const Matrix& r); //@} //@{ /** Multiplication by a vector represented by a pointer and array size */ template Matrix timesVector(const Num2* data, unsigned dataLen) const; template void timesVector(const Num2* data, unsigned dataLen, Num3* result, unsigned resultLen) const; //@} /** // Multiplication by a vector represented by a pointer and array size. // This function is useful when only one element of the result is needed // (or elements are needed one at a time). */ template Numeric timesVector(unsigned rowNumber, const Num2* data, unsigned dataLen) const; //@{ /** Multiplication by a row on the left */ template Matrix rowMultiply(const Num2* data, unsigned dataLen) const; template void rowMultiply(const Num2* data, unsigned dataLen, Num3* result, unsigned resultLen) const; //@} /** // Multiplication by a row on the left. This function is useful when // only one element of the result is needed (or elements are needed // one at a time). */ template Numeric rowMultiply(unsigned columnNumber, const Num2* data, unsigned dataLen) const; /** Bilinear form with a vector (v^T M v) */ template Numeric bilinear(const Num2* data, unsigned dataLen) const; /** Bilinear form with a matrix (r^T M r) */ Matrix bilinear(const Matrix& r) const; /** Bilinear form with a transposed matrix (r M r^T) */ Matrix bilinearT(const Matrix& r) const; /** // Solution of a linear system of equations M*x = rhs. // The matrix must be square. "false" is returned in case // the matrix is degenerate. */ bool solveLinearSystem(const Numeric* rhs, unsigned lenRhs, Numeric* solution) const; /** // Solution of linear systems of equations M*X = RHS. // The matrix M (this one) must be square. "false" is // returned in case this matrix is degenerate. */ bool solveLinearSystems(const Matrix& RHS, Matrix* X) const; /** // Solution of an underdetermined linear system of equations M*x = rhs // in the minimum norm sense (M is this matrix). The norm is defined by // sqrt(x^T*V^(-1)*x). The result is given by A*rhs, where A is the // right inverse of M (i.e., M*A = I) minimizing the norm. V plays // the role of the inverse Gram matrix. It must be symmetric and // positive-definite. // // This method can also be used for solving constrained least squares // problems. For example, minimization of (y - y0)^T*V^(-1)*(y - y0) // under the condition M*y = rhs can be reduced to solving an // underdetermined linear system using x = y - y0 and // M*x = rhs' = rhs - M*y0. The application code can call this function // assuming that V is the covariance matrix of y and using rhs' as the // "rhs" argument. Subsequently, y can be calculated as x + y0. // The covariance matrix of y can be subsequently calculated by // V_y = (I - P)*V*(I - P)^T, where P = A*M. Operator (I - P) actually // projects the x (i.e., y - y0) space onto the space of solutions of // the equation M*x = 0. // // "false" is returned in case of failure. */ bool underdeterminedLinearSystem(const Numeric* rhs, unsigned lenRhs, const Matrix& V, Numeric* solution, unsigned lenSolution, Numeric* resultNormSquared = 0, Matrix* A = 0) const; /** // Solution of a linear overdetermined system of equations M*x = rhs // in the least squares sense. The "lenRhs" parameter should be equal // to the number of matrix rows, and it should exceed "lenSolution" // ("lenSolution" should be equal to the number of columns). "false" // is returned in case of failure. */ bool linearLeastSquares(const Numeric* rhs, unsigned lenRhs, Numeric* solution, unsigned lenSolution) const; /** // Solution of a linear system of equations M*x = rhs1 in the least // squares sense, subject to the constraint B*x = rhs2. Number of // equations (i.e., dimensionality of rhs1) plus the number of // constraints (i.e., dimensionality of rhs2) should exceed the // dimensionality of x. "false" is returned in case of failure. If // the chi-square of the result is desired, make the "resultChiSquare" // argument point to some valid location. Same goes for the result // covariance matrix and all subsequent arguments which can be used // to extract various intermediate results: // // unconstrainedSolution -- Array to store the unconstrained solution. // Must have at least "lenSol" elements. // // unconstrainedCovmat -- Covariance matrix of the unconstrained // solution. // // projectionMatrix -- Projection matrix used to obtain a part // of the constrained solution from the // unconstrained one. // // A -- Matrix multiplied by rhs2 to get the second // part of the constrained solution. // // Note that requesting the result covariance matrix or any other // subsequent output switches the code to a different and slower // "textbook" algorithm instead of an optimized algorithm from Lapack. */ bool constrainedLeastSquares(const Numeric* rhs1, unsigned lenRhs1, const Matrix& B, const Numeric* rhs2, unsigned lenRhs2, Numeric* solution, unsigned lenSol, Numeric* resultChiSquare = 0, Matrix* resultCovarianceMatrix = 0, Numeric* unconstrainedSolution = 0, Matrix* unconstrainedCovmat = 0, Matrix* projectionMatrix = 0, Matrix* A = 0) const; /** // Weighted least squares problem. The inverse covariance matrix // must be symmetric and positive definite. If the chi-square of the // result is desired, make the "resultChiSquare" argument point to // some valid location. Same goes for the result covariance matrix. // "false" is returned in case of failure. */ bool weightedLeastSquares(const Numeric* rhs, unsigned lenRhs, const Matrix& inverseCovarianceMatrix, Numeric* solution, unsigned lenSolution, Numeric* resultChiSquare = 0, Matrix* resultCovarianceMatrix = 0) const; /** Return transposed matrix */ Matrix T() const; /** Return the product of this matrix transposed with this, M^T*M */ Matrix TtimesThis() const; /** Return the product of this matrix with its transpose, M*M^T */ Matrix timesT() const; /** Return the product of this matrix with the transposed argument */ Matrix timesT(const Matrix& r) const; /** Return the product of this matrix transposed with the argument */ Matrix Ttimes(const Matrix& r) const; /** // "directSum" method constructs a block-diagonal matrix with // this matrix and the argument matrix inside the diagonal blocks. // Block which contains this matrix is the top left one. */ Matrix directSum(const Matrix& added) const; /** Construct a symmetric matrix from this one (symmetrize) */ Matrix symmetrize() const; /** Construct an antisymmetric matrix from this one (antisymmetrize) */ Matrix antiSymmetrize() const; /** // Matrix outer product. The number of rows of the result equals // the product of the numbers of rows of this matrix and the argument // matrix. The number of columns of the result is the product of the // numbers of columns. */ Matrix outer(const Matrix& r) const; /** Minimum value for any row/column */ Numeric minValue() const; /** Maximum value for any row/column */ Numeric maxValue() const; /** Maximum absolute value for any row/column */ Numeric maxAbsValue() const; /** Row and column of the smallest matrix element */ std::pair argmin() const; /** Row and column of the largest matrix element */ std::pair argmax() const; /** // Row and column of the matrix element with the largest // absolute value */ std::pair argmaxAbs() const; /** Frobenius norm of this matrix */ double frobeniusNorm() const; //@{ /** Calculate the trace (matrix must be square) */ Numeric tr() const; inline Numeric sp() const {return tr();} //@} //@{ /** // Calculate the trace of a product of this matrix with another. // This works faster if only the trace is of interest but not the // product itself. */ Numeric productTr(const Matrix& r) const; inline Numeric productSp(const Matrix& r) const {return productTr(r);} //@} /** Calculate the determinant (matrix must be square) */ Numeric det() const; /** // Inverse of a real symmetric positive-definite matrix. Use this // function (it has better numerical precision than symInv()) if // you know that your matrix is symmetric and positive definite. // For example, a non-singular covariance matrix calculated for // some dataset has these properties. */ Matrix symPDInv() const; /** // Invert real symmetric positive-definite matrix using // eigendecomposition. Use this function if you know that your // matrix is supposed to be symmetric and positive-definite but // might not be due to numerical round-offs. The arguments of // this method are as folows: // // tol -- Eigenvalues whose ratios to the largest eigenvalue // are below this parameter will be either truncated // or extended (per value of the next argument). // This parameter must lie in (0, 1). // // extend -- If "true", the inverse of small eigenvalues will // be replaced by the inverse of the largest eigenvalue // divided by the tolerance parameter. If "false", // inverse of such eigenvalues will be set to 0. // // m -- LAPACK method to use for calculating the // eigendecomposition. */ Matrix symPDEigenInv(double tol, bool extend=true, EigenMethod m=EIGEN_SIMPLE) const; /** // Inverse of a real symmetric matrix. If the matrix is not symmetric, // its symmetrized version will be used internally. Use this function // (it has better numerical precision than inv()) if you know that // your matrix is symmetric. */ Matrix symInv() const; /** // Inverse of a general matrix (which must be square and non-singular) */ Matrix inv() const; /** + // Left pseudoinverse of a real, full rank matrix. This + // pseudoinverse can be calculated when the number of columns + // does not exceed the number of rows. + */ + Matrix leftInv() const; + + /** + // Right pseudoinverse of a real, full rank matrix. This + // pseudoinverse can be calculated when the number of rows + // does not exceed the number of columns. + */ + Matrix rightInv() const; + + /** // Eigenvalues and eigenvectors of a real tridiagonal symmetric // matrix. The "eigenvectors" pointer can be NULL if only eigenvalues // are needed. If it is not NULL, the _rows_ of that matrix will be // set to the computed eigenvectors. // // The "m" parameter specifies the LAPACK method to use for // calculating the eigendecomposition. */ void tdSymEigen(Numeric* eigenvalues, unsigned lenEigenvalues, Matrix* eigenvectors=0, EigenMethod m=EIGEN_RRR) const; /** // Eigenvalues and eigenvectors of a real symmetric matrix. // The "eigenvectors" pointer can be NULL if only eigenvalues // are needed. If it is not NULL, the _rows_ of that matrix // will be set to the computed eigenvectors. // // The "m" parameter specifies the LAPACK method to use for // calculating the eigendecomposition. */ void symEigen(Numeric* eigenvalues, unsigned lenEigenvalues, Matrix* eigenvectors=0, EigenMethod m=EIGEN_SIMPLE) const; /** // Eigenvalues and eigenvectors of a real general matrix. // Either "rightEigenvectors" or "leftEigenvectors" or both // can be NULL if they are not needed. // // If an eigenvalue is real, the corresponding row of the // output eigenvector matrix is set to the computed eigenvector. // If an eigenvalue is complex, all such eigenvalues come in // complex conjugate pairs. Corresponding eigenvectors make up // such a pair as well. The matrix row which corresponds to the // first eigenvalue in the pair is set to the real part of the pair // and the matrix row which corresponds to the second eigenvalue // in the pair is set to the imaginary part of the pair. For the // right eigenvectors, the first complex eigenvector in the pair // needs to be constructed by adding the imaginary part. The second // eigenvector has the imaginary part subtracted. For the left // eigenvectors the order is reversed: the first eigenvector needs // to be constructed by subtracting the imaginary part, and the // second by adding it. */ void genEigen(typename GeneralizedComplex::type *eigenvalues, unsigned lenEigenvalues, Matrix* rightEigenvectors = 0, Matrix* leftEigenvectors = 0) const; /** // Singular value decomposition of a matrix. The length of the // singular values buffer must be at least min(nrows, ncolumns). // // On exit, matrices U and V will contain singular vectors of // the current matrix, row-by-row. The decomposition is then // // *this = U^T * diag(singularValues, nrows, ncolumns) * V */ void svd(Numeric* singularValues, unsigned lenSingularValues, Matrix* U, Matrix* V, SvdMethod m = SVD_D_AND_C) const; /** // Calculate entropy-based effective rank. It is assumed that // the matrix is symmetric positive semidefinite. // // The "tol" parameter specifies the eigenvalue cutoff: // all eigenvalues equal to tol*maxEigenvalue or smaller // will be ignored. This parameter must be non-negative. // // The "m" parameter specifies the LAPACK method to use for // calculating eigenvalues. // // If "eigenSum" pointer is not NULL, *eigenSum will be filled // on exit with the sum of eigenvalues above the cutoff divided // by the largest eigenvalue. */ double symPSDefEffectiveRank(double tol = 0.0, EigenMethod m = EIGEN_D_AND_C, double* eigenSum = 0) const; /** // Calculate a function of this matrix. It is assumed that the // matrix is symmetric and real. The calculation is performed by // calculating the eignenbasis, evaluating the function of the // eigenvalues, and then transforming back to the original basis. */ template Matrix symFcn(const Functor& fcn) const; /** Power of a matrix. Matrix must be square. */ Matrix pow(unsigned degree) const; /** // The following function derives a correlation matrix // out of a covariance matrix. In the returned matrix // all diagonal elements will be set to 1. */ Matrix covarToCorr() const; /** // The following function derives a covariance matrix // out of a correlation matrix and of the covariances. */ Matrix corrToCovar(const double* variances, unsigned lenVariances) const; //@{ /** Method related to "geners" I/O */ inline gs::ClassId classId() const {return gs::ClassId(*this);} bool write(std::ostream& of) const; //@} static const char* classname(); static inline unsigned version() {return 1;} static void restore(const gs::ClassId& id, std::istream& in, Matrix* m); private: void initialize(unsigned nrows, unsigned ncols); void calcDiagonal(); void makeCoarse(unsigned n, unsigned m, Matrix* result) const; void invertDiagonal(Matrix* result) const; Numeric local_[Len]; Numeric* data_; unsigned nrows_; unsigned ncols_; unsigned len_; bool isDiagonal_; bool diagonalityKnown_; #ifdef SWIG public: inline std::vector mainDiagonal() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::mainDiagonal: uninitialized matrix"); const unsigned len = std::min(nrows_, ncols_); std::vector d(len); for (unsigned i=0; isetData(data, dataLen);} inline void setRow2(unsigned row, const Numeric* data, unsigned dataLen) {this->setRow(row, data, dataLen);} inline void setColumn2(unsigned col, const Numeric* data, unsigned dataLen) {this->setColumn(col, data, dataLen);} inline Matrix timesVector2(const Numeric* data, unsigned dataLen) const {return timesVector(data, dataLen);} inline Matrix rowMultiply2(const Numeric* data, unsigned dataLen) const {return rowMultiply(data, dataLen);} inline Numeric bilinear2(const Numeric* data, unsigned dataLen) const {return bilinear(data, dataLen);} inline std::vector solveLinearSystem2( const Numeric* data, unsigned dataLen) const { std::vector result(dataLen); if (!solveLinearSystem(data, dataLen, dataLen ? &result[0] : 0)) result.clear(); return result; } inline std::pair symPSDefEffectiveRank2( double tol = 0.0, EigenMethod m = EIGEN_D_AND_C) const { double eigenSum = 0; double erunk = this->symPSDefEffectiveRank(tol, m, &eigenSum); return std::pair(erunk, eigenSum); } inline std::vector symEigenValues() const { std::vector result(nrows_); symEigen(&result[0], nrows_, 0); return result; } inline Matrix rmul2(Numeric r) const {return *this * r;} // Perhaps, some day SWIG will be able to wrap the following. Currently, // its type deduction for GeneralizedComplex::type fails. // // inline std::vector::type> // genEigenValues() const // { // std::vector::type> r(nrows_); // genEigen(&r[0], nrows_); // return r; // } #endif // SWIG }; /** Utility for making square diagonal matrices from the given array */ template Matrix diag(const Numeric* data, unsigned dataLen); /** // Utility for making rectangular diagonal matrices from the given array. // The length of the array must be at least min(nrows, ncols). */ template Matrix diag(const Numeric* data, unsigned nrows, unsigned ncols); } #include template std::ostream& operator<<(std::ostream& os, const npstat::Matrix& m); template inline npstat::Matrix operator*(const Numeric& l, const npstat::Matrix& r) { return r*l; } #include "npstat/nm/Matrix.icc" #endif // NPSTAT_MATRIX_HH_