1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
#include "test/util/RunMultiQ.h"
#include <iomanip>
#include <iostream>
namespace formfactorTest {
const auto maindirs = {C3({1, 0, 0}), C3({0, 1, 0}), C3({0, 0, 1}), C3({1, 1, 0}),
C3({1, 0, 1}), C3({0, 1, 1}), C3({1, 1, 1})};
const auto sidedirs = {C3({1, 0, 0}), C3({0, 1, 0}), C3({0, 0, 1}), C3({1, 1, 0}),
C3({1, 0, 1}), C3({1, 0, 1}), C3({1, 1, 1})};
const auto realmags = {1e-19, 1e-17, 1e-15, 1e-12, 1e-9, 1e-6, 1e-3, .03, 1., 1e1, 1e2, 1e3, 1e4};
const auto sidemags = {-1e-15, +1e-14, -1e-11, 1e-7, -1e-3, .1, 1., 1.4142, 1.7321};
const auto imagrels = {0., -4e-16, +8e-16, -5e-11, 3e-7, -2e-3, .01, .1};
int run_test_for_many_q(std::function<complex_t(C3)> fff0, std::function<complex_t(C3)> fff1,
double qmag_min, double qmag_max, double eps, bool real_only)
{
auto evaluate = [&](C3 q) -> std::tuple<complex_t, complex_t, double, double, double> {
const complex_t f0 = fff0(q);
const complex_t f1 = fff1(q);
const double avge = (std::abs(f0) + std::abs(f1)) / 2;
const double abserr = real_only
? std::abs(f0) - std::abs(f1)
: std::max(fabs(real(f0) - real(f1)), fabs(imag(f0) - imag(f1)));
const double deviation = abserr / avge * std::min(1., eps * avge / 1e-16);
return {f0, f1, avge, abserr, deviation};
};
double max_deviation = 0;
C3 q_at_max;
int failures = 0;
for (const C3& q_maindir : maindirs) {
for (const C3& q_sidedir : sidedirs) {
for (const double qrealmag : realmags) {
for (const double qsidemag : sidemags) {
for (const double qimagrel : imagrels) {
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_or_throw();
const auto [f0, f1, avge, abserr, deviation] = evaluate(q);
if (deviation > eps) {
++failures;
if (deviation > max_deviation) {
max_deviation = deviation;
q_at_max = q;
}
}
}
}
}
}
}
if (failures) {
const auto [f0, f1, avge, abserr, deviation] = evaluate(q_at_max);
std::cout << "Deviation at q = " << q_at_max << ":\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 = " << deviation
<< ", limit = " << eps << "\n";
}
return failures;
}
} // namespace formfactorTest
|
