mirror/glm
1
0
Fork 0
mirror of https://github.com/g-truc/glm.git synced 2024-11-10 04:31:47 +00:00
Code Issues Packages Projects Releases Wiki Activity

Added further details on the comparison issue with covariance matrices on some VMs.

Browse Source 
Also corrected some code style guide, and changed `nullptr` to `GLM_NULLPTR` for better compatibility.
Tests are now executed in blocks of related tests, and only inbetween blocks the tests will exit.
This commit is contained in:
SGrottel 2021-05-10 21:32:01 +02:00
parent d0d7945141
commit a0ccbcc63d
2 changed files with 119 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

119
glm/gtx/pca.inl
View File

@ -29,15 +29,16 @@ namespace glm {
glm::mat<D, D, T, Q> m(0);
size_t cnt = 0;
for (I i = b; i != e; i++)
for(I i = b; i != e; i++)
{
vec<D, T, Q> 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<T>(v[x] * v[y]);
cnt++;
}
if (cnt > 0) m /= static_cast<T>(cnt);
if(cnt > 0)
m /= static_cast<T>(cnt);
return m;
}
@ -50,15 +51,16 @@ namespace glm {
glm::vec<D, T, Q> 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<T>(v[x] * v[y]);
cnt++;
}
if (cnt > 0) m /= static_cast<T>(cnt);
if(cnt > 0)
m /= static_cast<T>(cnt);
return m;
}
@ -77,12 +79,12 @@ namespace glm {
static const T epsilon = static_cast<T>(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<T>(1) + absb);
}
if (glm::equal<T>(absb, 0, epsilon)) return static_cast<T>(0);
if(glm::equal<T>(absb, 0, epsilon)) return static_cast<T>(0);
absa /= absb;
absa *= absa;
return absb * glm::sqrt(static_cast<T>(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<T>(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<T>(scale, 0, epsilon)) {
if(glm::equal<T>(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<T>(d[i - 1], 0, epsilon)) {
for (j = 1; j <= l; j++) {
if(!glm::equal<T>(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<T>(glm::abs(e[m - 1]) + dd, dd, epsilon)) break;
if(glm::equal<T>(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<T>(r, 0, epsilon)) {
if(glm::equal<T>(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<T>(r, 0, epsilon) && (i >= l)) continue;
if(glm::equal<T>(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

74
test/gtx/gtx_pca.cpp
View File

@ -32,8 +32,7 @@ bool matrixEpsilonEqual(glm::mat<D, D, T, Q> const& a, glm::mat<D, D, T, Q> cons
template<typename T>
T failReport(T line)
{
printf("I:Failed in line %d\n", static_cast<int>(line));
fprintf(stderr, "E:Failed in line %d\n", static_cast<int>(line));
fprintf(stderr, "Failed in line %d\n", static_cast<int>(line));
return line;
}
@ -211,12 +210,19 @@ namespace _1aga
template<glm::length_t D, typename T, glm::qualifier Q>
int checkCovarMat(glm::mat<D, D, T, Q> 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<T>(0.000001)))
{
fprintf(stderr, "E: %.15lf != %.15lf ; diff: %.20lf\n",
static_cast<double>(covarMat[y][x]),
static_cast<double>(expectedCovarData[x * 4 + y]),
static_cast<double>(covarMat[y][x] - expectedCovarData[x * 4 + y])
);
return failReport(__LINE__);
}
return 0;
}
@ -305,8 +311,8 @@ namespace _1aga
glm::vec<D, T, Q> const& evals,
glm::mat<D, D, T, Q> const& evecs)
{
const T* expectedEvals = nullptr;
const T* expectedEvecs = nullptr;
const T* expectedEvals = GLM_NULLPTR;
const T* expectedEvecs = GLM_NULLPTR;
getExpectedEigenvaluesEigenvectorsDataPtr<D, T>(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;
}