diff --git a/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp b/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp
index 200a679b637e2227d2dcc805cc55352806f58bbf..b7eef2bd641335d30e297759f2756aba654945e1 100644
--- a/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp
+++ b/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp
@@ -2,6 +2,7 @@
#include "Tests/GTestWrapper/google_test.h"
#include "Tests/Unit/Numeric/MultiQTest.h"
#include <ff/Face.h> // for diagnostic
+#include <iomanip>
//! Compare form factor for particle shapes A and B, where A is given special
//! parameter values so that it coincides with the more symmetric B.
@@ -10,27 +11,31 @@ class FFSpecializationTest : public testing::Test {
protected:
void run_test(IFormFactor* p0, IFormFactor* p1, double eps, double qmag1, double qmag2)
{
- formfactorTest::run_test_for_many_q([&](C3 q) { test_ff_eq(q, p0, p1, eps); }, qmag1,
- qmag2);
+ int failures = formfactorTest::run_test2_for_many_q(
+ [&](C3 q, bool report) -> double {
+ return test_ff_eq(q, p0, p1, eps, report); }, qmag1,
+ qmag2, eps);
+ EXPECT_EQ(failures, 0);
}
private:
- void test_ff_eq(C3 q, IFormFactor* p0, IFormFactor* p1, double eps)
+ double test_ff_eq(C3 q, IFormFactor* p0, IFormFactor* p1, double eps, bool report)
{
const complex_t f0 = p0->formfactor(q);
const complex_t f1 = p1->formfactor(q);
const double avge = (std::abs(f0) + std::abs(f1)) / 2;
- const double precision = std::max(1e-16, eps * avge);
- EXPECT_NEAR(real(f0), real(f1), precision) << "q=" << q << "\n"
-#ifdef ALGORITHM_DIAGNOSTIC
- << polyhedralDiagnosis.message() << "\n"
-#endif
- ;
- EXPECT_NEAR(imag(f0), imag(f1), precision) << "q=" << q << "\n"
-#ifdef ALGORITHM_DIAGNOSTIC
- << polyhedralDiagnosis.message() << "\n"
-#endif
- ;
+ const double abserr = std::max(fabs(real(f0) - real(f1)), fabs(imag(f0) - imag(f1)));
+ const double result = abserr / avge * std::min(1., eps * avge / 1e-16);
+ if (report) {
+ std::cout << "Deviation at q = " << q << ":\n";
+ std::cout << " Re(f0) = " << std::setprecision(16)
+ << real(f0) << ", Im(f0) = " << imag(f0) << "\n";
+ std::cout << " Re(f1) = " << real(f1) << ", Im(f1) = " << imag(f1) << "\n";
+ std::cout << " abs dev = " << std::setprecision(8)
+ << abserr << ", rel dev = " << abserr / avge
+ << ", score = " << result << ", limit = " << eps << "\n";
+ }
+ return result;
}
};
@@ -100,7 +105,7 @@ TEST_F(FFSpecializationTest, HorizontalCylinderSlicedVsUnsliced)
const double R = .3, L = 3;
HorizontalCylinder p0(R, L);
HorizontalCylinder p1(R, L, -R, +R * (1 - 3e-16));
- run_test(&p0, &p1, 4e-8, 1e-99, 100);
+ run_test(&p0, &p1, 3e-8, 1e-99, 100);
}
//*********** spheroids ***************
@@ -110,7 +115,7 @@ TEST_F(FFSpecializationTest, HemiEllipsoidAsTruncatedSphere)
const double R = 1.07;
HemiEllipsoid p0(R, R, R);
TruncatedSphere p1(R, R, 0);
- run_test(&p0, &p1, 1e-10, 1e-99, 5e2);
+ run_test(&p0, &p1, 6e-12, 1e-99, 5e2);
}
TEST_F(FFSpecializationTest, EllipsoidalCylinderAsCylinder)
@@ -118,7 +123,7 @@ TEST_F(FFSpecializationTest, EllipsoidalCylinderAsCylinder)
const double R = .8, H = 1.2;
EllipsoidalCylinder p0(R, R, H);
Cylinder p1(R, H);
- run_test(&p0, &p1, 1e-11, 1e-99, 50);
+ run_test(&p0, &p1, 5e-12, 1e-99, 50);
}
TEST_F(FFSpecializationTest, TruncatedSphereAsSphere)
@@ -126,7 +131,7 @@ TEST_F(FFSpecializationTest, TruncatedSphereAsSphere)
const double R = 1.;
TruncatedSphere p0(R, 2 * R, 0);
Sphere p1(R);
- run_test(&p0, &p1, 1e-11, .02, 5e1);
+ run_test(&p0, &p1, 1e-12, .02, 5e1);
}
TEST_F(FFSpecializationTest, SpheroidAsSphere)
diff --git a/Tests/Unit/Numeric/MultiQTest.cpp b/Tests/Unit/Numeric/MultiQTest.cpp
index d6bf03804f0ba4adb2f31f77d29410bb5be01cd2..61cc6aa0769739bf3752350daff30510410476a1 100644
--- a/Tests/Unit/Numeric/MultiQTest.cpp
+++ b/Tests/Unit/Numeric/MultiQTest.cpp
@@ -39,4 +39,40 @@ void run_test_for_many_q(std::function<void(C3)> run_one_test, double qmag_min,
}
}
+int run_test2_for_many_q(std::function<double(C3, bool)> run_one_test, double qmag_min,
+ double qmag_max, double eps, bool real_only)
+{
+ ParamGenerator<std::tuple<C3, C3, double, double, double>> gen = qlist;
+ double max_deviation = 0;
+ C3 q_at_max;
+ int failures = 0;
+ for (auto it : gen) {
+ const C3 q_maindir = std::get<0>(it);
+ const C3 q_sidedir = std::get<1>(it);
+ const double qrealmag = std::get<2>(it);
+ const double qsidemag = std::get<3>(it);
+ const double qimagrel = std::get<4>(it);
+ if (real_only && qimagrel)
+ continue;
+ const complex_t qmag(qrealmag, qrealmag * qimagrel);
+ if (std::abs(qmag) <= qmag_min || std::abs(qmag) >= qmag_max)
+ continue;
+ if (q_maindir == q_sidedir)
+ continue;
+ const C3 q = qmag * (q_maindir + qsidemag * q_sidedir).unit();
+
+ double deviation = run_one_test(q, false); // callback passed as argument
+ if (deviation > eps) {
+ ++failures;
+ if (deviation > max_deviation) {
+ max_deviation = deviation;
+ q_at_max = q;
+ }
+ }
+ }
+ if (failures)
+ run_one_test(q_at_max, true); // report worst case
+ return failures;
+}
+
} // namespace formfactorTest
diff --git a/Tests/Unit/Numeric/MultiQTest.h b/Tests/Unit/Numeric/MultiQTest.h
index 3bfad52ef2e478a078e2226d5c68211458d0fd57..152d2999a874d9300ee82ade7b80984bda427527 100644
--- a/Tests/Unit/Numeric/MultiQTest.h
+++ b/Tests/Unit/Numeric/MultiQTest.h
@@ -14,6 +14,10 @@ namespace formfactorTest {
void run_test_for_many_q(std::function<void(C3)> run_one_test, double qmag_min, double qmag_max,
bool real_only = false);
+//! Runs callback function "run_one_test(q)" for a huge number of different vectors q.
+int run_test2_for_many_q(std::function<double(C3, bool)> run_one_test, double qmag_min,
+ double qmag_max, double eps, bool real_only = false);
+
} // namespace formfactorTest