diff --git a/test/basis.cxx b/test/basis.cxx index 1de96c3..c521832 100644 --- a/test/basis.cxx +++ b/test/basis.cxx @@ -5,10 +5,13 @@ #include #include +#include // for Approx #include #include +using Catch::Approx; + TEST_CASE("FiniteTempBasis consistency tests", "[basis]") { SECTION("Basic consistency") { double beta = 1.0; @@ -112,6 +115,91 @@ TEST_CASE("FiniteTempBasis consistency tests", "[basis]") { auto basis = sparseir::FiniteTempBasis(beta, omega_max, epsilon, kernel); auto sve = basis.sve_result; auto s = sve->s; + //REQUIRE(s.size() == 32); + + std::vector s_ref = {0.5242807065966564, 0.361040299707525, 0.1600617039313169, 0.06192139783088188, 0.019641646995563183, 0.005321140031657106, 0.001245435134907047, 0.0002553808249508306, 4.6445392784931696e-5, 7.57389586542119e-6, 1.1180101601552092e-6, 1.506251988966008e-7, 1.8653991892840962e-8, 2.136773728637427e-9, 2.276179221544401e-10, 2.2655690134240947e-11, 2.115880115422964e-12, 1.861108037178489e-13, 1.5466716841180263e-14, 1.218212516630768e-15, 5.590657287253601e-16, 4.656548094642833e-16, 4.552528808131262e-16, 4.341440592462354e-16, 3.744993780121804e-16, 3.549006072192367e-16, 3.277985748785467e-16, 3.2621304578629284e-16, 3.2046691517654354e-16, 3.097729851576022e-16, 2.4973730182025644e-16, 2.476022474231314e-16}; + Eigen::VectorXd s_double = s.template cast(); + REQUIRE(std::fabs(s_double[0] - s_ref[0]) < 1e-10); + REQUIRE(std::fabs(s_double[1] - s_ref[1]) < 1e-10); + REQUIRE(std::fabs(s_double[2] - s_ref[2]) < 1e-10); + REQUIRE(std::fabs(s_double[3] - s_ref[3]) < 1e-10); + REQUIRE(std::fabs(s_double[4] - s_ref[4]) < 1e-10); + REQUIRE(std::fabs(s_double[5] - s_ref[5]) < 1e-10); + REQUIRE(std::fabs(s_double[6] - s_ref[6]) < 1e-10); + REQUIRE(std::fabs(s_double[7] - s_ref[7]) < 1e-10); + REQUIRE(std::fabs(s_double[8] - s_ref[8]) < 1e-10); + REQUIRE(std::fabs(s_double[9] - s_ref[9]) < 1e-10); + REQUIRE(std::fabs(s_double[10] - s_ref[10]) < 1e-10); + REQUIRE(std::fabs(s_double[11] - s_ref[11]) < 1e-10); + REQUIRE(std::fabs(s_double[12] - s_ref[12]) < 1e-10); + REQUIRE(std::fabs(s_double[13] - s_ref[13]) < 1e-10); + REQUIRE(std::fabs(s_double[14] - s_ref[14]) < 1e-10); + REQUIRE(std::fabs(s_double[15] - s_ref[15]) < 1e-10); + REQUIRE(std::fabs(s_double[16] - s_ref[16]) < 1e-10); + REQUIRE(std::fabs(s_double[17] - s_ref[17]) < 1e-10); + REQUIRE(std::fabs(s_double[18] - s_ref[18]) < 1e-10); + REQUIRE(std::fabs(s_double[19] - s_ref[19]) < 1e-10); + //REQUIRE(std::fabs(s_double[20] - s_ref[20]) < 1e-10); + // REQUIRE(std::fabs(s_double[21] - s_ref[21]) < 1e-10); + // REQUIRE(std::fabs(s_double[22] - s_ref[22]) < 1e-10); + + REQUIRE(sve->u[0].data.rows() == 10); + REQUIRE(sve->u[0].data.cols() == 32); + + std::vector u_knots_ref = { + -1.0, -0.9768276289532026, -0.9502121116288913, -0.9196860690044226, + -0.8847415486995369, -0.8448386704449569, -0.7994218020611462, -0.7479461808659303, + -0.6899180675604202, -0.6249508554943133, -0.552837354044473, -0.4736340017820308, + -0.38774586460365346, -0.2959932285976203, -0.19963497739688743, -0.10032604651986517, + 0.0, 0.10032604651986517, 0.19963497739688743, 0.2959932285976203, + 0.38774586460365346, 0.4736340017820308, 0.552837354044473, 0.6249508554943133, + 0.6899180675604202, 0.7479461808659303, 0.7994218020611462, 0.8448386704449569, + 0.8847415486995369, 0.9196860690044226, 0.9502121116288913, 0.9768276289532026, + 1.0 + }; + Eigen::VectorXd u_knots_ref_eigen = Eigen::Map(u_knots_ref.data(), u_knots_ref.size()); + REQUIRE(sve->u[0].knots.isApprox(u_knots_ref_eigen)); + + std::vector v_knots_ref = { + -1.0, -0.9833147686254275, -0.9470082310185116, -0.8938959515018162, -0.8283053538395936, + -0.7548706158857138, + -0.6778753393916265, -0.600858151891138, -0.5264593296556019, -0.45644270870032133, + -0.39184906182935686, -0.3331494756803358, -0.2804096832724622, -0.23343248554812435, + -0.19185635090170117, -0.15524305920516734, -0.12312152382089525, + -0.0950206581112576, -0.070491286445028, -0.04911709970058231, -0.03050369976269751, + -0.014178372359576086, 0.0, 0.014178372359576086, 0.03050369976269751, 0.04911709970058231, + 0.070491286445028, 0.0950206581112576, 0.12312152382089525, 0.15524305920516734, 0.19185635090170117, + 0.23343248554812435, 0.2804096832724622, 0.3331494756803358, 0.39184906182935686, 0.45644270870032133, + 0.5264593296556019, 0.600858151891138, 0.6778753393916265, 0.7548706158857138, 0.8283053538395936, + 0.8938959515018162, 0.9470082310185116, 0.9833147686254275, 1.0 + }; + Eigen::VectorXd v_knots_ref_eigen = Eigen::Map(v_knots_ref.data(), v_knots_ref.size()); + REQUIRE(sve->v[0].knots.isApprox(v_knots_ref_eigen)); + + REQUIRE(sve->u[1].xmin == -1.0); + REQUIRE(sve->u[1].xmax == 1.0); + + REQUIRE(sve->v[1].xmin == -1.0); + REQUIRE(sve->v[1].xmax == 1.0); + + REQUIRE(sve->u[0].l == 0); + REQUIRE(sve->u[1].l == 1); + REQUIRE(sve->u[2].l == 2); + + REQUIRE(sve->v[0].l == 0); + REQUIRE(sve->v[1].l == 1); + REQUIRE(sve->v[2].l == 2); + + //REQUIRE(sve->u[1].symm == 1); + //REQUIRE(sve->u[2].symm == -1); + //REQUIRE(sve->u[3].symm == 1); + //REQUIRE(sve->u[4].symm == -1); + + //REQUIRE(sve->v[1].symm == 1); + //REQUIRE(sve->v[2].symm == -1); + //REQUIRE(sve->v[3].symm == 1); + //REQUIRE(sve->v[4].symm == -1); + //std::cout << "Singular values: " << s.transpose() << std::endl; /* int L = 10; @@ -121,5 +209,4 @@ TEST_CASE("FiniteTempBasis consistency tests", "[basis]") { REQUIRE(pts_100.size() == 24); */ } - }