diff --git a/glm/gtx/pca.inl b/glm/gtx/pca.inl index c2eea458..d5a24b78 100644 --- a/glm/gtx/pca.inl +++ b/glm/gtx/pca.inl @@ -29,15 +29,16 @@ namespace glm { glm::mat m(0); size_t cnt = 0; - for (I i = b; i != e; i++) + for(I i = b; i != e; i++) { vec const& v = *i; - for (length_t x = 0; x < D; ++x) - for (length_t y = 0; y < D; ++y) + for(length_t x = 0; x < D; ++x) + for(length_t y = 0; y < D; ++y) m[x][y] += static_cast(v[x] * v[y]); cnt++; } - if (cnt > 0) m /= static_cast(cnt); + if(cnt > 0) + m /= static_cast(cnt); return m; } @@ -50,15 +51,16 @@ namespace glm { glm::vec v; size_t cnt = 0; - for (I i = b; i != e; i++) + for(I i = b; i != e; i++) { v = *i - c; - for (length_t x = 0; x < D; ++x) - for (length_t y = 0; y < D; ++y) + for(length_t x = 0; x < D; ++x) + for(length_t y = 0; y < D; ++y) m[x][y] += static_cast(v[x] * v[y]); cnt++; } - if (cnt > 0) m /= static_cast(cnt); + if(cnt > 0) + m /= static_cast(cnt); return m; } @@ -77,12 +79,12 @@ namespace glm { static const T epsilon = static_cast(0.0000001); T absa = glm::abs(a); T absb = glm::abs(b); - if (absa > absb) { + if(absa > absb) { absb /= absa; absb *= absb; return absa * glm::sqrt(static_cast(1) + absb); } - if (glm::equal(absb, 0, epsilon)) return static_cast(0); + if(glm::equal(absb, 0, epsilon)) return static_cast(0); absa /= absb; absa *= absa; return absb * glm::sqrt(static_cast(1) + absa); @@ -105,28 +107,33 @@ namespace glm { T d[D]; // diagonal elements T e[D]; // off-diagonal elements - for (length_t r = 0; r < D; r++) { - for (length_t c = 0; c < D; c++) { + for(length_t r = 0; r < D; r++) + for(length_t c = 0; c < D; c++) a[(r) * D + (c)] = covarMat[c][r]; - } - } // 1. Householder reduction. length_t l, k, j, i; T scale, hh, h, g, f; static const T epsilon = static_cast(0.0000001); - for (i = D; i >= 2; i--) { + for(i = D; i >= 2; i--) + { l = i - 1; h = scale = 0; - if (l > 1) { - for (k = 1; k <= l; k++) { + if(l > 1) + { + for(k = 1; k <= l; k++) + { scale += glm::abs(a[(i - 1) * D + (k - 1)]); } - if (glm::equal(scale, 0, epsilon)) { + if(glm::equal(scale, 0, epsilon)) + { e[i - 1] = a[(i - 1) * D + (l - 1)]; - } else { - for (k = 1; k <= l; k++) { + } + else + { + for(k = 1; k <= l; k++) + { a[(i - 1) * D + (k - 1)] /= scale; h += a[(i - 1) * D + (k - 1)] * a[(i - 1) * D + (k - 1)]; } @@ -136,50 +143,63 @@ namespace glm { h -= f * g; a[(i - 1) * D + (l - 1)] = f - g; f = 0; - for (j = 1; j <= l; j++) { + for(j = 1; j <= l; j++) + { a[(j - 1) * D + (i - 1)] = a[(i - 1) * D + (j - 1)] / h; g = 0; - for (k = 1; k <= j; k++) { + for(k = 1; k <= j; k++) + { g += a[(j - 1) * D + (k - 1)] * a[(i - 1) * D + (k - 1)]; } - for (k = j + 1; k <= l; k++) { + for(k = j + 1; k <= l; k++) + { g += a[(k - 1) * D + (j - 1)] * a[(i - 1) * D + (k - 1)]; } e[j - 1] = g / h; f += e[j - 1] * a[(i - 1) * D + (j - 1)]; } hh = f / (h + h); - for (j = 1; j <= l; j++) { + for(j = 1; j <= l; j++) + { f = a[(i - 1) * D + (j - 1)]; e[j - 1] = g = e[j - 1] - hh * f; - for (k = 1; k <= j; k++) { + for(k = 1; k <= j; k++) + { a[(j - 1) * D + (k - 1)] -= (f * e[k - 1] + g * a[(i - 1) * D + (k - 1)]); } } } - } else { + } + else + { e[i - 1] = a[(i - 1) * D + (l - 1)]; } d[i - 1] = h; } d[0] = 0; e[0] = 0; - for (i = 1; i <= D; i++) { + for(i = 1; i <= D; i++) + { l = i - 1; - if (!glm::equal(d[i - 1], 0, epsilon)) { - for (j = 1; j <= l; j++) { + if(!glm::equal(d[i - 1], 0, epsilon)) + { + for(j = 1; j <= l; j++) + { g = 0; - for (k = 1; k <= l; k++) { + for(k = 1; k <= l; k++) + { g += a[(i - 1) * D + (k - 1)] * a[(k - 1) * D + (j - 1)]; } - for (k = 1; k <= l; k++) { + for(k = 1; k <= l; k++) + { a[(k - 1) * D + (j - 1)] -= g * a[(k - 1) * D + (i - 1)]; } } } d[i - 1] = a[(i - 1) * D + (i - 1)]; a[(i - 1) * D + (i - 1)] = 1; - for (j = 1; j <= l; j++) { + for(j = 1; j <= l; j++) + { a[(j - 1) * D + (i - 1)] = a[(i - 1) * D + (j - 1)] = 0; } } @@ -189,20 +209,27 @@ namespace glm { T s, r, p, dd, c, b; const length_t MAX_ITER = 30; - for (i = 2; i <= D; i++) { + for(i = 2; i <= D; i++) + { e[i - 2] = e[i - 1]; } e[D - 1] = 0; - for (l = 1; l <= D; l++) { + for(l = 1; l <= D; l++) + { iter = 0; - do { - for (m = l; m <= D - 1; m++) { + do + { + for(m = l; m <= D - 1; m++) + { dd = glm::abs(d[m - 1]) + glm::abs(d[m - 1 + 1]); - if (glm::equal(glm::abs(e[m - 1]) + dd, dd, epsilon)) break; + if(glm::equal(glm::abs(e[m - 1]) + dd, dd, epsilon)) + break; } - if (m != l) { - if (iter++ == MAX_ITER) { + if(m != l) + { + if(iter++ == MAX_ITER) + { return 0; // Too many iterations in FindEigenvalues } g = (d[l - 1 + 1] - d[l - 1]) / (2 * e[l - 1]); @@ -210,11 +237,13 @@ namespace glm { g = d[m - 1] - d[l - 1] + e[l - 1] / (g + transferSign(r, g)); s = c = 1; p = 0; - for (i = m - 1; i >= l; i--) { + for(i = m - 1; i >= l; i--) + { f = s * e[i - 1]; b = c * e[i - 1]; e[i - 1 + 1] = r = pythag(f, g); - if (glm::equal(r, 0, epsilon)) { + if(glm::equal(r, 0, epsilon)) + { d[i - 1 + 1] -= p; e[m - 1] = 0; break; @@ -225,18 +254,20 @@ namespace glm { r = (d[i - 1] - g) * s + 2 * c * b; d[i - 1 + 1] = g + (p = s * r); g = c * r - b; - for (k = 1; k <= D; k++) { + for(k = 1; k <= D; k++) + { f = a[(k - 1) * D + (i - 1 + 1)]; a[(k - 1) * D + (i - 1 + 1)] = s * a[(k - 1) * D + (i - 1)] + c * f; a[(k - 1) * D + (i - 1)] = c * a[(k - 1) * D + (i - 1)] - s * f; } } - if (glm::equal(r, 0, epsilon) && (i >= l)) continue; + if(glm::equal(r, 0, epsilon) && (i >= l)) + continue; d[l - 1] -= p; e[l - 1] = g; e[m - 1] = 0; } - } while (m != l); + } while(m != l); } // 3. output diff --git a/test/gtx/gtx_pca.cpp b/test/gtx/gtx_pca.cpp index 8a6b56b7..6b741964 100644 --- a/test/gtx/gtx_pca.cpp +++ b/test/gtx/gtx_pca.cpp @@ -32,8 +32,7 @@ bool matrixEpsilonEqual(glm::mat const& a, glm::mat cons template T failReport(T line) { - printf("I:Failed in line %d\n", static_cast(line)); - fprintf(stderr, "E:Failed in line %d\n", static_cast(line)); + fprintf(stderr, "Failed in line %d\n", static_cast(line)); return line; } @@ -211,12 +210,19 @@ namespace _1aga template int checkCovarMat(glm::mat const& covarMat) { - const T* expectedCovarData = nullptr; + const T* expectedCovarData = GLM_NULLPTR; getExpectedCovarDataPtr(expectedCovarData); for(glm::length_t x = 0; x < D; ++x) for(glm::length_t y = 0; y < D; ++y) if(!glm::equal(covarMat[y][x], expectedCovarData[x * 4 + y], static_cast(0.000001))) + { + fprintf(stderr, "E: %.15lf != %.15lf ; diff: %.20lf\n", + static_cast(covarMat[y][x]), + static_cast(expectedCovarData[x * 4 + y]), + static_cast(covarMat[y][x] - expectedCovarData[x * 4 + y]) + ); return failReport(__LINE__); + } return 0; } @@ -305,8 +311,8 @@ namespace _1aga glm::vec const& evals, glm::mat const& evecs) { - const T* expectedEvals = nullptr; - const T* expectedEvecs = nullptr; + const T* expectedEvals = GLM_NULLPTR; + const T* expectedEvecs = GLM_NULLPTR; getExpectedEigenvaluesEigenvectorsDataPtr(expectedEvals, expectedEvecs); for(int i = 0; i < D; ++i) @@ -441,8 +447,7 @@ int testCovar(glm::length_t dataSize, unsigned int randomEngineSeed) mat covarMat = glm::computeCovarianceMatrix(testData.data(), testData.size(), center); if(_1aga::checkCovarMat(covarMat)) { - fprintf(stdout, "I:Reconstructed covarMat:\n%s\n", glm::to_string(covarMat).c_str()); - fprintf(stderr, "E:Reconstructed covarMat:\n%s\n", glm::to_string(covarMat).c_str()); + fprintf(stderr, "Reconstructed covarMat:\n%s\n", glm::to_string(covarMat).c_str()); return failReport(__LINE__); } @@ -636,6 +641,7 @@ int rndTest(unsigned int randomEngineSeed) int main() { + int error(0); // A small smoke test to fail early with most problems if(smokeTest()) @@ -643,54 +649,62 @@ int main() // test sorting utility. if(testEigenvalueSort<2, float, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvalueSort<2, double, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvalueSort<3, float, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvalueSort<3, double, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvalueSort<4, float, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvalueSort<4, double, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); + if (error != 0) + return error; // Note: the random engine uses a fixed seed to create consistent and reproducible test data // test covariance matrix computation from different data sources if(testCovar<2, float, glm::defaultp>(100, 12345) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testCovar<2, double, glm::defaultp>(100, 42) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testCovar<3, float, glm::defaultp>(100, 2021) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testCovar<3, double, glm::defaultp>(100, 815) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testCovar<4, float, glm::defaultp>(100, 3141) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testCovar<4, double, glm::defaultp>(100, 174) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); + if (error != 0) + return error; // test PCA eigen vector reconstruction if(testEigenvectors<2, float, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvectors<2, double, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvectors<3, float, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(testEigenvectors<3, double, glm::defaultp>() != 0) - return failReport(__LINE__); - if (testEigenvectors<4, float, glm::defaultp>() != 0) - return failReport(__LINE__); - if (testEigenvectors<4, double, glm::defaultp>() != 0) - return failReport(__LINE__); + error = failReport(__LINE__); + if(testEigenvectors<4, float, glm::defaultp>() != 0) + error = failReport(__LINE__); + if(testEigenvectors<4, double, glm::defaultp>() != 0) + error = failReport(__LINE__); + if(error != 0) + return error; // Final tests with randomized data #if GLM_HAS_CXX11_STL == 1 if(rndTest(12345) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); if(rndTest(42) != 0) - return failReport(__LINE__); + error = failReport(__LINE__); + if (error != 0) + return error; #endif // GLM_HAS_CXX11_STL == 1 - return 0; + return error; }