From 6878d2ebe91f94e5c2bc0fe33834083f6dce85b3 Mon Sep 17 00:00:00 2001 From: Gil Date: Thu, 15 Nov 2018 09:56:45 +0000 Subject: [PATCH] Refactored linalg operations tests ( addresses #4383) (#4405) * Refactored linalg operations tests * Changed tests where `int` types are required to work with both `int` and `float` * Did not add `complex128_t`, `bool` or `char` tests * added original SE test * solver tests will have a lower threshold, as these seem to change quite a bit depending on machine and lapack/blas installations * automatic tolerance detection --- .../operations/Eigen3_operations_unittest.cc | 1608 +++++++++-------- 1 file changed, 866 insertions(+), 742 deletions(-) diff --git a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc index 4e01a70c147..3494e165c06 100644 --- a/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc +++ b/tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc @@ -11,189 +11,248 @@ using namespace shogun; using namespace linalg; using namespace Eigen; -TEST(LinalgBackendEigen, SGVector_add) +// Tolerance values for tests +template +T get_epsilon() { - const float64_t alpha = 0.3; - const float64_t beta = -1.5; + return std::numeric_limits::epsilon() * 100; +} +template <> +floatmax_t get_epsilon() +{ + return 1e-13; +} + +template +class LinalgBackendEigenAllTypesTest : public ::testing::Test +{ +}; +template +class LinalgBackendEigenNonComplexTypesTest : public ::testing::Test +{ +}; +template +class LinalgBackendEigenRealTypesTest : public ::testing::Test +{ +}; +template +class LinalgBackendEigenNonIntegerTypesTest : public ::testing::Test +{ +}; + +// TODO: make global definitions +// Definition of the 4 groups of Shogun types +// (shogun/mathematics/linalg/LinalgBackendBase.h) +// TODO: add bool and complex128_t types +typedef ::testing::Types< + int8_t, int16_t, uint16_t, int32_t, uint32_t, int64_t, uint64_t, float32_t, + float64_t, floatmax_t, char> + AllTypes; +typedef ::testing::Types< + int8_t, int16_t, uint16_t, int32_t, uint32_t, int64_t, uint64_t, float32_t, + float64_t, floatmax_t> + NonComplexTypes; +typedef ::testing::Types RealTypes; +// TODO: add complex128_t type +typedef ::testing::Types NonIntegerTypes; - SGVector A(9); - SGVector B(9); +TYPED_TEST_CASE(LinalgBackendEigenAllTypesTest, AllTypes); +TYPED_TEST_CASE(LinalgBackendEigenNonComplexTypesTest, NonComplexTypes); +TYPED_TEST_CASE(LinalgBackendEigenRealTypesTest, RealTypes); +TYPED_TEST_CASE(LinalgBackendEigenNonIntegerTypesTest, NonIntegerTypes); + +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_add) +{ + const TypeParam alpha = 1; + const TypeParam beta = 2; + + SGVector A(9); + SGVector B(9); for (index_t i = 0; i < 9; ++i) { A[i] = i; - B[i] = 0.5*i; + B[i] = 2 * i; } auto result = add(A, B, alpha, beta); for (index_t i = 0; i < 9; ++i) - EXPECT_NEAR(alpha*A[i]+beta*B[i], result[i], 1e-15); + EXPECT_NEAR( + alpha * A[i] + beta * B[i], result[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_add) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_add) { - const float64_t alpha = 0.3; - const float64_t beta = -1.5; + const TypeParam alpha = 1; + const TypeParam beta = 2; const index_t nrows = 2, ncols = 3; - SGMatrix A(nrows, ncols); - SGMatrix B(nrows, ncols); + SGMatrix A(nrows, ncols); + SGMatrix B(nrows, ncols); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) { A[i] = i; - B[i] = 0.5*i; + B[i] = 2 * i; } auto result = add(A, B, alpha, beta); - for (index_t i = 0; i < nrows*ncols; ++i) - EXPECT_NEAR(alpha*A[i]+beta*B[i], result[i], 1e-15); + for (index_t i = 0; i < nrows * ncols; ++i) + EXPECT_NEAR( + alpha * A[i] + beta * B[i], result[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_add_in_place) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_add_in_place) { - const float64_t alpha = 0.3; - const float64_t beta = -1.5; + const TypeParam alpha = 1; + const TypeParam beta = 2; - SGVector A(9), B(9), C(9); + SGVector A(9), B(9), C(9); for (index_t i = 0; i < 9; ++i) { A[i] = i; - B[i] = 0.5*i; + B[i] = 2 * i; C[i] = i; } add(A, B, A, alpha, beta); for (index_t i = 0; i < 9; ++i) - EXPECT_NEAR(alpha*C[i]+beta*B[i], A[i], 1e-15); + EXPECT_NEAR(alpha * C[i] + beta * B[i], A[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_add_in_place) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_add_in_place) { - const float64_t alpha = 0.3; - const float64_t beta = -1.5; + const TypeParam alpha = 1; + const TypeParam beta = 2; const index_t nrows = 2, ncols = 3; - SGMatrix A(nrows, ncols); - SGMatrix B(nrows, ncols); - SGMatrix C(nrows, ncols); + SGMatrix A(nrows, ncols); + SGMatrix B(nrows, ncols); + SGMatrix C(nrows, ncols); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) { A[i] = i; - B[i] = 0.5*i; + B[i] = 3 * i; C[i] = i; } add(A, B, A, alpha, beta); - for (index_t i = 0; i < nrows*ncols; ++i) - EXPECT_NEAR(alpha*C[i]+beta*B[i], A[i], 1e-15); + for (index_t i = 0; i < nrows * ncols; ++i) + EXPECT_NEAR(alpha * C[i] + beta * B[i], A[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_add_col_vec_allocated) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_add_col_vec_allocated) { - const float64_t alpha = 0.7; - const float64_t beta = -1.2; + const TypeParam alpha = 1; + const TypeParam beta = 2; const index_t nrows = 2, ncols = 3; const index_t col = 1; - SGMatrix A(nrows, ncols); - SGVector b(nrows); - SGVector result(nrows); + SGMatrix A(nrows, ncols); + SGVector b(nrows); + SGVector result(nrows); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) A[i] = i; for (index_t i = 0; i < nrows; ++i) - b[i] = 0.5*i; + b[i] = 3 * i; add_col_vec(A, col, b, result, alpha, beta); for (index_t i = 0; i < nrows; ++i) - EXPECT_NEAR(result[i], alpha*A.get_element(i, col)+beta*b[i], 1e-15); + EXPECT_NEAR( + result[i], alpha * A.get_element(i, col) + beta * b[i], + get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_add_col_vec_in_place) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_add_col_vec_in_place) { - const float64_t alpha = 0.6; - const float64_t beta = -1.3; + const TypeParam alpha = 0; + const TypeParam beta = 3; const index_t nrows = 2, ncols = 3; const index_t col = 1; - SGMatrix A(nrows, ncols); - SGVector b(nrows); + SGMatrix A(nrows, ncols); + SGVector b(nrows); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) A[i] = i; for (index_t i = 0; i < nrows; ++i) - b[i] = 0.5*i; + b[i] = 2 * i; add_col_vec(A, col, b, b, alpha, beta); for (index_t i = 0; i < nrows; ++i) - EXPECT_NEAR(b[i], alpha*A.get_element(i, col)+beta*0.5*i, 1e-15); + EXPECT_NEAR( + b[i], alpha * A.get_element(i, col) + beta * 2 * i, + get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_add_col_vec_allocated) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_add_col_vec_allocated) { - const float64_t alpha = 0.8; - const float64_t beta = -1.4; + const TypeParam alpha = 0; + const TypeParam beta = 2; const index_t nrows = 2, ncols = 3; const index_t col = 1; - SGMatrix A(nrows, ncols); - SGVector b(nrows); - SGMatrix result(nrows, ncols); + SGMatrix A(nrows, ncols); + SGVector b(nrows); + SGMatrix result(nrows, ncols); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) A[i] = i; for (index_t i = 0; i < nrows; ++i) - b[i] = 0.5*i; + b[i] = 3 * i; add_col_vec(A, col, b, result, alpha, beta); for (index_t i = 0; i < nrows; ++i) EXPECT_NEAR( result.get_element(i, col), - alpha * A.get_element(i, col) + beta * b[i], 1e-15); + alpha * A.get_element(i, col) + beta * b[i], + get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_add_col_vec_in_place) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_add_col_vec_in_place) { - const float64_t alpha = 0.9; - const float64_t beta = -1.7; + const TypeParam alpha = 1; + const TypeParam beta = 2; const index_t nrows = 2, ncols = 3; const index_t col = 1; - SGMatrix A(nrows, ncols); - SGVector b(nrows); + SGMatrix A(nrows, ncols); + SGVector b(nrows); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) A[i] = i; for (index_t i = 0; i < nrows; ++i) - b[i] = 0.5*i; + b[i] = 3 * i; add_col_vec(A, col, b, A, alpha, beta); for (index_t i = 0; i < nrows; ++i) for (index_t j = 0; j < ncols; ++j) { - float64_t a = i+j*nrows; + TypeParam a = i + j * nrows; if (j == col) - EXPECT_NEAR(A.get_element(i, j), alpha*a+beta*b[i], 1e-15); + EXPECT_NEAR( + A.get_element(i, j), alpha * a + beta * b[i], + get_epsilon()); else - EXPECT_EQ(A.get_element(i,j), a); + EXPECT_EQ(A.get_element(i, j), a); } } -TEST(LinalgBackendEigen, add_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, add_diag) { - SGMatrix A1(2, 3); - SGVector b1(2); + SGMatrix A1(2, 3); + SGVector b1(2); A1(0, 0) = 1; A1(0, 1) = 2; @@ -205,29 +264,32 @@ TEST(LinalgBackendEigen, add_diag) b1[0] = 1; b1[1] = 2; - add_diag(A1, b1, 0.5, 2.0); + const TypeParam alpha = 1.0; + const TypeParam beta = 2.0; + + add_diag(A1, b1, alpha, beta); - EXPECT_NEAR(A1(0, 0), 2.5, 1e-15); - EXPECT_NEAR(A1(0, 1), 2, 1e-15); - EXPECT_NEAR(A1(0, 2), 3, 1e-15); - EXPECT_NEAR(A1(1, 0), 4, 1e-15); - EXPECT_NEAR(A1(1, 1), 6.5, 1e-15); - EXPECT_NEAR(A1(1, 2), 6, 1e-15); + EXPECT_NEAR(A1(0, 0), 3, get_epsilon()); + EXPECT_NEAR(A1(0, 1), 2, get_epsilon()); + EXPECT_NEAR(A1(0, 2), 3, get_epsilon()); + EXPECT_NEAR(A1(1, 0), 4, get_epsilon()); + EXPECT_NEAR(A1(1, 1), 9, get_epsilon()); + EXPECT_NEAR(A1(1, 2), 6, get_epsilon()); // test error cases - SGMatrix A2(2, 2); - SGVector b2(3); - SGMatrix A3; - SGVector b3; + SGMatrix A2(2, 2); + SGVector b2(3); + SGMatrix A3; + SGVector b3; EXPECT_THROW(add_diag(A2, b2), ShogunException); EXPECT_THROW(add_diag(A2, b3), ShogunException); EXPECT_THROW(add_diag(A3, b2), ShogunException); EXPECT_THROW(add_diag(A3, b3), ShogunException); } -TEST(LinalgBackendEigen, add_ridge) +TYPED_TEST(LinalgBackendEigenAllTypesTest, add_ridge) { - SGMatrix A1(2, 3); + SGMatrix A1(2, 3); A1(0, 0) = 1; A1(0, 1) = 2; @@ -236,32 +298,34 @@ TEST(LinalgBackendEigen, add_ridge) A1(1, 1) = 5; A1(1, 2) = 6; - add_ridge(A1, 1.0); + const TypeParam alpha = 1.0; - EXPECT_NEAR(A1(0, 0), 2, 1e-15); - EXPECT_NEAR(A1(0, 1), 2, 1e-15); - EXPECT_NEAR(A1(0, 2), 3, 1e-15); - EXPECT_NEAR(A1(1, 0), 4, 1e-15); - EXPECT_NEAR(A1(1, 1), 6, 1e-15); - EXPECT_NEAR(A1(1, 2), 6, 1e-15); + add_ridge(A1, alpha); + + EXPECT_NEAR(A1(0, 0), 2, get_epsilon()); + EXPECT_NEAR(A1(0, 1), 2, get_epsilon()); + EXPECT_NEAR(A1(0, 2), 3, get_epsilon()); + EXPECT_NEAR(A1(1, 0), 4, get_epsilon()); + EXPECT_NEAR(A1(1, 1), 6, get_epsilon()); + EXPECT_NEAR(A1(1, 2), 6, get_epsilon()); // test error cases - SGMatrix A2; - EXPECT_THROW(add_ridge(A2, 1.0), ShogunException); + SGMatrix A2; + EXPECT_THROW(add_ridge(A2, alpha), ShogunException); } -TEST(LinalgBackendEigen, add_vector) +TYPED_TEST(LinalgBackendEigenAllTypesTest, add_vector) { - const float64_t alpha = 0.7; - const float64_t beta = -1.2; + const TypeParam alpha = 1; + const TypeParam beta = 2; const index_t nrows = 2, ncols = 3; - SGMatrix A(nrows, ncols); - SGMatrix result(nrows, ncols); - SGVector b(nrows); + SGMatrix A(nrows, ncols); + SGMatrix result(nrows, ncols); + SGVector b(nrows); for (index_t i = 0; i < nrows; ++i) - b[i] = 0.5 * i; + b[i] = 3 * i; for (index_t j = 0; j < ncols; ++j) for (index_t i = 0; i < nrows; ++i) { @@ -272,110 +336,112 @@ TEST(LinalgBackendEigen, add_vector) add_vector(A, b, A, alpha, beta); for (index_t i = 0; i < nrows * ncols; ++i) - EXPECT_NEAR(A[i], result[i], 1e-15); + EXPECT_NEAR(A[i], result[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_add_scalar) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_add_scalar) { const index_t n = 4; - float64_t s = -0.3; + TypeParam s = -0.3; - SGVector a(n); + SGVector a(n); for (index_t i = 0; i < (index_t)a.size(); ++i) a[i] = i; - SGVector orig = a.clone(); + SGVector orig = a.clone(); add_scalar(a, s); for (index_t i = 0; i < (index_t)a.size(); ++i) - EXPECT_NEAR(a[i], orig[i] + s, 1e-15); + EXPECT_NEAR(a[i], orig[i] + s, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_add_scalar) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_add_scalar) { const index_t r = 4, c = 3; - float64_t s = 0.4; + TypeParam s = 0.4; - SGMatrix a(r, c); + SGMatrix a(r, c); for (index_t i = 0; i < (index_t)a.size(); ++i) a[i] = i; - SGMatrix orig = a.clone(); + SGMatrix orig = a.clone(); add_scalar(a, s); for (index_t i = 0; i < (index_t)a.size(); ++i) - EXPECT_NEAR(a[i], orig[i] + s, 1e-15); + EXPECT_NEAR(a[i], orig[i] + s, get_epsilon()); } -TEST(LinalgBackendEigen, center_matrix) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_center_matrix) { const index_t n = 3; - float64_t data[] = {0.8192343, 0.13191962, 0.50888604, - 0.16857468, 0.24107738, 0.89455301, - 0.40657379, 0.07902286, 0.24319651}; - float64_t result[] = {0.25587541, -0.1173183, -0.13855711, - -0.34283925, 0.04378442, 0.29905482, - 0.08696383, 0.07353387, -0.16049771}; + TypeParam data[] = {0.5, 0.3, 0.4, 0.4, 0.5, 0.3, 0.3, 0.4, 0.5}; + TypeParam result[] = {0.1, -0.1, 0.0, 0.0, 0.1, -0.1, -0.1, 0.0, 0.1}; - SGMatrix m(data, n, n, false); + SGMatrix m(data, n, n, false); center_matrix(m); for (index_t i = 0; i < (index_t)m.size(); ++i) - EXPECT_NEAR(m[i], result[i], 1e-8); + EXPECT_NEAR(m[i], result[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_cholesky_llt_lower) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_cholesky_llt_lower) { - const index_t size=2; - SGMatrix m(size, size); + const index_t size = 2; + SGMatrix m(size, size); + // need to adapt the Eigen::Matrix to TypeParam type + typedef Matrix Mxx; - m(0,0)=2.0; - m(0,1)=1.0; - m(1,0)=1.0; - m(1,1)=2.5; + m(0, 0) = 2.0; + m(0, 1) = 1.0; + m(1, 0) = 1.0; + m(1, 1) = 2.5; - //lower triangular cholesky decomposition - SGMatrix L = cholesky_factor(m); + // lower triangular cholesky decomposition + SGMatrix L = cholesky_factor(m); - Map map_A(m.matrix,m.num_rows,m.num_cols); - Map map_L(L.matrix,L.num_rows,L.num_cols); - EXPECT_NEAR((map_A-map_L*map_L.transpose()).norm(), - 0.0, 1E-15); + Map map_A(m.matrix, m.num_rows, m.num_cols); + Map map_L(L.matrix, L.num_rows, L.num_cols); + EXPECT_NEAR((map_A - map_L * map_L.transpose()).norm(), 0.0, 1E-6); EXPECT_EQ(m.num_rows, L.num_rows); EXPECT_EQ(m.num_cols, L.num_cols); } -TEST(LinalgBackendEigen, SGMatrix_cholesky_llt_upper) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_cholesky_llt_upper) { - const index_t size=2; - SGMatrix m(size, size); + typedef Matrix Mxx; - m(0,0)=2.0; - m(0,1)=1.0; - m(1,0)=1.0; - m(1,1)=2.5; + const index_t size = 2; + SGMatrix m(size, size); + + m(0, 0) = 2.0; + m(0, 1) = 1.0; + m(1, 0) = 1.0; + m(1, 1) = 2.5; - //upper triangular cholesky decomposition - SGMatrix U = cholesky_factor(m,false); + // upper triangular cholesky decomposition + SGMatrix U = cholesky_factor(m, false); - Map map_A(m.matrix,m.num_rows,m.num_cols); - Map map_U(U.matrix,U.num_rows,U.num_cols); - EXPECT_NEAR((map_A-map_U.transpose()*map_U).norm(), - 0.0, 1E-15); + Map map_A(m.matrix, m.num_rows, m.num_cols); + Map map_U(U.matrix, U.num_rows, U.num_cols); + EXPECT_NEAR((map_A - map_U.transpose() * map_U).norm(), 0.0, 1E-6); EXPECT_EQ(m.num_rows, U.num_rows); EXPECT_EQ(m.num_cols, U.num_cols); } -TEST(LinalgBackendEigen, SGMatrix_cholesky_rank_update_upper) +TYPED_TEST(LinalgBackendEigenRealTypesTest, SGMatrix_cholesky_rank_update_upper) { + typedef Matrix Mxx; + typedef Matrix Vx; + const index_t size = 2; - SGMatrix A(size, size); - SGMatrix U(size, size); - SGVector b(size); - Map A_eig(A.matrix, size, size); - Map U_eig(U.matrix, size, size); - Map b_eig(b.vector, size); + TypeParam alpha = 1; + SGMatrix A(size, size); + SGMatrix U(size, size); + SGVector b(size); + Map A_eig(A.matrix, size, size); + Map U_eig(U.matrix, size, size); + Map b_eig(b.vector, size); U(0, 0) = 2.0; U(0, 1) = 1.0; @@ -385,28 +451,36 @@ TEST(LinalgBackendEigen, SGMatrix_cholesky_rank_update_upper) A_eig = U_eig.transpose() * U_eig; auto A2 = A.clone(); - Map A2_eig(A2.matrix, A2.num_rows, A2.num_cols); + Map A2_eig(A2.matrix, A2.num_rows, A2.num_cols); A2(0, 0) += b[0] * b[0]; A2(0, 1) += b[0] * b[1]; A2(1, 0) += b[1] * b[0]; A2(1, 1) += b[1] * b[1]; - cholesky_rank_update(U, b, 1.0, false); - EXPECT_NEAR((A2_eig - U_eig.transpose() * U_eig).norm(), 0.0, 1e-14); + cholesky_rank_update(U, b, alpha, false); + EXPECT_NEAR( + (A2_eig - U_eig.transpose() * U_eig).norm(), 0.0, + get_epsilon()); - cholesky_rank_update(U, b, -1., false); - EXPECT_NEAR((A_eig - U_eig.transpose() * U_eig).norm(), 0.0, 1e-14); + cholesky_rank_update(U, b, -alpha, false); + EXPECT_NEAR( + (A_eig - U_eig.transpose() * U_eig).norm(), 0.0, + get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_cholesky_rank_update_lower) +TYPED_TEST(LinalgBackendEigenRealTypesTest, SGMatrix_cholesky_rank_update_lower) { + typedef Matrix Mxx; + typedef Matrix Vx; + const index_t size = 2; - SGMatrix A(size, size); - SGMatrix L(size, size); - SGVector b(size); - Map A_eig(A.matrix, size, size); - Map L_eig(L.matrix, size, size); - Map b_eig(b.vector, size); + TypeParam alpha = 1; + SGMatrix A(size, size); + SGMatrix L(size, size); + SGVector b(size); + Map A_eig(A.matrix, size, size); + Map L_eig(L.matrix, size, size); + Map b_eig(b.vector, size); L(0, 0) = 2.0; L(1, 0) = 1.0; @@ -416,23 +490,27 @@ TEST(LinalgBackendEigen, SGMatrix_cholesky_rank_update_lower) A_eig = L_eig * L_eig.transpose(); auto A2 = A.clone(); - Map A2_eig(A2.matrix, A2.num_rows, A2.num_cols); + Map A2_eig(A2.matrix, A2.num_rows, A2.num_cols); A2(0, 0) += b[0] * b[0]; A2(0, 1) += b[0] * b[1]; A2(1, 0) += b[1] * b[0]; A2(1, 1) += b[1] * b[1]; - cholesky_rank_update(L, b); - EXPECT_NEAR((A2_eig - L_eig * L_eig.transpose()).norm(), 0.0, 1e-14); + cholesky_rank_update(L, b, alpha); + EXPECT_NEAR( + (A2_eig - L_eig * L_eig.transpose()).norm(), 0.0, + get_epsilon()); - cholesky_rank_update(L, b, -1.); - EXPECT_NEAR((A_eig - L_eig * L_eig.transpose()).norm(), 0.0, 1e-14); + cholesky_rank_update(L, b, -alpha); + EXPECT_NEAR( + (A_eig - L_eig * L_eig.transpose()).norm(), 0.0, + get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_cholesky_ldlt_lower) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_cholesky_ldlt_lower) { const index_t size = 3; - SGMatrix m(size, size); + SGMatrix m(size, size); m(0, 0) = 0.0; m(0, 1) = 0.0; m(0, 2) = 0.0; @@ -443,60 +521,60 @@ TEST(LinalgBackendEigen, SGMatrix_cholesky_ldlt_lower) m(2, 1) = 2.0; m(2, 2) = 3.0; - SGMatrix L(size, size); - SGVector d(size); + SGMatrix L(size, size); + SGVector d(size); SGVector p(size); linalg::ldlt_factor(m, L, d, p); - EXPECT_NEAR(d[0], 3.0, 1e-15); - EXPECT_NEAR(d[1], -0.333333333333333, 1e-15); - EXPECT_NEAR(d[2], 0.0, 1e-15); + EXPECT_NEAR(d[0], 3.0, get_epsilon()); + EXPECT_NEAR(d[1], -0.333333333333333, get_epsilon()); + EXPECT_NEAR(d[2], 0.0, get_epsilon()); - EXPECT_NEAR(L(0, 0), 1.0, 1e-15); - EXPECT_NEAR(L(0, 1), 0.0, 1e-15); - EXPECT_NEAR(L(0, 2), 0.0, 1e-15); - EXPECT_NEAR(L(1, 0), 0.666666666666666, 1e-15); - EXPECT_NEAR(L(1, 1), 1.0, 1e-15); - EXPECT_NEAR(L(1, 2), 0.0, 1e-15); - EXPECT_NEAR(L(2, 0), 0.0, 1e-15); - EXPECT_NEAR(L(2, 1), 0.0, 1e-15); - EXPECT_NEAR(L(2, 2), 1.0, 1e-15); + EXPECT_NEAR(L(0, 0), 1.0, get_epsilon()); + EXPECT_NEAR(L(0, 1), 0.0, get_epsilon()); + EXPECT_NEAR(L(0, 2), 0.0, get_epsilon()); + EXPECT_NEAR(L(1, 0), 0.666666666666666, get_epsilon()); + EXPECT_NEAR(L(1, 1), 1.0, get_epsilon()); + EXPECT_NEAR(L(1, 2), 0.0, get_epsilon()); + EXPECT_NEAR(L(2, 0), 0.0, get_epsilon()); + EXPECT_NEAR(L(2, 1), 0.0, get_epsilon()); + EXPECT_NEAR(L(2, 2), 1.0, get_epsilon()); EXPECT_EQ(p[0], 2); EXPECT_EQ(p[1], 1); EXPECT_EQ(p[2], 2); } -TEST(LinalgBackendEigen, SGMatrix_cholesky_solver) +TYPED_TEST(LinalgBackendEigenRealTypesTest, SGMatrix_cholesky_solver) { - const index_t size=2; - SGMatrix A(size, size); - A(0,0)=2.0; - A(0,1)=1.0; - A(1,0)=1.0; - A(1,1)=2.5; + const index_t size = 2; + SGMatrix A(size, size); + A(0, 0) = 2.0; + A(0, 1) = 1.0; + A(1, 0) = 1.0; + A(1, 1) = 2.5; - SGVector b(size); + SGVector b(size); b[0] = 10; b[1] = 13; - SGVector x_ref(size); + SGVector x_ref(size); x_ref[0] = 3; x_ref[1] = 4; - SGMatrix L = cholesky_factor(A); - SGVector x_cal = cholesky_solver(L, b); + SGMatrix L = cholesky_factor(A); + SGVector x_cal = cholesky_solver(L, b); - EXPECT_NEAR(x_ref[0], x_cal[0], 1E-15); - EXPECT_NEAR(x_ref[1], x_cal[1], 1E-15); + EXPECT_NEAR(x_ref[0], x_cal[0], get_epsilon()); + EXPECT_NEAR(x_ref[1], x_cal[1], get_epsilon()); EXPECT_EQ(x_ref.size(), x_cal.size()); } -TEST(LinalgBackendEigen, SGMatrix_ldlt_solver) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_ldlt_solver) { const index_t size = 3; - SGMatrix A(size, size); + SGMatrix A(size, size); A(0, 0) = 0.0; A(0, 1) = 0.0; A(0, 2) = 0.0; @@ -507,86 +585,88 @@ TEST(LinalgBackendEigen, SGMatrix_ldlt_solver) A(2, 1) = 2.0; A(2, 2) = 3.0; - SGVector b(size); + SGVector b(size); b[0] = 0.0; b[1] = 5.0; b[2] = 11.0; - SGVector x_ref(size), x(size); + SGVector x_ref(size), x(size); x_ref[0] = 0.0; x_ref[1] = 7.0; x_ref[2] = -1.0; - SGMatrix L(size, size); - SGVector d(size); + SGMatrix L(size, size); + SGVector d(size); SGVector p(size); linalg::ldlt_factor(A, L, d, p, true); x = linalg::ldlt_solver(L, d, p, b, true); for (auto i : range(size)) - EXPECT_NEAR(x[i], x_ref[i], 1e-15); + EXPECT_NEAR(x[i], x_ref[i], get_epsilon()); linalg::ldlt_factor(A, L, d, p, false); x = linalg::ldlt_solver(L, d, p, b, false); for (auto i : range(size)) - EXPECT_NEAR(x[i], x_ref[i], 1e-15); + EXPECT_NEAR(x[i], x_ref[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_cross_entropy) + +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_cross_entropy) { - SGMatrix A(4, 3); - SGMatrix B(4, 3); + SGMatrix A(4, 3); + SGMatrix B(4, 3); - int32_t size = A.num_rows * A.num_cols; - for (float64_t i = 0; i < size; ++i) + uint32_t size = A.num_rows * A.num_cols; + for (TypeParam i = 0; i < size; ++i) { A[i] = i / size; B[i] = (i / size) * 0.5; } float64_t ref = 0; - for (int32_t i = 0; i < size; i++) - ref += A[i] * std::log(B[i] + 1e-30); + for (uint32_t i = 0; i < size; i++) + ref += A[i] * std::log(B[i] + 1e-15); ref *= -1; auto result = linalg::cross_entropy(A, B); - EXPECT_NEAR(ref, result, 1e-15); + EXPECT_NEAR(ref, result, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_pinv_psd) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_pinv_psd) { - float64_t A_data[] = {2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0}; + TypeParam A_data[] = {2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0}; // inverse generated by scipy pinv - float64_t scipy_result_data[] = {0.75, 0.5, 0.25, 0.5, 1.0, + TypeParam scipy_result_data[] = {0.75, 0.5, 0.25, 0.5, 1.0, 0.5, 0.25, 0.5, 0.75}; - SGMatrix A(A_data, 3, 3, false); - SGMatrix result(scipy_result_data, 3, 3, false); + SGMatrix A(A_data, 3, 3, false); + SGMatrix result(scipy_result_data, 3, 3, false); - SGMatrix identity_matrix(3, 3); + SGMatrix identity_matrix(3, 3); linalg::identity(identity_matrix); // using symmetric eigen solver - SGMatrix A_pinvh(3, 3); + SGMatrix A_pinvh(3, 3); linalg::pinvh(A, A_pinvh); // using singular value decomposition - SGMatrix A_pinv(3, 3); + SGMatrix A_pinv(3, 3); linalg::pinv(A, A_pinv); - SGMatrix I_check = linalg::matrix_prod(A, A_pinvh); + SGMatrix I_check = linalg::matrix_prod(A, A_pinvh); for (auto i : range(3)) { for (auto j : range(3)) { - EXPECT_NEAR(identity_matrix(i, j), I_check(i, j), 1e-14); - EXPECT_NEAR(result(i, j), A_pinvh(i, j), 1e-14); - EXPECT_NEAR(result(i, j), A_pinv(i, j), 1e-14); + EXPECT_NEAR( + identity_matrix(i, j), I_check(i, j), get_epsilon()); + EXPECT_NEAR(result(i, j), A_pinvh(i, j), get_epsilon()); + EXPECT_NEAR(result(i, j), A_pinv(i, j), get_epsilon()); } } - // no memory errors + // no memory errors EXPECT_NO_THROW(linalg::pinvh(A, A)); } -TEST(LinalgBackendEigen, SGMatrix_pinv_2x4) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_pinv_2x4) { - SGMatrix A(2, 4); + SGMatrix A(2, 4); A(0, 0) = 1; A(0, 1) = 1; A(0, 2) = 1; @@ -596,39 +676,40 @@ TEST(LinalgBackendEigen, SGMatrix_pinv_2x4) A(1, 2) = 7; A(1, 3) = 9; - SGMatrix identity_matrix(2, 2); + SGMatrix identity_matrix(2, 2); linalg::identity(identity_matrix); - SGMatrix A_pinverse(4, 2); + SGMatrix A_pinverse(4, 2); linalg::pinv(A, A_pinverse); - SGMatrix I_check = linalg::matrix_prod(A, A_pinverse); + SGMatrix I_check = linalg::matrix_prod(A, A_pinverse); for (auto i : range(2)) { for (auto j : range(2)) { - EXPECT_NEAR(identity_matrix(i, j), I_check(i, j), 1e-14); + EXPECT_NEAR( + identity_matrix(i, j), I_check(i, j), get_epsilon()); } } // compare result with scipy pinv - EXPECT_NEAR(A_pinverse(0, 0), 2.0, 1e-14); - EXPECT_NEAR(A_pinverse(0, 1), -0.25, 1e-14); + EXPECT_NEAR(A_pinverse(0, 0), 2.0, get_epsilon()); + EXPECT_NEAR(A_pinverse(0, 1), -0.25, get_epsilon()); - EXPECT_NEAR(A_pinverse(1, 0), 0.25, 1e-14); - EXPECT_NEAR(A_pinverse(1, 1), 0.0, 1e-14); + EXPECT_NEAR(A_pinverse(1, 0), 0.25, get_epsilon()); + EXPECT_NEAR(A_pinverse(1, 1), 0.0, get_epsilon()); - EXPECT_NEAR(A_pinverse(2, 0), 0.25, 1e-14); - EXPECT_NEAR(A_pinverse(2, 1), 0.0, 1e-14); + EXPECT_NEAR(A_pinverse(2, 0), 0.25, get_epsilon()); + EXPECT_NEAR(A_pinverse(2, 1), 0.0, get_epsilon()); - EXPECT_NEAR(A_pinverse(3, 0), -1.5, 1e-14); - EXPECT_NEAR(A_pinverse(3, 1), 0.25, 1e-14); + EXPECT_NEAR(A_pinverse(3, 0), -1.5, get_epsilon()); + EXPECT_NEAR(A_pinverse(3, 1), 0.25, get_epsilon()); // incorrect dimension EXPECT_THROW(linalg::pinv(A, A), ShogunException); } -TEST(LinalgBackendEigen, SGMatrix_pinv_4x2) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_pinv_4x2) { - SGMatrix A(4, 2); + SGMatrix A(4, 2); A(0, 0) = 2.0; A(0, 1) = -0.25; A(1, 0) = 0.25; @@ -638,60 +719,61 @@ TEST(LinalgBackendEigen, SGMatrix_pinv_4x2) A(3, 0) = -1.5; A(3, 1) = 0.25; - SGMatrix identity_matrix(2, 2); + SGMatrix identity_matrix(2, 2); linalg::identity(identity_matrix); - SGMatrix A_pinverse(2, 4); + SGMatrix A_pinverse(2, 4); linalg::pinv(A, A_pinverse); - SGMatrix I_check = linalg::matrix_prod(A_pinverse, A); + SGMatrix I_check = linalg::matrix_prod(A_pinverse, A); for (auto i : range(2)) { for (auto j : range(2)) { - EXPECT_NEAR(identity_matrix(i, j), I_check(i, j), 1e-14); + EXPECT_NEAR( + identity_matrix(i, j), I_check(i, j), get_epsilon()); } } // compare with results from scipy - EXPECT_NEAR(A_pinverse(0, 0), 1.0, 1e-14); - EXPECT_NEAR(A_pinverse(0, 1), 1.0, 1e-14); - EXPECT_NEAR(A_pinverse(0, 2), 1.0, 1e-14); - EXPECT_NEAR(A_pinverse(0, 3), 1.0, 1e-14); + EXPECT_NEAR(A_pinverse(0, 0), 1.0, get_epsilon()); + EXPECT_NEAR(A_pinverse(0, 1), 1.0, get_epsilon()); + EXPECT_NEAR(A_pinverse(0, 2), 1.0, get_epsilon()); + EXPECT_NEAR(A_pinverse(0, 3), 1.0, get_epsilon()); - EXPECT_NEAR(A_pinverse(1, 0), 5.0, 1e-14); - EXPECT_NEAR(A_pinverse(1, 1), 7.0, 1e-14); - EXPECT_NEAR(A_pinverse(1, 2), 7.0, 1e-14); - EXPECT_NEAR(A_pinverse(1, 3), 9.0, 1e-14); + EXPECT_NEAR(A_pinverse(1, 0), 5.0, get_epsilon()); + EXPECT_NEAR(A_pinverse(1, 1), 7.0, get_epsilon()); + EXPECT_NEAR(A_pinverse(1, 2), 7.0, get_epsilon()); + EXPECT_NEAR(A_pinverse(1, 3), 9.0, get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_dot) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_dot) { const index_t size = 3; - SGVector a(size), b(size); + SGVector a(size), b(size); a.range_fill(0); b.range_fill(0); auto result = dot(a, b); - EXPECT_NEAR(result, 5, 1E-15); + EXPECT_NEAR(result, 5, get_epsilon()); } -TEST(LinalgBackendEigen, eigensolver) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, eigensolver) { const index_t n = 4; - float64_t data[] = {0.09987322, 0.80575314, 0.79068641, 0.69989667, + TypeParam data[] = {0.09987322, 0.80575314, 0.79068641, 0.69989667, 0.62323516, 0.16837367, 0.85027625, 0.60165948, 0.04898732, 0.96701123, 0.51683275, 0.51116495, 0.18277926, 0.6179262, 0.43745891, 0.63685464}; - float64_t result_eigenvectors[] = { - -0.63494074, 0.75831593, -0.14014031, 0.04656076, - 0.82257205, -0.28671857, -0.44196422, -0.21409185, - -0.005932, -0.20233724, -0.52285555, 0.82803776, - -0.23930111, -0.56199714, -0.57298901, -0.54642272}; - float64_t result_eigenvalues[] = {-0.6470538, -0.19125664, 0.16205101, + TypeParam result_eigenvectors[] = { + -0.63494074, 0.75831593, -0.1401403109, 0.04656076, + 0.82257205, -0.286718557, -0.44196422, -0.214091861, + -0.005932, -0.20233723, -0.52285555, 0.82803776, + -0.23930111, -0.56199714, -0.57298901, -0.54642272}; + TypeParam result_eigenvalues[] = {-0.6470538, -0.19125664, 0.16205101, 2.0981937}; - SGMatrix m(data, n, n, false); - SGMatrix eigenvectors(n, n); - SGVector eigenvalues(n); + SGMatrix m(data, n, n, false); + SGMatrix eigenvectors(n, n); + SGVector eigenvalues(n); eigen_solver(m, eigenvalues, eigenvectors); @@ -699,35 +781,35 @@ TEST(LinalgBackendEigen, eigensolver) for (index_t i = 0; i < n; ++i) { index_t idx = args[i]; - EXPECT_NEAR(eigenvalues[idx], result_eigenvalues[i], 1e-7); + EXPECT_NEAR(eigenvalues[idx], result_eigenvalues[i], 1e-6); auto s = CMath::sign(eigenvectors[idx * n] * result_eigenvectors[i * n]); for (index_t j = 0; j < n; ++j) EXPECT_NEAR( eigenvectors[idx * n + j], s * result_eigenvectors[i * n + j], - 1e-7); + 1e-6); } } -TEST(LinalgBackendEigen, eigensolver_symmetric) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, eigensolver_symmetric) { const index_t n = 4; - float64_t data[] = {0.09987322, 0.80575314, 0.04898732, 0.69989667, + TypeParam data[] = {0.09987322, 0.80575314, 0.04898732, 0.69989667, 0.80575314, 0.16837367, 0.96701123, 0.6179262, 0.04898732, 0.96701123, 0.51683275, 0.43745891, 0.69989667, 0.6179262, 0.43745891, 0.63685464}; - float64_t result_eigenvectors[] = { + TypeParam result_eigenvectors[] = { -0.54618542, 0.69935447, -0.45219663, 0.09001671, -0.56171388, -0.41397154, 0.17642953, 0.69424612, -0.46818396, 0.16780603, 0.73247599, -0.46489119, 0.40861077, 0.55800718, 0.47735703, 0.542029037}; - float64_t result_eigenvalues[] = {-1.00663298, -0.18672196, 0.42940933, + TypeParam result_eigenvalues[] = {-1.00663298, -0.18672196, 0.42940933, 2.18587989}; - SGMatrix m(data, n, n, false); - SGMatrix eigenvectors(n, n); - SGVector eigenvalues(n); + SGMatrix m(data, n, n, false); + SGMatrix eigenvectors(n, n); + SGVector eigenvalues(n); eigen_solver(m, eigenvalues, eigenvectors); @@ -735,95 +817,104 @@ TEST(LinalgBackendEigen, eigensolver_symmetric) for (index_t i = 0; i < n; ++i) { index_t idx = args[i]; - EXPECT_NEAR(eigenvalues[idx], result_eigenvalues[i], 1e-7); + EXPECT_NEAR(eigenvalues[idx], result_eigenvalues[i], 1e-5); auto s = CMath::sign(eigenvectors[idx * n] * result_eigenvectors[i * n]); for (index_t j = 0; j < n; ++j) EXPECT_NEAR( eigenvectors[idx * n + j], s * result_eigenvectors[i * n + j], - 1e-7); + 1e-5); } } -TEST(LinalgBackendEigen, SGMatrix_elementwise_product) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_elementwise_product) { - const auto m = 3; - SGMatrix A(m, m); - SGMatrix B(m, m); + const auto m = 2; + SGMatrix A(m, m); + SGMatrix B(m, m); for (auto i : range(m * m)) { A[i] = i; - B[i] = 0.5*i; + B[i] = 2 * i; } auto result = element_prod(A, B); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(i, j) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(i, j) * B(i, j), get_epsilon()); result = element_prod(A, B, true, false); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(i, j), get_epsilon()); result = element_prod(A, B, false, true); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(i, j), get_epsilon()); result = element_prod(A, B, true, true); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(j, i), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(j, i), get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_elementwise_product_in_place) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, SGMatrix_elementwise_product_in_place) { - const auto m = 3; - SGMatrix A(m, m); - SGMatrix B(m, m); - SGMatrix result(m, m); + const auto m = 2; + SGMatrix A(m, m); + SGMatrix B(m, m); + SGMatrix result(m, m); for (auto i : range(m * m)) { A[i] = i; - B[i] = 0.5*i; + B[i] = 2 * i; } element_prod(A, B, result); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(i, j) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(i, j) * B(i, j), get_epsilon()); element_prod(A, B, result, true, false); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(i, j), get_epsilon()); element_prod(A, B, result, false, true); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(i, j), get_epsilon()); element_prod(A, B, result, true, true); for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(j, i), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(j, i), get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_block_elementwise_product) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_block_elementwise_product) { const index_t nrows = 2; const index_t ncols = 3; - SGMatrix A(nrows,ncols); - SGMatrix B(ncols,nrows); + SGMatrix A(nrows, ncols); + SGMatrix B(ncols, nrows); for (auto i : range(nrows)) for (auto j : range(ncols)) @@ -842,7 +933,8 @@ TEST(LinalgBackendEigen, SGMatrix_block_elementwise_product) for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(i, j) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(i, j) * B(i, j), get_epsilon()); result = element_prod(A_block, B_block, true, false); @@ -851,7 +943,8 @@ TEST(LinalgBackendEigen, SGMatrix_block_elementwise_product) for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(i, j), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(i, j), get_epsilon()); result = element_prod(A_block, B_block, false, true); @@ -860,7 +953,8 @@ TEST(LinalgBackendEigen, SGMatrix_block_elementwise_product) for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(i, j) * B(j, i), 1E-15); + EXPECT_NEAR( + result(i, j), A(i, j) * B(j, i), get_epsilon()); result = element_prod(A_block, B_block, true, true); @@ -869,430 +963,446 @@ TEST(LinalgBackendEigen, SGMatrix_block_elementwise_product) for (auto i : range(m)) for (auto j : range(m)) - EXPECT_NEAR(result(i, j), A(j, i) * B(j, i), 1E-15); + EXPECT_NEAR( + result(i, j), A(j, i) * B(j, i), get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_elementwise_product) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_elementwise_product) { const index_t len = 4; - SGVector a(len); - SGVector b(len); - SGVector c(len); + SGVector a(len); + SGVector b(len); + SGVector c(len); for (index_t i = 0; i < len; ++i) { a[i] = i; - b[i] = 0.5 * i; + b[i] = 2 * i; } c = element_prod(a, b); for (index_t i = 0; i < len; ++i) - EXPECT_NEAR(a[i] * b[i], c[i], 1e-15); + EXPECT_NEAR(a[i] * b[i], c[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_elementwise_product_in_place) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, SGVector_elementwise_product_in_place) { const index_t len = 4; - SGVector a(len); - SGVector b(len); - SGVector c(len); + SGVector a(len); + SGVector b(len); + SGVector c(len); for (index_t i = 0; i < len; ++i) { a[i] = i; - b[i] = 0.5 * i; + b[i] = 2 * i; c[i] = i; } element_prod(a, b, a); for (index_t i = 0; i < len; ++i) - EXPECT_NEAR(c[i] * b[i], a[i], 1e-15); + EXPECT_NEAR(c[i] * b[i], a[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_exponent) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGVector_exponent) { const index_t len = 4; - SGVector a(len); - a[0] = -2.4; - a[1] = 0; - a[2] = 0.5; - a[3] = 3.9; + SGVector a(len); + a[0] = 0; + a[1] = 1; + a[2] = 2; + a[3] = 3; auto result = exponent(a); - EXPECT_NEAR(result[0], 0.090717953289413, 1E-15); - EXPECT_NEAR(result[1], 1.0, 1E-15); - EXPECT_NEAR(result[2], 1.648721270700128, 1E-15); - EXPECT_NEAR(result[3], 49.40244910553017, 1E-15); + EXPECT_NEAR(result[0], 1.0, get_epsilon()); + EXPECT_NEAR(result[1], 2.718281828459045, get_epsilon()); + EXPECT_NEAR(result[2], 7.3890560989306495, get_epsilon()); + EXPECT_NEAR(result[3], 20.085536923187664, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_exponent) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_exponent) { const index_t n = 2; - SGMatrix a(n, n); - a[0] = -2.4; - a[1] = 0; - a[2] = 0.5; - a[3] = 3.9; + SGMatrix a(n, n); + a[0] = 0; + a[1] = 1; + a[2] = 2; + a[3] = 3; auto result = exponent(a); - EXPECT_NEAR(result[0], 0.090717953289413, 1E-15); - EXPECT_NEAR(result[1], 1.0, 1E-15); - EXPECT_NEAR(result[2], 1.648721270700128, 1E-15); - EXPECT_NEAR(result[3], 49.40244910553017, 1E-15); + EXPECT_NEAR(result[0], 1.0, get_epsilon()); + EXPECT_NEAR(result[1], 2.718281828459045, get_epsilon()); + EXPECT_NEAR(result[2], 7.3890560989306495, get_epsilon()); + EXPECT_NEAR(result[3], 20.085536923187664, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_identity) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_identity) { const index_t n = 4; - SGMatrix A(n, n); + SGMatrix A(n, n); identity(A); for (index_t i = 0; i < n; ++i) for (index_t j = 0; j < n; ++j) - EXPECT_EQ(A.get_element(i, j), (i==j)); + EXPECT_EQ(A.get_element(i, j), (i == j)); } -TEST(LinalgBackendEigen, logistic) +// TODO: write test for int types +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, logistic) { - SGMatrix A(3,3); - SGMatrix B(3,3); + SGMatrix A(3, 3); + SGMatrix B(3, 3); - range_fill(A, 0.0); + for (index_t i = 0; i < 9; ++i) + A[i] = i; B.zero(); linalg::logistic(A, B); for (index_t i = 0; i < 9; ++i) - EXPECT_NEAR(1.0 / (1 + std::exp(-1 * A[i])), B[i], 1e-15); + EXPECT_NEAR( + 1.0 / (1 + std::exp(-1 * A[i])), B[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_SGVector_matrix_prod) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_SGVector_matrix_prod) { - const index_t rows=4; - const index_t cols=3; + const index_t rows = 4; + const index_t cols = 3; - SGMatrix A(rows, cols); - SGVector b(cols); + SGMatrix A(rows, cols); + SGVector b(cols); for (index_t i = 0; i < cols; ++i) { for (index_t j = 0; j < rows; ++j) A(j, i) = i * rows + j; - b[i]=0.5 * i; + b[i] = 2 * i; } auto x = matrix_prod(A, b); - float64_t ref[] = {10, 11.5, 13, 14.5}; + TypeParam ref[] = {40, 46, 52, 58}; EXPECT_EQ(x.vlen, A.num_rows); - for (index_t i = 0; i < cols; ++i) - EXPECT_NEAR(x[i], ref[i], 1e-15); + for (index_t i = 0; i < rows; ++i) + EXPECT_NEAR(x[i], ref[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_SGVector_matrix_prod_transpose) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_matrix_prod_transpose) { - const index_t rows=4; - const index_t cols=3; + const index_t rows = 4; + const index_t cols = 3; - SGMatrix A(cols, rows); - SGVector b(cols); + SGMatrix A(cols, rows); + SGVector b(cols); for (index_t i = 0; i < cols; ++i) { for (index_t j = 0; j < rows; ++j) A(i, j) = i * cols + j; - b[i] = 0.5 * i; + b[i] = 2 * i; } auto x = matrix_prod(A, b, true); - float64_t ref[] = {7.5, 9, 10.5, 14.5}; + TypeParam ref[] = {30, 36, 42}; EXPECT_EQ(x.vlen, A.num_cols); for (index_t i = 0; i < cols; ++i) - EXPECT_NEAR(x[i], ref[i], 1e-15); + EXPECT_NEAR(x[i], ref[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_SGVector_matrix_prod_in_place) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, SGMatrix_SGVector_matrix_prod_in_place) { - const index_t rows=4; - const index_t cols=3; + const index_t rows = 4; + const index_t cols = 3; - SGMatrix A(rows, cols); - SGVector b(cols); - SGVector x(rows); + SGMatrix A(rows, cols); + SGVector b(cols); + SGVector x(rows); - for (index_t i = 0; i()); } -TEST(LinalgBackendEigen, SGMatrix_SGVector_matrix_prod_in_place_transpose) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, + SGMatrix_SGVector_matrix_prod_in_place_transpose) { - const index_t rows=4; - const index_t cols=3; + const index_t rows = 4; + const index_t cols = 3; - SGMatrix A(cols, rows); - SGVector b(cols); - SGVector x(rows); + SGMatrix A(cols, rows); + SGVector b(cols); + SGVector x(rows); for (index_t i = 0; i < cols; ++i) { for (index_t j = 0; j < rows; ++j) A(i, j) = i * cols + j; - b[i] = 0.5 * i; + b[i] = 2 * i; } matrix_prod(A, b, x, true); - float64_t ref[] = {7.5, 9, 10.5, 14.5}; + TypeParam ref[] = {30, 36, 42}; for (index_t i = 0; i < cols; ++i) - EXPECT_NEAR(x[i], ref[i], 1e-15); + EXPECT_NEAR(x[i], ref[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_matrix_product) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim1, dim2); - SGMatrix B(dim2, dim3); + const index_t dim1 = 2, dim2 = 4, dim3 = 2; + SGMatrix A(dim1, dim2); + SGMatrix B(dim2, dim3); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5*i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; auto cal = linalg::matrix_prod(A, B); - float64_t ref[] = {14, 17, 38, 49, 62, 81}; + TypeParam ref[] = {28, 34, 76, 98}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product_transpose_A) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_matrix_product_transpose_A) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim2, dim1); - SGMatrix B(dim2, dim3); + const index_t dim1 = 2, dim2 = 3, dim3 = 3; + SGMatrix A(dim2, dim1); + SGMatrix B(dim2, dim3); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5*i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; auto cal = linalg::matrix_prod(A, B, true); - float64_t ref[] = {7, 19, 19, 63, 31, 107}; + TypeParam ref[] = {5, 14, 14, 50, 23, 86}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product_transpose_B) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_matrix_product_transpose_B) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim1, dim2); - SGMatrix B(dim3, dim2); + const index_t dim1 = 2, dim2 = 3, dim3 = 3; + SGMatrix A(dim1, dim2); + SGMatrix B(dim3, dim2); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5*i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; auto cal = linalg::matrix_prod(A, B, false, true); - float64_t ref[] = {42, 51, 48, 59, 54, 67}; + TypeParam ref[] = {30, 39, 36, 48, 42, 57}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product_transpose_A_B) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, SGMatrix_matrix_product_transpose_A_B) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim2, dim1); - SGMatrix B(dim3, dim2); + const index_t dim1 = 2, dim2 = 3, dim3 = 3; + SGMatrix A(dim2, dim1); + SGMatrix B(dim3, dim2); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5*i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; auto cal = linalg::matrix_prod(A, B, true, true); - float64_t ref[] = {21, 57, 24, 68, 27, 79}; + TypeParam ref[] = {15, 42, 18, 54, 21, 66}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product_in_place) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_matrix_product_in_place) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim1, dim2); - SGMatrix B(dim2, dim3); - SGMatrix cal(dim1, dim3); + const index_t dim1 = 2, dim2 = 3, dim3 = 3; + SGMatrix A(dim1, dim2); + SGMatrix B(dim2, dim3); + SGMatrix cal(dim1, dim3); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5 * i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; cal.zero(); linalg::matrix_prod(A, B, cal); - float64_t ref[] = {14, 17, 38, 49, 62, 81}; + TypeParam ref[] = {10, 13, 28, 40, 46, 67}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product_in_place_transpose_A) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, + SGMatrix_matrix_product_in_place_transpose_A) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim2, dim1); - SGMatrix B(dim2, dim3); - SGMatrix cal(dim1, dim3); + const index_t dim1 = 2, dim2 = 3, dim3 = 3; + SGMatrix A(dim2, dim1); + SGMatrix B(dim2, dim3); + SGMatrix cal(dim1, dim3); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5*i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; cal.zero(); linalg::matrix_prod(A, B, cal, true); - float64_t ref[] = {7, 19, 19, 63, 31, 107}; + TypeParam ref[] = {5, 14, 14, 50, 23, 86}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product_in_place_transpose_B) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, + SGMatrix_matrix_product_in_place_transpose_B) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim1, dim2); - SGMatrix B(dim3, dim2); - SGMatrix cal(dim1, dim3); + const index_t dim1 = 2, dim2 = 3, dim3 = 3; + SGMatrix A(dim1, dim2); + SGMatrix B(dim3, dim2); + SGMatrix cal(dim1, dim3); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5*i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; cal.zero(); linalg::matrix_prod(A, B, cal, false, true); - float64_t ref[] = {42, 51, 48, 59, 54, 67}; + TypeParam ref[] = {30, 39, 36, 48, 42, 57}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGMatrix_matrix_product_in_place_transpose_A_B) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, + SGMatrix_matrix_product_in_place_transpose_A_B) { - const index_t dim1 = 2, dim2 = 4, dim3 = 3; - SGMatrix A(dim2, dim1); - SGMatrix B(dim3, dim2); - SGMatrix cal(dim1, dim3); + const index_t dim1 = 2, dim2 = 3, dim3 = 3; + SGMatrix A(dim2, dim1); + SGMatrix B(dim3, dim2); + SGMatrix cal(dim1, dim3); - for (index_t i = 0; i < dim1*dim2; ++i) + for (index_t i = 0; i < dim1 * dim2; ++i) A[i] = i; - for (index_t i = 0; i < dim2*dim3; ++i) - B[i] = 0.5*i; + for (index_t i = 0; i < dim2 * dim3; ++i) + B[i] = i; cal.zero(); linalg::matrix_prod(A, B, cal, true, true); - float64_t ref[] = {21, 57, 24, 68, 27, 79}; + TypeParam ref[] = {15, 42, 18, 54, 21, 66}; EXPECT_EQ(dim1, cal.num_rows); EXPECT_EQ(dim3, cal.num_cols); - for (index_t i = 0; i < dim1*dim3; ++i) + for (index_t i = 0; i < dim1 * dim3; ++i) EXPECT_EQ(ref[i], cal[i]); } -TEST(LinalgBackendEigen, SGVector_max) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_max) { - SGVector A(9); + SGVector A(9); - float64_t a[] = {1, 2, 5, 8, 3, 1, 0, -1, 4}; + TypeParam a[] = {1, 2, 5, 8, 3, 1, 0, 2, 4}; for (index_t i = 0; i < A.size(); ++i) A[i] = a[i]; - EXPECT_NEAR(8, max(A), 1e-15); + EXPECT_NEAR(8, max(A), get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_max) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_max) { - const index_t nrows = 2, ncols = 3; - SGMatrix A(nrows, ncols); + const index_t nrows = 3, ncols = 3; + SGMatrix A(nrows, ncols); - float64_t a[] = {1, 2, 5, 8, 3, 1, 0, -1, 12}; + TypeParam a[] = {1, 2, 5, 8, 3, 1, 0, 2, 12}; - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) A[i] = a[i]; - EXPECT_NEAR(8, max(A), 1e-15); + EXPECT_NEAR(12, max(A), get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_mean) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_mean) { - const index_t size = 6; - SGVector vec(size); + const index_t size = 9; + SGVector vec(size); vec.range_fill(0); auto result = mean(vec); - EXPECT_NEAR(result, 2.5, 1E-15); + EXPECT_NEAR(result, 4, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_mean) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_mean) { - const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + const index_t nrows = 3, ncols = 3; + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows * ncols; ++i) mat[i] = i; auto result = mean(mat); - EXPECT_NEAR(result, 2.5, 1E-15); + EXPECT_NEAR(result, 4, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_multiply_by_logistic_derivative) +TYPED_TEST( + LinalgBackendEigenAllTypesTest, SGMatrix_multiply_by_logistic_derivative) { - SGMatrix A(3, 3); - SGMatrix B(3, 3); + SGMatrix A(3, 3); + SGMatrix B(3, 3); - for (float64_t i = 0; i < 9; ++i) + for (TypeParam i = 9; i < 9; i += 9) { A[i] = i / 9; B[i] = i; @@ -1301,15 +1411,17 @@ TEST(LinalgBackendEigen, SGMatrix_multiply_by_logistic_derivative) linalg::multiply_by_logistic_derivative(A, B); for (index_t i = 0; i < 9; ++i) - EXPECT_NEAR(i * A[i] * (1.0 - A[i]), B[i], 1e-15); + EXPECT_NEAR(i * A[i] * (1.0 - A[i]), B[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_multiply_by_rectified_linear_derivative) +TYPED_TEST( + LinalgBackendEigenNonIntegerTypesTest, + SGMatrix_multiply_by_rectified_linear_derivative) { - SGMatrix A(3, 3); - SGMatrix B(3, 3); + SGMatrix A(3, 3); + SGMatrix B(3, 3); - for (float64_t i = 0; i < 9; ++i) + for (TypeParam i = 0; i < 9; ++i) { A[i] = i * 0.5 - 0.5; B[i] = i; @@ -1318,184 +1430,190 @@ TEST(LinalgBackendEigen, SGMatrix_multiply_by_rectified_linear_derivative) multiply_by_rectified_linear_derivative(A, B); for (index_t i = 0; i < 9; ++i) - EXPECT_NEAR(i * (A[i] != 0), B[i], 1e-15); + EXPECT_NEAR(i * (A[i] != 0), B[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_norm) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_norm) { - const index_t n = 5; - SGVector v(n); - float64_t gt = 0; + const index_t n = 24; + SGVector v(n); + TypeParam gt = 0; for (index_t i = 0; i < n; ++i) { v[i] = i; gt += i * i; } + gt = std::sqrt(gt); auto result = norm(v); - EXPECT_NEAR(result, gt, 1E-15); + EXPECT_NEAR(result, gt, get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_qr_solver) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGVector_qr_solver) { const index_t n = 3; - float64_t data_A[] = {0.02800922, 0.99326012, 0.15204902, + TypeParam data_A[] = {0.02800922, 0.99326012, 0.15204902, 0.30492837, 0.39708534, 0.40466969, 0.36415317, 0.04407589, 0.9095746}; - float64_t data_b[] = {0.39461571, 0.6816856, 0.43323709}; - float64_t result[] = {0.07135206, 1.56393127, -0.23141312}; + TypeParam data_b[] = {0.39461571, 0.6816856, 0.43323709}; + TypeParam result[] = {0.07135206, 1.56393127, -0.23141312}; - SGMatrix A(data_A, n, n, false); - SGVector b(data_b, n, false); + SGMatrix A(data_A, n, n, false); + SGVector b(data_b, n, false); auto x = qr_solver(A, b); for (index_t i = 0; i < x.size(); ++i) - EXPECT_NEAR(x[i], result[i], 1E-7); + EXPECT_NEAR(x[i], result[i], 1e-6); } -TEST(LinalgBackendEigen, SGMatrix_qr_solver) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_qr_solver) { const index_t n = 3, m = 2; - float64_t data_A[] = {0.02800922, 0.99326012, 0.15204902, + TypeParam data_A[] = {0.02800922, 0.99326012, 0.15204902, 0.30492837, 0.39708534, 0.40466969, 0.36415317, 0.04407589, 0.9095746}; - float64_t data_B[] = {0.76775073, 0.88471312, 0.34795225, + TypeParam data_B[] = {0.76775073, 0.88471312, 0.34795225, 0.94311546, 0.59630347, 0.65820143}; - float64_t result[] = {-0.73834587, 4.22750496, -1.37484721, + TypeParam result[] = {-0.73834587, 4.22750496, -1.37484721, -1.14718091, 4.49142548, -1.08282992}; - SGMatrix A(data_A, n, n, false); - SGMatrix B(data_B, n, m, false); + SGMatrix A(data_A, n, n, false); + SGMatrix B(data_B, n, m, false); auto X = qr_solver(A, B); for (index_t i = 0; i < (index_t)X.size(); ++i) - EXPECT_NEAR(X[i], result[i], 1E-7); + EXPECT_NEAR(X[i], result[i], 1e-5); } -TEST(LinalgBackendEigen, SGVector_range_fill) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_range_fill) { const index_t size = 5; - SGVector vec(size); - range_fill(vec, 1); + SGVector vec(size); + TypeParam start = 1; + range_fill(vec, start); for (index_t i = 0; i < size; ++i) - EXPECT_NEAR(vec[i], i + 1, 1E-15); + EXPECT_NEAR(vec[i], i + 1, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_range_fill) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_range_fill) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); - range_fill(mat, 1); + SGMatrix mat(nrows, ncols); + TypeParam start = 1; + range_fill(mat, start); - for (index_t i = 0; i < nrows*ncols; ++i) - EXPECT_NEAR(mat[i], i + 1, 1E-15); + for (index_t i = 0; i < nrows * ncols; ++i) + EXPECT_NEAR(mat[i], i + 1, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_rectified_linear) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_rectified_linear) { - SGMatrix A(3, 3); - SGMatrix B(3, 3); - - range_fill(A, -5.0); + SGMatrix A(3, 3); + SGMatrix B(3, 3); + TypeParam start = 1; + range_fill(A, start); linalg::rectified_linear(A, B); for (index_t i = 0; i < 9; ++i) - EXPECT_NEAR(CMath::max(0.0, A[i]), B[i], 1e-15); + EXPECT_NEAR( + CMath::max(static_cast(0.0), A[i]), B[i], + get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_scale) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_scale) { const index_t size = 5; - const float64_t alpha = 0.3; - SGVector a(size); + const TypeParam alpha = 2; + SGVector a(size); a.range_fill(0); auto result = scale(a, alpha); for (index_t i = 0; i < size; ++i) - EXPECT_NEAR(alpha * a[i], result[i], 1e-15); + EXPECT_NEAR(alpha * a[i], result[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_scale) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_scale) { - const float64_t alpha = 0.3; + const TypeParam alpha = 2; const index_t nrows = 2, ncols = 3; - SGMatrix A(nrows, ncols); + SGMatrix A(nrows, ncols); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) A[i] = i; auto result = scale(A, alpha); - for (index_t i = 0; i < nrows*ncols; ++i) - EXPECT_NEAR(alpha*A[i], result[i], 1e-15); + for (index_t i = 0; i < nrows * ncols; ++i) + EXPECT_NEAR(alpha * A[i], result[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_scale_in_place) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_scale_in_place) { const index_t size = 5; - const float64_t alpha = 0.3; - SGVector a(size); + const TypeParam alpha = 2; + SGVector a(size); a.range_fill(0); scale(a, a, alpha); for (index_t i = 0; i < size; ++i) - EXPECT_NEAR(alpha * i, a[i], 1e-15); + EXPECT_NEAR(alpha * i, a[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_scale_in_place) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_scale_in_place) { - const float64_t alpha = 0.3; + const TypeParam alpha = 2; const index_t nrows = 2, ncols = 3; - SGMatrix A(nrows, ncols); + SGMatrix A(nrows, ncols); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) A[i] = i; scale(A, A, alpha); - for (index_t i = 0; i < nrows*ncols; ++i) - EXPECT_NEAR(alpha*i, A[i], 1e-15); + for (index_t i = 0; i < nrows * ncols; ++i) + EXPECT_NEAR(alpha * i, A[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_set_const) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_set_const) { const index_t size = 5; - const float64_t value = 2; - SGVector a(size); + const TypeParam value = 2; + SGVector a(size); set_const(a, value); for (index_t i = 0; i < size; ++i) - EXPECT_NEAR(a[i], value, 1E-15); + EXPECT_NEAR(a[i], value, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_set_const) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_set_const) { const index_t nrows = 2, ncols = 3; - const float64_t value = 2; - SGMatrix a(nrows, ncols); + const TypeParam value = 2; + SGMatrix a(nrows, ncols); set_const(a, value); - for (index_t i = 0; i < nrows*ncols; ++i) - EXPECT_NEAR(a[i], value, 1E-15); + for (index_t i = 0; i < nrows * ncols; ++i) + EXPECT_NEAR(a[i], value, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_softmax) +// TODO: extend to all types +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_softmax) { - SGMatrix A(4, 3); - SGMatrix ref(4, 3); + SGMatrix A(4, 3); + SGMatrix ref(4, 3); - for (float64_t i = 0; i < 12; ++i) + for (TypeParam i = 0; i < 12; ++i) A[i] = i / 12; for (index_t i = 0; i < 12; ++i) @@ -1503,7 +1621,7 @@ TEST(LinalgBackendEigen, SGMatrix_softmax) for (index_t j = 0; j < ref.num_cols; ++j) { - float64_t sum = 0; + TypeParam sum = 0; for (index_t i = 0; i < ref.num_rows; ++i) sum += ref(i, j); @@ -1514,72 +1632,51 @@ TEST(LinalgBackendEigen, SGMatrix_softmax) linalg::softmax(A); for (index_t i = 0; i < 12; ++i) - EXPECT_NEAR(ref[i], A[i], 1e-15); + EXPECT_NEAR(ref[i], A[i], get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_squared_error) -{ - SGMatrix A(4, 3); - SGMatrix B(4, 3); - - int32_t size = A.num_rows * A.num_cols; - for (float64_t i = 0; i < size; ++i) - { - A[i] = i / size; - B[i] = (i / size) * 0.5; - } - - float64_t ref = 0; - for (index_t i = 0; i < size; i++) - ref += CMath::pow(A[i] - B[i], 2); - ref *= 0.5; - - auto result = linalg::squared_error(A, B); - EXPECT_NEAR(ref, result, 1e-15); -} - -TEST(LinalgBackendEigen, SGVector_sum) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_sum) { const index_t size = 10; - SGVector vec(size); + SGVector vec(size); vec.range_fill(0); auto result = sum(vec); - EXPECT_NEAR(result, 45, 1E-15); + EXPECT_NEAR(result, 45, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_sum) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_sum) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows * ncols; ++i) mat[i] = i; auto result = sum(mat); - EXPECT_NEAR(result, 15, 1E-15); + EXPECT_NEAR(result, 15, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_sum_no_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_sum_no_diag) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows * ncols; ++i) mat[i] = i; auto result = sum(mat, true); - EXPECT_NEAR(result, 12, 1E-15); + EXPECT_NEAR(result, 12, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_symmetric_with_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_symmetric_with_diag) { const index_t n = 3; - SGMatrix mat(n, n); - mat.set_const(1.0); + SGMatrix mat(n, n); + mat.set_const(1); for (index_t i = 0; i < n; ++i) for (index_t j = i + 1; j < n; ++j) @@ -1588,14 +1685,14 @@ TEST(LinalgBackendEigen, SGMatrix_symmetric_with_diag) mat(j, i) = mat(i, j); } - EXPECT_NEAR(sum_symmetric(mat), 39.0, 1E-15); + EXPECT_NEAR(sum_symmetric(mat), 39, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_symmetric_no_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_symmetric_no_diag) { const index_t n = 3; - SGMatrix mat(n, n); - mat.set_const(1.0); + SGMatrix mat(n, n); + mat.set_const(1); for (index_t i = 0; i < n; ++i) for (index_t j = i + 1; j < n; ++j) @@ -1604,13 +1701,13 @@ TEST(LinalgBackendEigen, SGMatrix_symmetric_no_diag) mat(j, i) = mat(i, j); } - EXPECT_NEAR(sum_symmetric(mat, true), 36.0, 1E-15); + EXPECT_NEAR(sum_symmetric(mat, true), 36, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_symmetric_exception) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_symmetric_exception) { const index_t n = 3; - SGMatrix mat(n, n + 1); + SGMatrix mat(n, n + 1); mat.set_const(1.0); for (index_t i = 0; i < n; ++i) @@ -1623,24 +1720,24 @@ TEST(LinalgBackendEigen, SGMatrix_symmetric_exception) EXPECT_THROW(sum_symmetric(mat), ShogunException); } -TEST(LinalgBackendEigen, SGMatrix_block_sum) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_block_sum) { const index_t n = 3; - SGMatrix mat(n, n); + SGMatrix mat(n, n); for (index_t i = 0; i < n; ++i) for (index_t j = 0; j < n; ++j) - mat(i, j)=i * 10 + j + 1; + mat(i, j) = i * 10 + j + 1; auto result = sum(linalg::block(mat, 0, 0, 2, 3)); - EXPECT_NEAR(result, 42.0, 1E-15); + EXPECT_NEAR(result, 42.0, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_symmetric_block_with_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_symmetric_block_with_diag) { const index_t n = 3; - SGMatrix mat(n, n); - mat.set_const(1.0); + SGMatrix mat(n, n); + mat.set_const(1); for (index_t i = 0; i < n; ++i) for (index_t j = i + 1; j < n; ++j) @@ -1649,15 +1746,15 @@ TEST(LinalgBackendEigen, SGMatrix_symmetric_block_with_diag) mat(j, i) = mat(i, j); } - float64_t sum = sum_symmetric(linalg::block(mat,1,1,2,2)); - EXPECT_NEAR(sum, 28.0, 1E-15); + TypeParam sum = sum_symmetric(linalg::block(mat, 1, 1, 2, 2)); + EXPECT_NEAR(sum, 28, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_symmetric_block_no_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_symmetric_block_no_diag) { const index_t n = 3; - SGMatrix mat(n, n); - mat.set_const(1.0); + SGMatrix mat(n, n); + mat.set_const(1); for (index_t i = 0; i < n; ++i) for (index_t j = i + 1; j < n; ++j) @@ -1666,11 +1763,11 @@ TEST(LinalgBackendEigen, SGMatrix_symmetric_block_no_diag) mat(j, i) = mat(i, j); } - float64_t sum = sum_symmetric(linalg::block(mat,1,1,2,2), true); - EXPECT_NEAR(sum, 26.0, 1E-15); + TypeParam sum = sum_symmetric(linalg::block(mat, 1, 1, 2, 2), true); + EXPECT_NEAR(sum, 26, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_symmetric_block_exception) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_symmetric_block_exception) { const index_t n = 3; SGMatrix mat(n, n); @@ -1683,45 +1780,45 @@ TEST(LinalgBackendEigen, SGMatrix_symmetric_block_exception) mat(j, i) = mat(i, j); } - EXPECT_THROW(sum_symmetric(linalg::block(mat,1,1,2,3)), ShogunException); + EXPECT_THROW( + sum_symmetric(linalg::block(mat, 1, 1, 2, 3)), ShogunException); } -TEST(LinalgBackendEigen, SGMatrix_colwise_sum) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_colwise_sum) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows * ncols; ++i) mat[i] = i; - SGVector result = colwise_sum(mat); + SGVector result = colwise_sum(mat); for (index_t j = 0; j < ncols; ++j) { - int32_t sum = 0; + TypeParam sum = 0; for (index_t i = 0; i < nrows; ++i) sum += mat(i, j); - EXPECT_NEAR(sum, result[j], 1E-15); + EXPECT_NEAR(sum, result[j], get_epsilon()); } } -TEST(LinalgBackendEigen, SGMatrix_colwise_sum_no_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_colwise_sum_no_diag) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows * ncols; ++i) mat[i] = i; - SGVector result = colwise_sum(mat, true); + SGVector result = colwise_sum(mat, true); - EXPECT_NEAR(result[0], 1, 1E-15); - EXPECT_NEAR(result[1], 2, 1E-15); - EXPECT_NEAR(result[2], 9, 1E-15); + EXPECT_NEAR(result[0], 1, get_epsilon()); + EXPECT_NEAR(result[1], 2, get_epsilon()); + EXPECT_NEAR(result[2], 9, get_epsilon()); } - -TEST(LinalgBackendEigen, SGMatrix_block_colwise_sum) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_block_colwise_sum) { const index_t nrows = 2, ncols = 3; SGMatrix mat(nrows, ncols); @@ -1734,50 +1831,50 @@ TEST(LinalgBackendEigen, SGMatrix_block_colwise_sum) for (index_t j = 0; j < ncols; ++j) { - float64_t sum = 0; + TypeParam sum = 0; for (index_t i = 0; i < nrows; ++i) sum += mat(i, j); - EXPECT_NEAR(sum, result[j], 1E-15); + EXPECT_NEAR(sum, result[j], get_epsilon()); } } -TEST(LinalgBackendEigen, SGMatrix_rowwise_sum) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_rowwise_sum) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows * ncols; ++i) mat[i] = i; - SGVector result = rowwise_sum(mat); + SGVector result = rowwise_sum(mat); for (index_t i = 0; i < nrows; ++i) { - int32_t sum = 0; + TypeParam sum = 0; for (index_t j = 0; j < ncols; ++j) sum += mat(i, j); - EXPECT_NEAR(sum, result[i], 1E-15); + EXPECT_NEAR(sum, result[i], get_epsilon()); } } -TEST(LinalgBackendEigen, SGMatrix_rowwise_sum_no_diag) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_rowwise_sum_no_diag) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows * ncols; ++i) mat[i] = i; - SGVector result = rowwise_sum(mat, true); + SGVector result = rowwise_sum(mat, true); - EXPECT_NEAR(result[0], 6, 1E-15); - EXPECT_NEAR(result[1], 6, 1E-15); + EXPECT_NEAR(result[0], 6, get_epsilon()); + EXPECT_NEAR(result[1], 6, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_block_rowwise_sum) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_block_rowwise_sum) { const index_t nrows = 2, ncols = 3; - SGMatrix mat(nrows, ncols); + SGMatrix mat(nrows, ncols); for (index_t i = 0; i < nrows; ++i) for (index_t j = 0; j < ncols; ++j) @@ -1787,29 +1884,29 @@ TEST(LinalgBackendEigen, SGMatrix_block_rowwise_sum) for (index_t i = 0; i < nrows; ++i) { - float64_t sum = 0; + TypeParam sum = 0; for (index_t j = 0; j < ncols; ++j) sum += mat(i, j); - EXPECT_NEAR(sum, result[i], 1E-15); + EXPECT_NEAR(sum, result[i], get_epsilon()); } } -TEST(LinalgBackendEigen, SGMatrix_svd_jacobi_thinU) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_svd_jacobi_thinU) { const index_t m = 5, n = 3; - float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, + TypeParam data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, 0.30786772, 0.25503552, 0.34367041, 0.66491478, 0.20488809, 0.5734351, 0.87179189, 0.07139643, 0.28540373, 0.06264684, 0.56204061}; - float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883}; - float64_t result_U[] = {-0.60700926, -0.16647013, -0.56501385, -0.26696629, + TypeParam result_s[] = {1.75382524, 0.56351367, 0.41124883}; + TypeParam result_U[] = {-0.60700926, -0.16647013, -0.56501385, -0.26696629, -0.46186125, -0.69145782, 0.29548428, 0.5718984, 0.31771648, -0.08101592, -0.27461424, 0.37170223, -0.12681555, -0.53830325, 0.69323293}; - SGMatrix A(data, m, n, false); - SGMatrix U(m, n); - SGVector s(n); + SGMatrix A(data, m, n, false); + SGMatrix U(m, n); + SGVector s(n); svd(A, s, U, true, SVDAlgorithm::Jacobi); @@ -1817,30 +1914,30 @@ TEST(LinalgBackendEigen, SGMatrix_svd_jacobi_thinU) { auto c = CMath::sign(U[i * m] * result_U[i * m]); for (index_t j = 0; j < m; ++j) - EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7); + EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-6); } for (index_t i = 0; i < (index_t)s.size(); ++i) - EXPECT_NEAR(s[i], result_s[i], 1e-7); + EXPECT_NEAR(s[i], result_s[i], 1e-6); } -TEST(LinalgBackendEigen, SGMatrix_svd_jacobi_fullU) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_svd_jacobi_fullU) { const index_t m = 5, n = 3; - float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, + TypeParam data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, 0.30786772, 0.25503552, 0.34367041, 0.66491478, 0.20488809, 0.5734351, 0.87179189, 0.07139643, 0.28540373, 0.06264684, 0.56204061}; - float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883}; - float64_t result_U[] = { + TypeParam result_s[] = {1.75382524, 0.56351367, 0.41124883}; + TypeParam result_U[] = { -0.60700926, -0.16647013, -0.56501385, -0.26696629, -0.46186125, -0.69145782, 0.29548428, 0.5718984, 0.31771648, -0.08101592, -0.27461424, 0.37170223, -0.12681555, -0.53830325, 0.69323293, -0.27809756, -0.68975171, -0.11662812, 0.38274703, 0.53554354, 0.025973184, 0.520631112, -0.56921636, 0.62571522, 0.11287970}; - SGMatrix A(data, m, n, false); - SGMatrix U(m, m); - SGVector s(n); + SGMatrix A(data, m, n, false); + SGMatrix U(m, m); + SGVector s(n); svd(A, s, U, false, SVDAlgorithm::Jacobi); @@ -1848,29 +1945,29 @@ TEST(LinalgBackendEigen, SGMatrix_svd_jacobi_fullU) { auto c = CMath::sign(U[i * m] * result_U[i * m]); for (index_t j = 0; j < m; ++j) - EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7); + EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-6); } for (index_t i = 0; i < (index_t)s.size(); ++i) - EXPECT_NEAR(s[i], result_s[i], 1e-7); + EXPECT_NEAR(s[i], result_s[i], 1e-6); } #if EIGEN_VERSION_AT_LEAST(3, 3, 0) -TEST(LinalgBackendEigen, SGMatrix_svd_bdc_thinU) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_svd_bdc_thinU) { const index_t m = 5, n = 3; - float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, + TypeParam data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, 0.30786772, 0.25503552, 0.34367041, 0.66491478, 0.20488809, 0.5734351, 0.87179189, 0.07139643, 0.28540373, 0.06264684, 0.56204061}; - float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883}; - float64_t result_U[] = {-0.60700926, -0.16647013, -0.56501385, -0.26696629, + TypeParam result_s[] = {1.75382524, 0.56351367, 0.41124883}; + TypeParam result_U[] = {-0.60700926, -0.16647013, -0.56501385, -0.26696629, -0.46186125, -0.69145782, 0.29548428, 0.5718984, 0.31771648, -0.08101592, -0.27461424, 0.37170223, -0.12681555, -0.53830325, 0.69323293}; - SGMatrix A(data, m, n, false); - SGMatrix U(m, n); - SGVector s(n); + SGMatrix A(data, m, n, false); + SGMatrix U(m, n); + SGVector s(n); svd(A, s, U, true, SVDAlgorithm::BidiagonalDivideConquer); @@ -1878,30 +1975,30 @@ TEST(LinalgBackendEigen, SGMatrix_svd_bdc_thinU) { auto c = CMath::sign(U[i * m] * result_U[i * m]); for (index_t j = 0; j < m; ++j) - EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7); + EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-6); } for (index_t i = 0; i < (index_t)s.size(); ++i) - EXPECT_NEAR(s[i], result_s[i], 1e-7); + EXPECT_NEAR(s[i], result_s[i], 1e-6); } -TEST(LinalgBackendEigen, SGMatrix_svd_bdc_fullU) +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_svd_bdc_fullU) { const index_t m = 5, n = 3; - float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, + TypeParam data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, 0.30786772, 0.25503552, 0.34367041, 0.66491478, 0.20488809, 0.5734351, 0.87179189, 0.07139643, 0.28540373, 0.06264684, 0.56204061}; - float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883}; - float64_t result_U[] = { + TypeParam result_s[] = {1.75382524, 0.56351367, 0.41124883}; + TypeParam result_U[] = { -0.60700926, -0.16647013, -0.56501385, -0.26696629, -0.46186125, -0.69145782, 0.29548428, 0.5718984, 0.31771648, -0.08101592, -0.27461424, 0.37170223, -0.12681555, -0.53830325, 0.69323293, -0.27809756, -0.68975171, -0.11662812, 0.38274703, 0.53554354, 0.025973184, 0.520631112, -0.56921636, 0.62571522, 0.11287970}; - SGMatrix A(data, m, n, false); - SGMatrix U(m, m); - SGVector s(n); + SGMatrix A(data, m, n, false); + SGMatrix U(m, m); + SGVector s(n); svd(A, s, U, false, SVDAlgorithm::BidiagonalDivideConquer); @@ -1909,165 +2006,171 @@ TEST(LinalgBackendEigen, SGMatrix_svd_bdc_fullU) { auto c = CMath::sign(U[i * m] * result_U[i * m]); for (index_t j = 0; j < m; ++j) - EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7); + EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-6); } for (index_t i = 0; i < (index_t)s.size(); ++i) - EXPECT_NEAR(s[i], result_s[i], 1e-7); + EXPECT_NEAR(s[i], result_s[i], 1e-6); } #endif -TEST(LinalgBackendEigen, SGMatrix_trace) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_trace) { const index_t n = 4; - SGMatrix A(n, n); - for (index_t i = 0; i < n*n; ++i) + SGMatrix A(n, n); + for (index_t i = 0; i < n * n; ++i) A[i] = i; - float64_t tr = 0; + TypeParam tr = 0; for (index_t i = 0; i < n; ++i) tr += A.get_element(i, i); - EXPECT_NEAR(trace(A), tr, 1e-15); + EXPECT_NEAR(trace(A), tr, get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_trace_dot) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_trace_dot) { - const index_t m = 2; - float64_t data_A[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194}; - float64_t data_B[] = {0.30786772, 0.25503552, 0.34367041, 0.66491478}; - - SGMatrix A(data_A, m, m, false); - SGMatrix B(data_B, m, m, false); + const index_t n = 2; + SGMatrix A(n, n), B(n, n); + for (index_t i = 0; i < n * n; ++i) + { + A[i] = i; + B[i] = i * 2; + } auto C = matrix_prod(A, B); auto tr = 0.0; - for (auto i : range(m)) + for (auto i : range(n)) tr += C(i, i); - EXPECT_NEAR(tr, trace_dot(A, B), 1e-15); + EXPECT_NEAR(tr, trace_dot(A, B), get_epsilon()); } -TEST(LinalgBackendEigen, SGMatrix_transpose_matrix) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_transpose_matrix) { const index_t m = 5, n = 3; - float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194, - 0.30786772, 0.25503552, 0.34367041, 0.66491478, - 0.20488809, 0.5734351, 0.87179189, 0.07139643, - 0.28540373, 0.06264684, 0.56204061}; - - SGMatrix A(data, m, n, false); + SGMatrix A(m, n); + linalg::range_fill(A, static_cast(1)); auto T = transpose_matrix(A); for (index_t i = 0; i < m; ++i) for (index_t j = 0; j < n; ++j) - EXPECT_NEAR(A.get_element(i, j), T.get_element(j, i), 1e-15); + EXPECT_NEAR( + A.get_element(i, j), T.get_element(j, i), + get_epsilon()); } -TEST(LinalgBackendEigen, SGVector_triangular_solver_lower) +TYPED_TEST( + LinalgBackendEigenNonIntegerTypesTest, SGVector_triangular_solver_lower) { const index_t n = 3; - float64_t data_L[] = {-0.92947874, -1.1432887, -0.87119086, + TypeParam data_L[] = {-0.92947874, -1.1432887, -0.87119086, 0., -0.27048649, -0.05919915, 0., 0., 0.11869106}; - float64_t data_b[] = {0.39461571, 0.6816856, 0.43323709}; - float64_t result[] = {-0.42455592, -0.72571316, 0.17192745}; + TypeParam data_b[] = {0.39461571, 0.6816856, 0.43323709}; + TypeParam result[] = {-0.42455592, -0.72571316, 0.17192745}; - SGMatrix L(data_L, n, n, false); - SGVector b(data_b, n, false); + SGMatrix L(data_L, n, n, false); + SGVector b(data_b, n, false); auto x = triangular_solver(L, b, true); for (index_t i = 0; i < (index_t)x.size(); ++i) - EXPECT_NEAR(x[i], result[i], 1E-6); + EXPECT_NEAR(x[i], result[i], 1e-6); } -TEST(LinalgBackendEigen, SGVector_triangular_solver_upper) +TYPED_TEST( + LinalgBackendEigenNonIntegerTypesTest, SGVector_triangular_solver_upper) { const index_t n = 3; - float64_t data_U[] = {-0.92947874, 0., 0., + TypeParam data_U[] = {-0.92947874, 0., 0., -1.1432887, -0.27048649, 0., -0.87119086, -0.05919915, 0.11869106}; - float64_t data_b[] = {0.39461571, 0.6816856, 0.43323709}; - float64_t result[] = {0.23681135, -3.31909306, 3.65012412}; + TypeParam data_b[] = {0.39461571, 0.6816856, 0.43323709}; + TypeParam result[] = {0.23681135, -3.31909306, 3.65012412}; - SGMatrix U(data_U, n, n, false); - SGVector b(data_b, n, false); + SGMatrix U(data_U, n, n, false); + SGVector b(data_b, n, false); auto x = triangular_solver(U, b, false); for (index_t i = 0; i < (index_t)x.size(); ++i) - EXPECT_NEAR(x[i], result[i], 1E-6); + EXPECT_NEAR(x[i], result[i], 1e-6); } -TEST(LinalgBackendEigen, SGMatrix_triangular_solver_lower) +TYPED_TEST( + LinalgBackendEigenNonIntegerTypesTest, SGMatrix_triangular_solver_lower) { const index_t n = 3, m = 2; - float64_t data_L[] = {-0.92947874, -1.1432887, -0.87119086, + TypeParam data_L[] = {-0.92947874, -1.1432887, -0.87119086, 0., -0.27048649, -0.05919915, 0., 0., 0.11869106}; - float64_t data_B[] = {0.76775073, 0.88471312, 0.34795225, + TypeParam data_B[] = {0.76775073, 0.88471312, 0.34795225, 0.94311546, 0.59630347, 0.65820143}; - float64_t result[] = {-0.82600139, 0.22050986, -3.02127745, + TypeParam result[] = {-0.82600139, 0.22050986, -3.02127745, -1.01467136, 2.08424024, -0.86262387}; - SGMatrix L(data_L, n, n, false); - SGMatrix B(data_B, n, m, false); + SGMatrix L(data_L, n, n, false); + SGMatrix B(data_B, n, m, false); auto X = triangular_solver(L, B, true); for (index_t i = 0; i < (index_t)X.size(); ++i) - EXPECT_NEAR(X[i], result[i], 1E-6); + EXPECT_NEAR(X[i], result[i], 1e-6); } -TEST(LinalgBackendEigen, SGMatrix_triangular_solver_upper) +TYPED_TEST( + LinalgBackendEigenNonIntegerTypesTest, SGMatrix_triangular_solver_upper) { const index_t n = 3, m = 2; - float64_t data_U[] = {-0.92947874, 0., 0., + TypeParam data_U[] = {-0.92947874, 0., 0., -1.1432887, -0.27048649, 0., -0.87119086, -0.05919915, 0.11869106}; - float64_t data_B[] = {0.76775073, 0.88471312, 0.34795225, + TypeParam data_B[] = {0.76775073, 0.88471312, 0.34795225, 0.94311546, 0.59630347, 0.65820143}; - float64_t result[] = {1.238677, -3.91243241, 2.9315793, + TypeParam result[] = {1.238677, -3.91243241, 2.9315793, -2.00784647, -3.41825732, 5.54550138}; - SGMatrix L(data_U, n, n, false); - SGMatrix B(data_B, n, m, false); + SGMatrix L(data_U, n, n, false); + SGMatrix B(data_B, n, m, false); auto X = triangular_solver(L, B, false); for (index_t i = 0; i < (index_t)X.size(); ++i) - EXPECT_NEAR(X[i], result[i], 1E-6); + EXPECT_NEAR(X[i], result[i], 1e-6); } -TEST(LinalgBackendEigen, SGVector_zero) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGVector_zero) { const index_t n = 16; - SGVector a(n); + SGVector a(n); zero(a); for (index_t i = 0; i < n; ++i) EXPECT_EQ(a[i], 0); } -TEST(LinalgBackendEigen, SGMatrix_zero) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_zero) { const index_t nrows = 3, ncols = 4; - SGMatrix A(nrows, ncols); + SGMatrix A(nrows, ncols); zero(A); - for (index_t i = 0; i < nrows*ncols; ++i) + for (index_t i = 0; i < nrows * ncols; ++i) EXPECT_EQ(A[i], 0); } -TEST(LinalgBackendEigen, SGMatrix_rank_update) +TYPED_TEST(LinalgBackendEigenAllTypesTest, SGMatrix_rank_update) { + typedef Matrix Mxx; + typedef Matrix Vxd; + const index_t size = 2; - SGMatrix A(size, size); - SGVector b(size); - Map A_eig(A.matrix, size, size); - Map b_eig(b.vector, size); + SGMatrix A(size, size); + SGVector b(size); + Map A_eig(A.matrix, size, size); + Map b_eig(b.vector, size); A(0, 0) = 2.0; A(1, 0) = 1.0; @@ -2077,15 +2180,36 @@ TEST(LinalgBackendEigen, SGMatrix_rank_update) b[1] = 3; auto A2 = A.clone(); - Map A2_eig(A2.matrix, size, size); + Map A2_eig(A2.matrix, size, size); A2(0, 0) += b[0] * b[0]; A2(0, 1) += b[0] * b[1]; A2(1, 0) += b[1] * b[0]; A2(1, 1) += b[1] * b[1]; - rank_update(A, b); - EXPECT_NEAR((A2_eig - A_eig).norm(), 0.0, 1e-14); + rank_update(A, b, static_cast(1)); + EXPECT_NEAR((A2_eig - A_eig).norm(), 0, get_epsilon()); + + rank_update(A, b, static_cast(-1)); + EXPECT_NEAR((A_eig - A_eig).norm(), 0, get_epsilon()); +} + +TYPED_TEST(LinalgBackendEigenNonIntegerTypesTest, SGMatrix_squared_error) +{ + SGMatrix A(4, 3); + SGMatrix B(4, 3); - rank_update(A, b, -1.); - EXPECT_NEAR((A_eig - A_eig).norm(), 0.0, 1e-14); + int32_t size = A.num_rows * A.num_cols; + for (float64_t i = 0; i < size; ++i) + { + A[i] = i / size; + B[i] = (i / size) * 0.5; + } + + float64_t ref = 0; + for (index_t i = 0; i < size; i++) + ref += std::pow(A[i] - B[i], 2); + ref *= 0.5; + + auto result = linalg::squared_error(A, B); + EXPECT_NEAR(ref, result, get_epsilon()); }