anulax
/
glm
1
0
Fork 0
Code Issues Pull Requests Packages project_board Releases Wiki Activity

Merge pull request #1077 from sgrottel/gtx-pca-test-rework

Browse Source 
Reworked test for GTX PCA #1077 #1071
master
Christophe ago%!(EXTRA string=4 years) committed by GitHub
parent 66062497b1 104c1214ee
commit 8efb7408b5
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 175 additions and 187 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 362
      test/gtx/gtx_pca.cpp

362
test/gtx/gtx_pca.cpp
Unescape View File

@ -17,33 +17,47 @@ GLM_INLINE GLM_CONSTEXPR float myEpsilon<float>() { return 0.00001f; }
template<>
GLM_INLINE GLM_CONSTEXPR double myEpsilon<double>() { return 0.000001; }
template<typename T>
T myEpsilon2();
template<>
GLM_INLINE GLM_CONSTEXPR float myEpsilon2<float>() { return 0.01f; }
template<>
GLM_INLINE GLM_CONSTEXPR double myEpsilon2<double>() { return 0.000001; }
template<glm::length_t D, typename T, glm::qualifier Q>
bool vectorEpsilonEqual(glm::vec<D, T, Q> const& a, glm::vec<D, T, Q> const& b)
bool vectorEpsilonEqual(glm::vec<D, T, Q> const& a, glm::vec<D, T, Q> const& b, T epsilon)
{
for (int c = 0; c < D; ++c)
if (!glm::epsilonEqual(a[c], b[c], myEpsilon<T>()))
if (!glm::epsilonEqual(a[c], b[c], epsilon))
{
fprintf(stderr, "failing vectorEpsilonEqual: [%d] %lf != %lf (~%lf)\n",
c,
static_cast<double>(a[c]),
static_cast<double>(b[c]),
static_cast<double>(epsilon)
);
return false;
}
return true;
}
template<glm::length_t D, typename T, glm::qualifier Q>
bool matrixEpsilonEqual(glm::mat<D, D, T, Q> const& a, glm::mat<D, D, T, Q> const& b)
bool matrixEpsilonEqual(glm::mat<D, D, T, Q> const& a, glm::mat<D, D, T, Q> const& b, T epsilon)
{
for (int c = 0; c < D; ++c)
for (int r = 0; r < D; ++r)
if (!glm::epsilonEqual(a[c][r], b[c][r], myEpsilon<T>()))
if (!glm::epsilonEqual(a[c][r], b[c][r], epsilon))
{
fprintf(stderr, "failing vectorEpsilonEqual: [%d][%d] %lf != %lf (~%lf)\n",
c, r,
static_cast<double>(a[c][r]),
static_cast<double>(b[c][r]),
static_cast<double>(epsilon)
);
return false;
}
return true;
}
template<typename T>
GLM_INLINE bool sameSign(T const& a, T const& b)
{
return ((a >= 0) && (b >= 0)) || ((a < 0) && (b < 0));
}
template<typename T>
T failReport(T line)
{
@ -200,153 +214,123 @@ namespace _1aga
outTestData[i][d] = static_cast<typename vec::value_type>(_1aga[i * 4 + d]);
}
void getExpectedCovarDataPtr(const double*& ptr)
{
static const double _1agaCovar4x4d[] = {
9.624340680272107, -0.000066573696146, -4.293213765684049, 0.018793741874528,
-0.000066573696146, 9.624439378684805, 5.351138726379443, -0.115692591458806,
-4.293213765684049, 5.351138726379443, 35.628485496346691, 0.908742392542202,
0.018793741874528, -0.115692591458806, 0.908742392542202, 1.097059718568909
};
ptr = _1agaCovar4x4d;
}
void getExpectedCovarDataPtr(const float*& ptr)
{
// note: the value difference to `_1agaCovar4x4d` is due to the numeric error propagation during computation of the covariance matrix.
static const float _1agaCovar4x4f[] = {
9.624336242675781f, -0.000066711785621f, -4.293214797973633f, 0.018793795257807f,
-0.000066711785621f, 9.624438285827637f, 5.351140022277832f, -0.115692682564259f,
-4.293214797973633f, 5.351140022277832f, 35.628479003906250f, 0.908742427825928f,
0.018793795257807f, -0.115692682564259f, 0.908742427825928f, 1.097059369087219f
};
ptr = _1agaCovar4x4f;
}
// All reference values computed separately using symbolic precision
// https://github.com/sgrottel/exp-pca-precision
// This applies to all functions named: `_1aga::expected*()`
template<glm::length_t D, typename T, glm::qualifier Q>
int checkCovarMat(glm::mat<D, D, T, Q> const& covarMat)
GLM_INLINE glm::dmat4 const& expectedCovarData()
{
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], myEpsilon<T>()))
{
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;
static const glm::dmat4 covar4x4d(
9.62434068027210898322, -0.00006657369614512471, -4.29321376568405099761, 0.01879374187452758846,
-0.00006657369614512471, 9.62443937868480681175, 5.35113872637944076871, -0.11569259145880574080,
-4.29321376568405099761, 5.35113872637944076871, 35.62848549634668415820, 0.90874239254220201545,
0.01879374187452758846, -0.11569259145880574080, 0.90874239254220201545, 1.09705971856890904803
);
return covar4x4d;
}
template<glm::length_t D, typename T> void getExpectedEigenvaluesEigenvectorsDataPtr(const T*& evals, const T*& evecs);
template<> void getExpectedEigenvaluesEigenvectorsDataPtr<2, float>(const float*& evals, const float*& evecs)
template<glm::length_t D>
GLM_INLINE glm::vec<D, double, glm::defaultp> const& expectedEigenvalues();
template<>
GLM_INLINE glm::dvec2 const& expectedEigenvalues<2>()
{
static const float expectedEvals[] = {
9.624471664428711f, 9.624302864074707f
};
static const float expectedEvecs[] = {
-0.443000972270966f, 0.896521151065826f,
0.896521151065826f, 0.443000972270966f
};
evals = expectedEvals;
evecs = expectedEvecs;
static const glm::dvec2 evals2(
9.62447289926297399961763301774251330057894539467032275382255,
9.62430715969394210015560961264297422776572580714373620309355
);
return evals2;
}
template<> void getExpectedEigenvaluesEigenvectorsDataPtr<2, double>(const double*& evals, const double*& evecs)
template<>
GLM_INLINE glm::dvec3 const& expectedEigenvalues<3>()
{
static const double expectedEvals[] = {
9.624472899262972, 9.624307159693940
};
static const double expectedEvecs[] = {
-0.449720461624363, 0.893169360421846,
0.893169360421846, 0.449720461624363
};
evals = expectedEvals;
evecs = expectedEvecs;
static const glm::dvec3 evals3(
37.3274494274683425233695502581182052836449738530676689472257,
9.62431434161498823505729817436585077939509766554969096873168,
7.92550178622027216422369326567668971675332732240052872097887
);
return evals3;
}
template<> void getExpectedEigenvaluesEigenvectorsDataPtr<3, float>(const float*& evals, const float*& evecs)
template<>
GLM_INLINE glm::dvec4 const& expectedEigenvalues<4>()
{
static const float expectedEvals[] = {
37.327442169189453f, 9.624311447143555f, 7.925499439239502f
};
static const float expectedEvecs[] = {
-0.150428697466850f, 0.187497511506081f, 0.970678031444550f,
0.779980957508087f, 0.625803351402283f, -0.000005212802080f,
0.607454538345337f, -0.757109522819519f, 0.240383237600327f
};
evals = expectedEvals;
evecs = expectedEvecs;
static const glm::dvec4 evals4(
37.3477389918792213596879452204499702406947817221901007885630,
9.62470688921105696017807313860277172063600080413412567999700,
7.94017075281634999342344275928070533134615133171969063657713,
1.06170863996588365446060186982477896078741484440002343404155
);
return evals4;
}
template<> void getExpectedEigenvaluesEigenvectorsDataPtr<3, double>(const double*& evals, const double*& evecs)
template<glm::length_t D>
GLM_INLINE glm::mat<D, D, double, glm::defaultp> const& expectedEigenvectors();
template<>
GLM_INLINE glm::dmat2 const& expectedEigenvectors<2>()
{
static const double expectedEvals[] = {
37.327449427468345, 9.624314341614987, 7.925501786220276
};
static const double expectedEvecs[] = {
-0.150428640509585, 0.187497426513576, 0.970678082149394,
0.779981605126846, 0.625802441381904, -0.000004919018357,
0.607453635908278, -0.757110308615089, 0.240383154173870
};
evals = expectedEvals;
evecs = expectedEvecs;
static const glm::dmat2 evecs2(
glm::dvec2(
-0.503510847492551904906870957742619139443409162857537237123308,
1
),
glm::dvec2(
1.98605453086051402895741763848787613048533838388005162794043,
1
)
);
return evecs2;
}
template<> void getExpectedEigenvaluesEigenvectorsDataPtr<4, float>(const float*& evals, const float*& evecs)
template<>
GLM_INLINE glm::dmat3 const& expectedEigenvectors<3>()
{
static const float expectedEvals[] = {
37.347740173339844f, 9.624703407287598f, 7.940164566040039f, 1.061712265014648f
};
static const float expectedEvecs[] = {
-0.150269940495491f, 0.187220811843872f, 0.970467865467072f, 0.023652425035834f,
0.779159665107727f, 0.626788496971130f, -0.000105984276161f, -0.006797631736845f,
0.608242213726044f, -0.755563497543335f, 0.238818943500519f, 0.046158745884895f,
-0.019251370802522f, 0.034755907952785f, -0.034024771302938f, 0.998630762100220f,
};
evals = expectedEvals;
evecs = expectedEvecs;
static const glm::dmat3 evecs3(
glm::dvec3(
-0.154972738414395866005286433008304444294405085038689821864654,
0.193161285869815165989799191097521722568079378840201629578695,
1
),
glm::dvec3(
-158565.112775416943154745839952575022429933119522746586149868,
-127221.506282351944358932458687410410814983610301927832439675,
1
),
glm::dvec3(
2.52702248596556806145700361724323960543858113426446460406536,
-3.14959802931313870497377546974185300816008580801457419079412,
1
)
);
return evecs3;
}
template<> void getExpectedEigenvaluesEigenvectorsDataPtr<4, double>(const double*& evals, const double*& evecs)
template<>
GLM_INLINE glm::dmat4 const& expectedEigenvectors<4>()
{
static const double expectedEvals[] = {
37.347738991879226, 9.624706889211053, 7.940170752816341, 1.061708639965897
};
static const double expectedEvecs[] = {
-0.150269954805403, 0.187220917596058, 0.970467838469868, 0.023652551509145,
0.779159831346545, 0.626788431871120, -0.000105940250315, -0.006797622027466,
0.608241962267880, -0.755563776664248, 0.238818902950296, 0.046158707986616,
-0.019251317755512, 0.034755849578017, -0.034024915369495, 0.998630924225204,
};
evals = expectedEvals;
evecs = expectedEvecs;
}
template<glm::length_t D, typename T, glm::qualifier Q>
int checkEigenvaluesEigenvectors(
glm::vec<D, T, Q> const& evals,
glm::mat<D, D, T, Q> const& evecs)
{
const T* expectedEvals = GLM_NULLPTR;
const T* expectedEvecs = GLM_NULLPTR;
getExpectedEigenvaluesEigenvectorsDataPtr<D, T>(expectedEvals, expectedEvecs);
for(int i = 0; i < D; ++i)
if(!glm::equal(evals[i], expectedEvals[i], myEpsilon<T>()))
return failReport(__LINE__);
for (int i = 0; i < D; ++i)
for (int d = 0; d < D; ++d)
if (!glm::equal(evecs[i][d], expectedEvecs[i * D + d], myEpsilon2<T>()))
{
fprintf(stderr, "E: %.15lf != %.15lf ; diff: %.20lf\n",
static_cast<double>(evecs[i][d]),
static_cast<double>(expectedEvecs[i * D + d]),
static_cast<double>(evecs[i][d] - expectedEvecs[i * D + d])
);
return failReport(__LINE__);
}
return 0;
static const glm::dmat4 evecs4(
glm::dvec4(
-6.35322390281037045217295803597357821705371650876122113027264,
7.91546394153385394517767054617789939529794642646629201212056,
41.0301543819240679808549819457450130787045236815736490549663,
1
),
glm::dvec4(
-114.622418941087829756565311692197154422302604224781253861297,
-92.2070185807065289900871215218752013659402949497379896153118,
0.0155846091025912430932734548933329458404665760587569100867246,
1
),
glm::dvec4(
13.1771887761559019483954743159026938257325190511642952175789,
-16.3688257459634877666638419310116970616615816436949741766895,
5.17386502341472097227408249233288958059579189051394773143190,
1
),
glm::dvec4(
-0.0192777078948229800494895064532553117703859768210647632969276,
0.0348034950916108873629241563077465542944938906271231198634442,
-0.0340715609308469289267379681032545422644143611273049912226126,
1
)
);
return evecs4;
}
} // namespace _1aga
@ -425,9 +409,9 @@ int testEigenvalueSort()
};
// initial sanity check
if(!vectorEpsilonEqual(refVal, refVal))
if(!vectorEpsilonEqual(refVal, refVal, myEpsilon<T>()))
return failReport(__LINE__);
if(!matrixEpsilonEqual(refVec, refVec))
if(!matrixEpsilonEqual(refVec, refVec, myEpsilon<T>()))
return failReport(__LINE__);
// Exhaustive search through all permutations
@ -443,9 +427,9 @@ int testEigenvalueSort()
glm::sortEigenvalues(testVal, testVec);
if (!vectorEpsilonEqual(testVal, refVal))
if (!vectorEpsilonEqual(testVal, refVal, myEpsilon<T>()))
return failReport(__LINE__);
if (!matrixEpsilonEqual(testVec, refVec))
if (!matrixEpsilonEqual(testVec, refVec, myEpsilon<T>()))
return failReport(__LINE__);
}
@ -473,7 +457,7 @@ int testCovar(
vec center = computeCenter(testData);
mat covarMat = glm::computeCovarianceMatrix(testData.data(), testData.size(), center);
if(_1aga::checkCovarMat(covarMat))
if(!matrixEpsilonEqual(covarMat, mat(_1aga::expectedCovarData()), myEpsilon<T>()))
{
fprintf(stderr, "Reconstructed covarMat:\n%s\n", glm::to_string(covarMat).c_str());
return failReport(__LINE__);
@ -505,27 +489,27 @@ int testCovar(
mat c3 = glm::computeCovarianceMatrix(testData.data(), testData.size(), center);
mat c4 = glm::computeCovarianceMatrix<D, T, Q>(testData.rbegin(), testData.rend(), center);
if(!matrixEpsilonEqual(c1, c2))
if(!matrixEpsilonEqual(c1, c2, myEpsilon<T>()))
return failReport(__LINE__);
if(!matrixEpsilonEqual(c1, c3))
if(!matrixEpsilonEqual(c1, c3, myEpsilon<T>()))
return failReport(__LINE__);
if(!matrixEpsilonEqual(c1, c4))
if(!matrixEpsilonEqual(c1, c4, myEpsilon<T>()))
return failReport(__LINE__);
#endif // GLM_HAS_CXX11_STL == 1
return 0;
}
// Computes eigenvalues and eigenvectors from well-known covariance matrix
template<glm::length_t D, typename T, glm::qualifier Q>
int testEigenvectors()
int testEigenvectors(T epsilon)
{
typedef glm::vec<D, T, Q> vec;
typedef glm::mat<D, D, T, Q> mat;
// test expected result with fixed data set
std::vector<vec> testData;
_1aga::fillTestData(testData);
vec center = computeCenter(testData);
mat covarMat = glm::computeCovarianceMatrix(testData.data(), testData.size(), center);
mat covarMat(_1aga::expectedCovarData());
vec eigenvalues;
mat eigenvectors;
unsigned int c = glm::findEigenvaluesSymReal(covarMat, eigenvalues, eigenvectors);
@ -533,16 +517,25 @@ int testEigenvectors()
return failReport(__LINE__);
glm::sortEigenvalues(eigenvalues, eigenvectors);
if(_1aga::checkEigenvaluesEigenvectors(eigenvalues, eigenvectors) != 0)
if (!vectorEpsilonEqual(eigenvalues, vec(_1aga::expectedEigenvalues<D>()), epsilon))
return failReport(__LINE__);
for (int i = 0; i < D; ++i)
{
vec act = glm::normalize(eigenvectors[i]);
vec exp = glm::normalize(_1aga::expectedEigenvectors<D>()[i]);
if (!sameSign(act[0], exp[0])) exp = -exp;
if (!vectorEpsilonEqual(act, exp, epsilon))
return failReport(__LINE__);
}
return 0;
}
/// A simple small smoke test:
/// - a uniformly sampled block
/// - reconstruct main axes
/// - check order of eigenvalues equals order of extends of block in direction of main axes
// A simple small smoke test:
// - a uniformly sampled block
// - reconstruct main axes
// - check order of eigenvalues equals order of extends of block in direction of main axes
int smokeTest()
{
using glm::vec3;
@ -579,11 +572,11 @@ int smokeTest()
std::swap(eVec[1], eVec[2]);
}
if(!vectorEpsilonEqual(glm::abs(eVec[0]), vec3(0, 1, 0)))
if(!vectorEpsilonEqual(glm::abs(eVec[0]), vec3(0, 1, 0), myEpsilon<float>()))
return failReport(__LINE__);
if(!vectorEpsilonEqual(glm::abs(eVec[1]), vec3(1, 0, 0)))
if(!vectorEpsilonEqual(glm::abs(eVec[1]), vec3(1, 0, 0), myEpsilon<float>()))
return failReport(__LINE__);
if(!vectorEpsilonEqual(glm::abs(eVec[2]), vec3(0, 0, 1)))
if(!vectorEpsilonEqual(glm::abs(eVec[2]), vec3(0, 0, 1), myEpsilon<float>()))
return failReport(__LINE__);
return 0;
@ -615,14 +608,9 @@ int rndTest(unsigned int randomEngineSeed)
glm::dvec3 z = glm::normalize(glm::cross(x, y));
y = glm::normalize(glm::cross(z, x));
//printf("\n");
//printf("x: %.10lf, %.10lf, %.10lf\n", x.x, x.y, x.z);
//printf("y: %.10lf, %.10lf, %.10lf\n", y.x, y.y, y.z);
//printf("z: %.10lf, %.10lf, %.10lf\n", z.x, z.y, z.z);
// generate input point data
std::vector<glm::dvec3> ptData;
static const int patters[] = {
static const int pattern[] = {
8, 0, 0,
4, 1, 2,
0, 2, 0,
@ -635,9 +623,9 @@ int rndTest(unsigned int randomEngineSeed)
for(int zs = 1; zs >= -1; zs -= 2)
ptData.push_back(
offset
+ x * static_cast<double>(patters[p * 3 + 0] * xs)
+ y * static_cast<double>(patters[p * 3 + 1] * ys)
+ z * static_cast<double>(patters[p * 3 + 2] * zs));
+ x * static_cast<double>(pattern[p * 3 + 0] * xs)
+ y * static_cast<double>(pattern[p * 3 + 1] * ys)
+ z * static_cast<double>(pattern[p * 3 + 2] * zs));
// perform PCA:
glm::dvec3 center = computeCenter(ptData);
@ -649,16 +637,14 @@ int rndTest(unsigned int randomEngineSeed)
return failReport(__LINE__);
glm::sortEigenvalues(evals, evecs);
//printf("\n");
//printf("evec0: %.10lf, %.10lf, %.10lf\n", evecs[0].x, evecs[0].y, evecs[0].z);
//printf("evec2: %.10lf, %.10lf, %.10lf\n", evecs[2].x, evecs[2].y, evecs[2].z);
//printf("evec1: %.10lf, %.10lf, %.10lf\n", evecs[1].x, evecs[1].y, evecs[1].z);
if(glm::length(glm::abs(x) - glm::abs(evecs[0])) > myEpsilon<double>())
if (!sameSign(evecs[0][0], x[0])) evecs[0] = -evecs[0];
if(!vectorEpsilonEqual(x, evecs[0], myEpsilon<double>()))
return failReport(__LINE__);
if(glm::length(glm::abs(y) - glm::abs(evecs[2])) > myEpsilon<double>())
if (!sameSign(evecs[2][0], y[0])) evecs[2] = -evecs[2];
if (!vectorEpsilonEqual(y, evecs[2], myEpsilon<double>()))
return failReport(__LINE__);
if(glm::length(glm::abs(z) - glm::abs(evecs[1])) > myEpsilon<double>())
if (!sameSign(evecs[1][0], z[0])) evecs[1] = -evecs[1];
if (!vectorEpsilonEqual(z, evecs[1], myEpsilon<double>()))
return failReport(__LINE__);
return 0;
@ -707,17 +693,19 @@ int main()
return error;
// test PCA eigen vector reconstruction
if(testEigenvectors<2, float, glm::defaultp>() != 0)
// Expected epsilon precision evaluated separately:
// https://github.com/sgrottel/exp-pca-precision
if(testEigenvectors<2, float, glm::defaultp>(0.002f) != 0)
error = failReport(__LINE__);
if(testEigenvectors<2, double, glm::defaultp>() != 0)
if(testEigenvectors<2, double, glm::defaultp>(0.00000000001) != 0)
error = failReport(__LINE__);
if(testEigenvectors<3, float, glm::defaultp>() != 0)
if(testEigenvectors<3, float, glm::defaultp>(0.00001f) != 0)
error = failReport(__LINE__);
if(testEigenvectors<3, double, glm::defaultp>() != 0)
if(testEigenvectors<3, double, glm::defaultp>(0.0000000001) != 0)
error = failReport(__LINE__);
if(testEigenvectors<4, float, glm::defaultp>() != 0)
if(testEigenvectors<4, float, glm::defaultp>(0.0001f) != 0)
error = failReport(__LINE__);
if(testEigenvectors<4, double, glm::defaultp>() != 0)
if(testEigenvectors<4, double, glm::defaultp>(0.0000001) != 0)
error = failReport(__LINE__);
if(error != 0)
return error;

Write Preview
Loading…
Cancel
Save