From a8a9fb6eb8b3790b84ad4e798fa838958ffe32e5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 16 Sep 2020 18:17:23 -0400 Subject: [PATCH] adapt unit test for Jacobi class --- unittest/utils/CMakeLists.txt | 4 + unittest/utils/test_math_eigen_impl.cpp | 766 ++++++++++++++++++++++++ 2 files changed, 770 insertions(+) create mode 100644 unittest/utils/test_math_eigen_impl.cpp diff --git a/unittest/utils/CMakeLists.txt b/unittest/utils/CMakeLists.txt index 5b5b931210..1a10613403 100644 --- a/unittest/utils/CMakeLists.txt +++ b/unittest/utils/CMakeLists.txt @@ -14,3 +14,7 @@ set_tests_properties(Utils PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_PO add_executable(test_fmtlib test_fmtlib.cpp) target_link_libraries(test_fmtlib PRIVATE lammps GTest::GMockMain GTest::GMock GTest::GTest) add_test(FmtLib test_fmtlib) + +add_executable(test_math_eigen_impl test_math_eigen_impl.cpp) +target_include_directories(test_math_eigen_impl PRIVATE ${LAMMPS_SOURCE_DIR}) +add_test(MathEigen test_math_eigen_impl 10 5) diff --git a/unittest/utils/test_math_eigen_impl.cpp b/unittest/utils/test_math_eigen_impl.cpp new file mode 100644 index 0000000000..a39a71417d --- /dev/null +++ b/unittest/utils/test_math_eigen_impl.cpp @@ -0,0 +1,766 @@ +// THIS FILE USED TO BE EASY TO READ until I added "#if defined" statements. +// (They were added to test for many different kinds of array formats.) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "math_eigen_impl.h" + +using std::cout; +using std::cerr; +using std::endl; +using std::setprecision; +using std::vector; +using std::array; +using namespace MathEigen; + + +// This code works with various types of C++ matrices (for example, +// double **, vector> array,5>). +// I use "#if defined" statements to test different matrix types. +// For some of these (eg. array,5>), the size of the matrix +// must be known at compile time. I specify that size now. +#if defined USE_ARRAY_OF_ARRAYS +const int NF=5; //(the array size must be known at compile time) +#elif defined USE_C_FIXED_SIZE_ARRAYS +const int NF=5; //(the array size must be known at compile time) +#endif + + +// @brief Are two numbers "similar"? +template +inline static bool Similar(Scalar a, Scalar b, + Scalar eps=1.0e-06, + Scalar ratio=1.0e-06, + Scalar ratio_denom=1.0) +{ + return ((std::abs(a-b)<=std::abs(eps)) + || + (std::abs(ratio_denom)*std::abs(a-b) + <= + std::abs(ratio)*0.5*(std::abs(a)+std::abs(b)))); +} + +/// @brief Are two vectors (containing n numbers) similar? +template +inline static bool SimilarVec(Vector a, Vector b, int n, + Scalar eps=1.0e-06, + Scalar ratio=1.0e-06, + Scalar ratio_denom=1.0) +{ + for (int i = 0; i < n; i++) + if (not Similar(a[i], b[i], eps, ratio, ratio_denom)) + return false; + return true; +} + +/// @brief Are two vectors (or their reflections) similar? +template +inline static bool SimilarVecUnsigned(Vector a, Vector b, int n, + Scalar eps=1.0e-06, + Scalar ratio=1.0e-06, + Scalar ratio_denom=1.0) +{ + if (SimilarVec(a, b, n, eps)) + return true; + else { + for (int i = 0; i < n; i++) + if (not Similar(a[i], -b[i], eps, ratio, ratio_denom)) + return false; + return true; + } +} + + +/// @brief Multiply two matrices A and B, store the result in C. (C = AB). + +template +void mmult(ConstMatrix A, // +void +SortRows(Vector eval, + Matrix evec, + int n, + bool sort_decreasing=true, + bool sort_abs=false) +{ + for (int i = 0; i < n-1; i++) { + int i_max = i; + for (int j = i+1; j < n; j++) { + if (sort_decreasing) { + if (sort_abs) { //sort by absolute value? + if (std::abs(eval[j]) > std::abs(eval[i_max])) + i_max = j; + } + else if (eval[j] > eval[i_max]) + i_max = j; + } + else { + if (sort_abs) { //sort by absolute value? + if (std::abs(eval[j]) < std::abs(eval[i_max])) + i_max = j; + } + else if (eval[j] < eval[i_max]) + i_max = j; + } + } + std::swap(eval[i], eval[i_max]); // sort "eval" + for (int k = 0; k < n; k++) + std::swap(evec[i][k], evec[i_max][k]); // sort "evec" + } +} + + + +/// @brief Generate a random orthonormal n x n matrix + +template +void GenRandOrth(Matrix R, + int n, + std::default_random_engine &rand_generator) +{ + std::normal_distribution gaussian_distribution(0,1); + std::vector v(n); + + for (int i = 0; i < n; i++) { + // Generate a vector, "v", in a random direction subject to the constraint + // that it is orthogonal to the first i-1 rows-vectors of the R matrix. + Scalar rsq = 0.0; + while (rsq == 0.0) { + // Generate a vector in a random direction + // (This works because we are using a normal (Gaussian) distribution) + for (int j = 0; j < n; j++) + v[j] = gaussian_distribution(rand_generator); + + //Now subtract from v, the projection of v onto the first i-1 rows of R. + //This will produce a vector which is orthogonal to these i-1 row-vectors. + //(They are already normalized and orthogonal to each other.) + for (int k = 0; k < i; k++) { + Scalar v_dot_Rk = 0.0; + for (int j = 0; j < n; j++) + v_dot_Rk += v[j] * R[k][j]; + for (int j = 0; j < n; j++) + v[j] -= v_dot_Rk * R[k][j]; + } + // check if it is linearly independent of the other vectors and non-zero + rsq = 0.0; + for (int j = 0; j < n; j++) + rsq += v[j]*v[j]; + } + // Now normalize the vector + Scalar r_inv = 1.0 / std::sqrt(rsq); + for (int j = 0; j < n; j++) + v[j] *= r_inv; + // Now copy this vector to the i'th row of R + for (int j = 0; j < n; j++) + R[i][j] = v[j]; + } //for (int i = 0; i < n; i++) +} //void GenRandOrth() + + + +/// @brief Generate a random symmetric n x n matrix, M. +/// This function generates random numbers for the eigenvalues ("evals_known") +/// as well as the eigenvectors ("evecs_known"), and uses them to generate M. +/// The "eval_magnitude_range" argument specifies the the base-10 logarithm +/// of the range of eigenvalues desired. The "n_degeneracy" argument specifies +/// the number of repeated eigenvalues desired (if any). +/// @returns This function does not return a value. However after it is +/// invoked, the M matrix will be filled with random numbers. +/// Additionally, the "evals" and "evecs" arguments will contain +/// the eigenvalues and eigenvectors (one eigenvector per row) +/// of the matrix. Later, they can be compared with the eigenvalues +/// and eigenvectors calculated by Jacobi::Diagonalize() + +template +void GenRandSymm(Matrix M, // random_real01; + std::normal_distribution gaussian_distribution(0, max_eval_size); + bool use_log_uniform_distribution = false; + if (min_eval_size > 0.0) + use_log_uniform_distribution = true; + #if defined USE_VECTOR_OF_VECTORS + vector > D(n, vector(n)); + vector > tmp(n, vector(n)); + #elif defined USE_ARRAY_OF_ARRAYS + array, NF> D; + array, NF> tmp; + #elif defined USE_C_FIXED_SIZE_ARRAYS + Scalar D[NF][NF], tmp[NF][NF]; + #else + #define USE_C_POINTER_TO_POINTERS + Scalar **D, **tmp; + Alloc2D(n, n, &D); + Alloc2D(n, n, &tmp); + #endif + + // Randomly generate the eigenvalues + for (int i = 0; i < n; i++) { + if (use_log_uniform_distribution) { + // Use a "log-uniform distribution" (a.k.a. "reciprocal distribution") + // (This is a way to specify numbers with a precise range of magnitudes.) + assert((min_eval_size > 0.0) && (max_eval_size > 0.0)); + Scalar log_min = std::log(std::abs(min_eval_size)); + Scalar log_max = std::log(std::abs(max_eval_size)); + Scalar log_eval = (log_min + random_real01(rand_generator)*(log_max-log_min)); + evals[i] = std::exp(log_eval); + // also consider both positive and negative eigenvalues: + if (random_real01(rand_generator) < 0.5) + evals[i] = -evals[i]; + } + else { + evals[i] = gaussian_distribution(rand_generator); + } + } + + // Does the user want us to force some of the eigenvalues to be the same? + if (n_degeneracy > 1) { + int *permutation = new int[n]; //a random permutation from 0...n-1 + for (int i = 0; i < n; i++) + permutation[i] = i; + std::shuffle(permutation, permutation+n, rand_generator); + for (int i = 1; i < n_degeneracy; i++) //set the first n_degeneracy to same + evals[permutation[i]] = evals[permutation[0]]; + delete [] permutation; + } + + // D is a diagonal matrix whose diagonal elements are the eigenvalues + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) + D[i][j] = ((i == j) ? evals[i] : 0.0); + + // Now randomly generate the (transpose of) the "evecs" matrix + GenRandOrth(evecs, n, rand_generator); //(will transpose it later) + + // Construct the test matrix, M, where M = Rt * D * R + + // Original code: + //mmult(evecs, D, tmp, n); // <--> tmp = Rt * D + // Unfortunately, C++ guesses the types incorrectly. Must manually specify: + // #ifdefs making the code ugly again: + #if defined USE_VECTOR_OF_VECTORS + mmult >&, const vector >&> + #elif defined USE_ARRAY_OF_ARRAYS + mmult,NF>&, const array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + mmult + #else + mmult + #endif + (evecs, D, tmp, n); + + for (int i = 0; i < n-1; i++) + for (int j = i+1; j < n; j++) + std::swap(evecs[i][j], evecs[j][i]); //transpose "evecs" + + // Original code: + //mmult(tmp, evecs, M, n); + // Unfortunately, C++ guesses the types incorrectly. Must manually specify: + // #ifdefs making the code ugly again: + #if defined USE_VECTOR_OF_VECTORS + mmult >&, const vector >&> + #elif defined USE_ARRAY_OF_ARRAYS + mmult,NF>&, const array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + mmult + #else + mmult + #endif + (tmp, evecs, M, n); + //at this point M = Rt*D*R (where "R"="evecs") + + #if defined USE_C_POINTER_TO_POINTERS + Dealloc2D(&D); + Dealloc2D(&tmp); + #endif +} // GenRandSymm() + + + +template +void TestJacobi(int n, //&, + vector >&, + const vector >& > ecalc(n); + + // allocate the matrix, eigenvalues, eigenvectors + vector > M(n, vector(n)); + vector > evecs(n, vector(n)); + vector > evecs_known(n, vector(n)); + vector evals(n); + vector evals_known(n); + vector test_evec(n); + + #elif defined USE_ARRAY_OF_ARRAYS + + n = NF; + cout << "Testing std::array (fixed size).\n" + "(Ignoring first argument, and setting matrix size to " << n << ")" << endl; + + Jacobi&, + array, NF>&, + const array, NF>&> ecalc(n); + + // allocate the matrix, eigenvalues, eigenvectors + array, NF> M; + array, NF> evecs; + array, NF> evecs_known; + array evals; + array evals_known; + array test_evec; + + #elif defined USE_C_FIXED_SIZE_ARRAYS + + n = NF; + cout << "Testing C fixed size arrays.\n" + "(Ignoring first argument, and setting matrix size to " << n << ")" << endl; + Jacobi ecalc(n); + + // allocate the matrix, eigenvalues, eigenvectors + Scalar M[NF][NF]; + Scalar evecs[NF][NF]; + Scalar evecs_known[NF][NF]; + Scalar evals[NF]; + Scalar evals_known[NF]; + Scalar test_evec[NF]; + + #else + + #define USE_C_POINTER_TO_POINTERS + + // Note: Normally, you would just use this to instantiate Jacobi: + // Jacobi ecalc(n); + // ------------------------- + // ..but since Jacobi manages its own memory using new and delete, I also want + // to test that the copy constructors, copy operators, and destructors work. + // The following lines do this: + Jacobi ecalc_test_mem1(n); + Jacobi ecalc_test_mem2(2); + // test the = operator + ecalc_test_mem2 = ecalc_test_mem1; + // test the copy constructor + Jacobi ecalc(ecalc_test_mem2); + // allocate the matrix, eigenvalues, eigenvectors + Scalar **M, **evecs, **evecs_known; + Alloc2D(n, n, &M); + Alloc2D(n, n, &evecs); + Alloc2D(n, n, &evecs_known); + Scalar *evals = new Scalar[n]; + Scalar *evals_known = new Scalar[n]; + Scalar *test_evec = new Scalar[n]; + + #endif + + + // -------------------------------------------------------------------- + // Now, generate random matrices and test Jacobi::Diagonalize() on them. + // -------------------------------------------------------------------- + + for(int imat = 0; imat < n_matrices; imat++) { + + // Create a randomly generated symmetric matrix. + //This function generates random numbers for the eigenvalues ("evals_known") + //as well as the eigenvectors ("evecs_known"), and uses them to generate M. + + #if defined USE_VECTOR_OF_VECTORS + GenRandSymm&, vector >&> + #elif defined USE_ARRAY_OF_ARRAYS + GenRandSymm&, array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + GenRandSymm + #else + GenRandSymm + #endif + (M, + n, + evals_known, + evecs_known, + rand_generator, + min_eval_size, + max_eval_size, + n_degeneracy); + + // Sort the matrix evals and eigenvector rows: + // Original code: + //SortRows(evals_known, evecs_known, n); + // Unfortunately, C++ guesses the types incorrectly. Must use #ifdefs again: + #if defined USE_VECTOR_OF_VECTORS + SortRows&, vector >&> + #elif defined USE_ARRAY_OF_ARRAYS + SortRows&, array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + SortRows + #else + SortRows + #endif + (evals_known, evecs_known, n); + + + if (n_matrices == 1) { + cout << "Eigenvalues (after sorting):\n"; + for (int i = 0; i < n; i++) + cout << evals_known[i] << " "; + cout << "\n"; + cout << "Eigenvectors (rows) which are known in advance:\n"; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) + cout << evecs_known[i][j] << " "; + cout << "\n"; + } + cout << " (The eigenvectors calculated by Jacobi::Diagonalize() should match these.)\n"; + } + + for (int i_test = 0; i_test < n_tests_per_matrix; i_test++) { + + if (test_code_coverage) { + + // test SORT_INCREASING_ABS_EVALS: + #if defined USE_VECTOR_OF_VECTORS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + vector >&, + const vector >& >::SORT_INCREASING_ABS_EVALS); + #elif defined USE_ARRAY_OF_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + array,NF>&, + const array,NF>&>::SORT_INCREASING_ABS_EVALS); + #elif defined USE_C_FIXED_SIZE_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::SORT_INCREASING_ABS_EVALS); + #else + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::SORT_INCREASING_ABS_EVALS); + #endif + + for (int i = 1; i < n; i++) + assert(std::abs(evals[i-1])<=std::abs(evals[i])); + + // test SORT_DECREASING_ABS_EVALS: + #if defined USE_VECTOR_OF_VECTORS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + vector >&, + const vector >& >::SORT_DECREASING_ABS_EVALS); + #elif defined USE_ARRAY_OF_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + array,NF>&, + const array,NF>&>::SORT_DECREASING_ABS_EVALS); + #elif defined USE_C_FIXED_SIZE_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::SORT_DECREASING_ABS_EVALS); + #else + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::SORT_DECREASING_ABS_EVALS); + #endif + + for (int i = 1; i < n; i++) + assert(std::abs(evals[i-1])>=std::abs(evals[i])); + + // test SORT_INCREASING_EVALS: + #if defined USE_VECTOR_OF_VECTORS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + vector >&, + const vector >& >::SORT_INCREASING_EVALS); + #elif defined USE_ARRAY_OF_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + array,NF>&, + const array,NF>&>::SORT_INCREASING_EVALS); + #elif defined USE_C_FIXED_SIZE_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::SORT_INCREASING_EVALS); + #else + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::SORT_INCREASING_EVALS); + #endif + for (int i = 1; i < n; i++) + assert(evals[i-1] <= evals[i]); + + // test DO_NOT_SORT + #if defined USE_VECTOR_OF_VECTORS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + vector >&, + const vector >& >::DO_NOT_SORT); + #elif defined USE_ARRAY_OF_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi&, + array,NF>&, + const array,NF>&>::DO_NOT_SORT); + #elif defined USE_C_FIXED_SIZE_ARRAYS + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::DO_NOT_SORT); + #else + ecalc.Diagonalize(M, + evals, + evecs, + Jacobi::DO_NOT_SORT); + #endif + + } //if (test_code_coverage) + + + // Now (finally) calculate the eigenvalues and eigenvectors + int n_sweeps = ecalc.Diagonalize(M, evals, evecs); + + if ((n_matrices == 1) && (i_test == 0)) { + cout <<"Jacobi::Diagonalize() ran for "< Σ_b M[a][b]*evecs[i][b] = evals[i]*evecs[i][b] (for all a) + for (int i = 0; i < n; i++) { + for (int a = 0; a < n; a++) { + test_evec[a] = 0.0; + for (int b = 0; b < n; b++) + test_evec[a] += M[a][b] * evecs[i][b]; + assert(Similar(test_evec[a], + evals[i] * evecs[i][a], + eps, // tolerance (absolute difference) + eps*max_eval_size, // tolerance ratio (numerator) + evals_known[i] // tolerance ration (denominator) + )); + } + } + + } //for (int i_test = 0; i_test < n_tests_per_matrix; i++) + + } //for(int imat = 0; imat < n_matrices; imat++) { + + #if defined USE_C_POINTER_TO_POINTERS + Dealloc2D(&M); + Dealloc2D(&evecs); + Dealloc2D(&evecs_known); + delete [] evals; + delete [] evals_known; + delete [] test_evec; + #endif + +} //TestJacobi() + + +int main(int argc, char **argv) { + int n_size = 2; + int n_matr = 1; + double emin = 0.0; + double emax = 1.0; + int n_tests = 1; + int n_degeneracy = 1; + unsigned seed = 0; + + if (argc <= 1) { + cerr << + "Error: This program requires at least 1 argument.\n" + "\n" + "Description: Run Jacobi::Diagonalize() on randomly generated matrices.\n" + "\n" + "Arguments: n_size [n_matr emin emax n_degeneracy n_tests seed eps]\n" + " n_size = the size of the matrices\n" + " (NOTE: The remaining arguments are optional.)\n" + " n_matr = the number of randomly generated matrices to test\n" + " emin = the smallest possible eigenvalue magnitude (eg. 1e-05)\n" + " emax = the largest possible eigenvalue magnitude (>0 eg. 1e+05)\n" + " (NOTE: If emin=0, a normal distribution is used centered at 0.\n" + " Otherwise a log-uniform distribution is used from emin to emax.)\n" + " n_degeneracy = the number of repeated eigenvalues (1 disables, default)\n" + " n_tests = the number of times the eigenvalues and eigenvectors\n" + " are calculated for EACH matrix. By default this is 1.\n" + " (Increase this to at least 20 if you plan to use this\n" + " program for benchmarking (speed testing), because the time\n" + " needed for generating a random matrix is not negligible.)\n" + " (IF THIS NUMBER IS 0, it will test CODE-COVERAGE instead.)\n" + " seed = the seed used by the random number \"rand_generator\".\n" + " (If this number is 0, which is the default, the system\n" + " clock is used to choose a random seed.)\n" + " eps = the tolerance. The difference between eigenvalues and their\n" + " true value, cannot exceed this (multiplied by the eigenvalue\n" + " of maximum magnitude). Similarly, the difference between\n" + " the eigenvectors after multiplication by the matrix and by\n" + " and after multiplication by the eigenvalue, cannot exceed\n" + " eps*maximum_eigenvalue/eigenvalue. The default value is\n" + " 1.0e-06 (which works well for double precision numbers).\n" + << endl; + return 1; + } + + n_size = std::stoi(argv[1]); + if (argc > 2) + n_matr = std::stoi(argv[2]); + if (argc > 3) + emin = std::stof(argv[3]); + if (argc > 4) + emax = std::stof(argv[4]); + if (argc > 5) + n_degeneracy = std::stoi(argv[5]); + if (argc > 6) + n_tests = std::stoi(argv[6]); + if (argc > 7) + seed = std::stoi(argv[7]); + double eps = 1.0e-06; + if (argc > 8) + eps = std::stof(argv[8]); + + TestJacobi(n_size, n_matr, emin, emax, n_tests, n_degeneracy, seed, eps); + + cout << "test passed\n" << endl; + return EXIT_SUCCESS; +}