From a6e8434c6650203db29af0c8d12d8814d57e10b8 Mon Sep 17 00:00:00 2001
From: Guillaume Bittoun <guillaume.bittoun@gmail.com>
Date: Mon, 8 Dec 2014 13:15:47 +0000
Subject: [PATCH] simplex_interpolator: adding tests for simplex base matrix
 check.

---
 src/lol/math/simplex_interpolator.h |  13 +++-
 src/t/math/simplex_interpolator.cpp | 104 ++++++++++++++++++++++++++--
 2 files changed, 110 insertions(+), 7 deletions(-)

diff --git a/src/lol/math/simplex_interpolator.h b/src/lol/math/simplex_interpolator.h
index 7dee0d98..ce38a0fa 100644
--- a/src/lol/math/simplex_interpolator.h
+++ b/src/lol/math/simplex_interpolator.h
@@ -21,12 +21,23 @@ class simplex_interpolator
 {
 public:
 
+    simplex_interpolator() :
+        m_samples()
+    {
+        this->InitBase();
+    }
+
     simplex_interpolator(arraynd<N, T> const & samples) :
         m_samples(samples)
     {
         this->InitBase();
     }
 
+    inline void SetSamples(arraynd<N, T> const & samples)
+    {
+        this->m_samples = samples;
+    }
+
     inline arraynd<N, T> const & GetSamples() const
     {
         return this->m_samples;
@@ -132,7 +143,7 @@ protected:
                 | b b a b … |                 | d d c d … |
                 | …         |                 | …         |
 
-            where a and b are computed below ↴
+            where a, b, c, d are computed below ↴
         */
 
         float b = (1.f - lol::sqrt(N + 1.f)) / lol::sqrt((float)N * N * N);
diff --git a/src/t/math/simplex_interpolator.cpp b/src/t/math/simplex_interpolator.cpp
index 3a4db250..4884ed88 100644
--- a/src/t/math/simplex_interpolator.cpp
+++ b/src/t/math/simplex_interpolator.cpp
@@ -24,12 +24,16 @@ class test_interpolator : public simplex_interpolator<N, T>
 {
 public:
 
+    test_interpolator() :
+        simplex_interpolator<N, T>()
+    {
+    }
+
     test_interpolator(arraynd<N, T> const & samples) :
         simplex_interpolator<N, T>(samples)
     {
     }
 
-
     void DumpMatrix()
     {
         std::cout << std::endl;
@@ -55,6 +59,35 @@ public:
 
         std::cout << std::endl;
     }
+
+    void DumpMatrix(float result[N][N])
+    {
+        for (int i = 0; i < N; ++i)
+        {
+            for (int j = 0; j < N; ++j)
+            {
+                result[i][j] = this->m_base[i][j];
+            }
+        }
+    }
+
+    void DumpCheckInverse(float result[N][N])
+    {
+        for (int i = 0; i < N; ++i)
+            for (int j = 0; j < N; ++j)
+                result[i][j] = 0;
+
+        for (int i = 0; i < N; ++i)
+        {
+            for (int j = 0; j < N; ++j)
+            {
+                for (int k = 0; k < N; ++k)
+                {
+                    result[i][j] += this->m_base[i][k] * this->m_base_inverse[k][j];
+                }
+            }
+        }
+    }
 };
 
 lolunit_declare_fixture(SimplexInterpolatorTest)
@@ -69,13 +102,72 @@ lolunit_declare_fixture(SimplexInterpolatorTest)
         test_interpolator<2, vec2> c({{ vec2(0) }});
     }
 
-    lolunit_declare_test(CoordinateMatrixTest)
+    template<int N>
+    void check_base_matrix()
     {
-        test_interpolator<2> s2({{1.f}});
-        s2.DumpMatrix();
+        test_interpolator<N> s;
 
-        test_interpolator<3> s3(arraynd<3, float>({{{1.f}}}));
-        s3.DumpMatrix();
+        float result[N][N];
+        s.DumpCheckInverse(result);
+
+        // Check base matrix inverse
+        for (int i = 0; i < N; ++i)
+        {
+            for (int j = 0; j < N; ++j)
+            {
+                if (i == j)
+                    lolunit_assert_doubles_equal(1, result[i][j], 1e-6);
+                else
+                    lolunit_assert_doubles_equal(0, result[i][j], 1e-6);
+            }
+        }
+
+        s.DumpMatrix(result);
+
+        // Check base vectors’ norms
+        for (int i = 0; i < N; ++i)
+        {
+            float norm = 0;
+
+            for (int j = 0; j < N; ++j)
+            {
+                norm += result[i][j] * result[i][j];
+            }
+
+            lolunit_assert_doubles_equal(1, norm, 1e-6);
+        }
+
+        // Check result of vectors x, y z => must have norm = 1
+        float vec_result[N];
+        for (int i = 0; i < N; ++i)
+        {
+            vec_result[i] = 0;
+            for (int j = 0; j < N; ++j)
+            {
+                vec_result[i] += result[i][j];
+            }
+        }
+
+        float norm = 0;
+        for (int i = 0 ; i < N ; ++i)
+        {
+            norm += vec_result[i] * vec_result[i];
+        }
+        lolunit_assert_doubles_equal(1, norm, 1e-6);
+    }
+
+    lolunit_declare_test(CoordinateMatrixTest)
+    {
+        check_base_matrix<1>();
+        check_base_matrix<2>();
+        check_base_matrix<3>();
+        check_base_matrix<4>();
+        check_base_matrix<5>();
+        check_base_matrix<6>();
+        check_base_matrix<7>();
+        check_base_matrix<8>();
+        check_base_matrix<9>();
+        check_base_matrix<10>();
     }
 
 #if 0