GCC Code Coverage Report

Directory: ./
File: include/observables.hpp
Date: 2025-11-05 16:49:12
Exec Total Coverage
Lines: 36 154 23.4%
Functions: 15 55 27.3%
Branches: 27 58 46.6%
Line Branch Exec Source
1 #pragma once
2
3 #include <H5Cpp.h>
4 #include <iostream>
5 #include <stdexcept>
6 #include <string>
7 #include <unordered_map>
8 #include <vector>
9 #include <yaml-cpp/yaml.h>
10
11 #include "multiindex.hpp"
12 #include "particle.hpp"
13 #include "analyzers.hpp"
14 #include "rng.hpp"
15 #include "system.hpp"
16 #include "types.hpp"
17
18 template <int d> class Histogram;
19
20 class Observables {
21 protected:
22 virtual void doSample(double weight) = 0;
23
24 public:
25 System *system;
26 int numSpeciesIndices;
27 int numSpecies;
28 std::vector<std::string> keys;
29 double totalWeight;
30 double tLast;
31 double every;
32
33 Observables(System *system, int numSpeciesIndices, const YAML::Node &config);
34 ✗ virtual ~Observables() = default;
35
36 virtual std::string type() = 0;
37 virtual std::vector<std::string> axisLabels() = 0;
38 virtual std::vector<int> shape() = 0;
39 virtual std::vector<double> coordinates(int axis) = 0;
40 virtual std::vector<double> &data(const std::string &key,
41 Eigen::Array<int, Eigen::Dynamic, 1> species) = 0;
42
43 ✗ virtual std::vector<std::string> analyzers() { return {}; };
44 void sample(double weight);
45 };
46
47 template <int d, int speciesIndices> class Fields : public Observables {
48 protected:
49 6 void customExtent(const YAML::Node &config, const std::string &key, int axis) {
50
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 6 if (config[key]) {
51 ✗ double eNew = YAML::parse<double>(config, key);
52 ✗ if (eNew < extent[axis]) {
53 ✗ extent[axis] = eNew;
54 ✗ dr = extent / bins.template cast<double>();
55 } else {
56 std::clog
57 << "WARNING: You have set the custom extent " << key << " to " << std::to_string(eNew)
58 ✗ << " in Observable " << this->type()
59 ✗ << ", which is above the maximum allowed value. This custom setting will be ignored";
60 }
61 }
62 }
63
64 virtual double dV(const Eigen::Array<int, d, 1> &index) = 0;
65
66 public:
67 Eigen::Array<int, d, 1> bins;
68 Eigen::Array<double, d, 1> extent;
69 Eigen::Array<double, d, 1> dr;
70 Multiindex<speciesIndices> speciesMultiindex;
71 std::unordered_map<std::string, std::vector<Histogram<d>>> histograms;
72
73 18 Fields(System *system, const Eigen::Array<int, d, 1> &bins,
74 const Eigen::Array<double, d, 1> &extent, const YAML::Node &config)
75 18 : Observables(system, speciesIndices, config), bins(bins), extent(extent),
76 18 dr(extent / bins.template cast<double>()),
77
2/3
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
 18 speciesMultiindex(Eigen::Array<int, speciesIndices, 1>::Constant(numSpecies)) {
78
7/12
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 18 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 18 times.
✓ Branch 15 taken 9 times.
✓ Branch 16 taken 9 times.
✗ Branch 17 not taken.
 126 for (auto &k : config["keys"]) {
79
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 std::string key = k.as<std::string>();
80
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 36 this->keys.push_back(key);
81
3/6
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
 72 this->histograms[key] = std::vector<Histogram<d>>(speciesMultiindex.linsize, Histogram<d>(bins));
82 }
83 }
84
85 18 virtual ~Fields() = default;
86
87 24 std::vector<int> shape() {
88 24 std::vector<int> shape;
89
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 11 times.
 78 for (int axis = 0; axis < d; axis++) {
90
1/2
✓ Branch 2 taken 27 times.
✗ Branch 3 not taken.
 54 shape.push_back(bins[axis]);
91 }
92 24 return shape;
93 }
94
95 8 std::vector<double> coordinates(int axis) {
96 8 std::vector<double> coords;
97 8 Eigen::Array<int, d, 1> index = Eigen::Array<int, d, 1>::Zero();
98 8 Eigen::Array<double, d, 1> real;
99
2/2
✓ Branch 3 taken 52 times.
✓ Branch 4 taken 4 times.
 112 for (; index[axis] < bins[axis]; ++index[axis]) {
100 104 real = index2real(index);
101
1/2
✓ Branch 2 taken 52 times.
✗ Branch 3 not taken.
 104 coords.push_back(real[axis]);
102 }
103
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 coords.push_back(extent[axis]);
104 8 return coords;
105 }
106
107
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 6 std::vector<double> &data(
108 const std::string &key,
109 Eigen::Array<int, Eigen::Dynamic, 1> species = Eigen::Array<int, speciesIndices, 1>::Zero()) {
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 6 if (species.size() != speciesIndices) {
111 ✗ throw std::invalid_argument("Wrong number of species given in data()");
112 }
113 6 Histogram<d> &histogram = histograms.at(key).at(speciesMultiindex.multi2lin(species));
114
2/2
✓ Branch 0 taken 746 times.
✓ Branch 1 taken 3 times.
 1498 for (int i = 0; i < histogram.multiindex.linsize; ++i) {
115 1494 double normFactor = dV(histogram.multiindex.lin2multi(i)) * totalWeight;
116 1492 histogram.valuesNormed[i] = histogram.values[i] / normFactor;
117 }
118 6 return histogram.valuesNormed;
119 }
120
121 11 Eigen::Array<int, d, 1> real2index(const Eigen::Array<double, d, 1> &real) {
122
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 11 return (real / dr).template cast<int>();
123 }
124
125 78 Eigen::Array<double, d, 1> index2real(const Eigen::Array<int, d, 1> &index) {
126 78 return index.template cast<double>() * dr;
127 }
128 };
129
130 class Scalar : public Fields<0, 0> {
131 public:
132 inline static std::unordered_map<std::string, double (*)(System *s)> functions = {
133 {"N", [](System *s) -> double { return s->particles.size(); }},
134 {"E", [](System *s) -> double { return s->E; }},
135 {"clusterSizeMax",
136 ✗ [](System *s) -> double {
137 ✗ PercolationAnalyzer *percolation =
138 ✗ static_cast<PercolationAnalyzer *>(s->analyzers["percolation"]);
139 ✗ return percolation->clusterSizeMax;
140 }},
141 };
142
143 double dV(const Eigen::Array<int, 0, 1> &index);
144 void doSample(double weight);
145
146 Scalar(System *system, const YAML::Node &config);
147
148 ✗ std::string type() { return "scalar"; }
149 ✗ std::vector<std::string> axisLabels() { return {}; }
150 std::vector<std::string> analyzers();
151 };
152
153 class Onebody : public Fields<DIM, 1> {
154 double dV_;
155 double dV(const Eigen::Array<int, DIM, 1> &index);
156 void doSample(double weight);
157
158 public:
159 inline static std::unordered_map<std::string, double (*)(System *s, const Particle &p)>
160 functions = {
161 {"1", [](System *s, const Particle &p) -> double { return 1.; }},
162 {"Eext", [](System *s, const Particle &p) -> double { return p.Eext; }},
163 {"Fextx", [](System *s, const Particle &p) -> double { return p.Fext[0]; }},
164 {"Fexty", [](System *s, const Particle &p) -> double { return p.Fext[1]; }},
165 {"Fextz", [](System *s, const Particle &p) -> double { return p.Fext[2]; }},
166 {"Fintx", [](System *s, const Particle &p) -> double { return p.Fint[0]; }},
167 {"Finty", [](System *s, const Particle &p) -> double { return p.Fint[1]; }},
168 {"Fintz", [](System *s, const Particle &p) -> double { return p.Fint[2]; }},
169 {"rRandx", [](System *s, const Particle &p) -> double { return p.rRand[0]; }},
170 {"rRandy", [](System *s, const Particle &p) -> double { return p.rRand[1]; }},
171 {"rRandz", [](System *s, const Particle &p) -> double { return p.rRand[2]; }},
172 {"N", [](System *s, const Particle &p) -> double { return s->particles.size(); }},
173 {"H", [](System *s, const Particle &p) -> double { return s->E; }},
174 {"clusterSizeMax",
175 ✗ [](System *s, const Particle &p) -> double {
176 ✗ PercolationAnalyzer *percolation = static_cast<PercolationAnalyzer *>(s->analyzers["percolation"]);
177 ✗ return percolation->clusterSizeMax;
178 }
179 },
180 };
181
182 Onebody(System *system, const YAML::Node &config);
183
184 ✗ std::string type() { return "onebody"; }
185 ✗ std::vector<std::string> axisLabels() { return {"x", "y", "z"}; }
186 std::vector<std::string> analyzers();
187 };
188
189 ✗ static Vec deltaF_deltar(System *s, Particle &p, const Vec &e_r, const double deltar = 1e-5) {
190 // Save current state
191 ✗ Vec r_before = p.r;
192 ✗ Vec F_before = p.Fint;
193 // Shift particle for finite difference
194 ✗ p.r += deltar * e_r;
195 ✗ s->applyPeriodic(p);
196 ✗ s->updateParticleState(p);
197 ✗ Vec F_after = p.Fint;
198 // Restore initial state
199 ✗ p.r = r_before;
200 ✗ s->updateParticleMeta(p);
201 ✗ p.Fint = F_before;
202 ✗ return (F_after - F_before) / deltar;
203 }
204
205 ✗ static Vec deltaF1_deltar0(System *s, Particle &p0, Particle &p1, const Vec &e_r0,
206 double deltar0 = 1e-5) {
207 // Save current state
208 ✗ Vec r0_before = p0.r;
209 ✗ Vec F1_before = p1.Fint;
210 // Shift particle 0 for finite difference
211 ✗ deltar0 *= rand() % 2 ? 1 : -1;
212 ✗ p0.r += deltar0 * e_r0;
213 ✗ s->applyPeriodic(p0);
214 ✗ s->updateParticleMeta(p0);
215 ✗ s->updateParticleState(p1);
216 ✗ Vec F1_after = p1.Fint;
217 // Restore initial state
218 ✗ p0.r = r0_before;
219 ✗ s->updateParticleMeta(p0);
220 ✗ p1.Fint = F1_before;
221 ✗ return (F1_after - F1_before) / deltar0;
222 }
223
224 ✗ static Vec deltaF1_deltar0_fullUpdate(System *s, Particle &p0, Particle &p1, const Vec &e_r0,
225 const double deltar0 = 1e-5) {
226 // Save current state
227 ✗ static std::vector<Particle> particlesBefore = s->particles;
228 ✗ Vec F1_before = p1.Fint;
229 // Shift particle 0 for finite difference
230 ✗ p0.r += deltar0 * e_r0;
231 ✗ s->applyPeriodic(p0);
232 ✗ s->updateState();
233 ✗ Vec F1_after = p1.Fint;
234 // Restore initial state
235 ✗ s->particles = particlesBefore;
236 ✗ s->updateParticleMeta(p0);
237 ✗ return (F1_after - F1_before) / deltar0;
238 }
239
240 ✗ inline static Vec e_tan_safe(const Vec &e_r, const Vec &F1, const Vec &F2 = Vec::Zero(),
241 const double tol = 1e-4) {
242 ✗ Vec e_tan = (F1 - F1.dot(e_r) * e_r).normalized();
243 ✗ if (std::abs(e_tan.norm() - 1.) > tol || std::abs(e_tan.dot(e_r)) > tol) {
244 ✗ e_tan = (F2 - F2.dot(e_r) * e_r).normalized();
245 ✗ if (std::abs(e_tan.norm() - 1.) > tol || std::abs(e_tan.dot(e_r)) > tol) {
246 ✗ if constexpr (DIM == 3) {
247 ✗ e_tan = e_r.cross(Vec{1, 0, 0}).normalized();
248 }
249 }
250 }
251 ✗ return e_tan;
252 }
253
254 ✗ inline static Vec e_tan_random(const Vec &e_r, const double tol = 1e-4) {
255 ✗ static RNG rng{};
256 Vec e_tan;
257 ✗ do {
258 ✗ Vec randVec = rng.vecNormal().normalized();
259 ✗ e_tan = (randVec - randVec.dot(e_r) * e_r).normalized();
260 ✗ } while (std::abs(e_tan.norm() - 1.) > tol || std::abs(e_tan.dot(e_r)) > tol);
261 ✗ return e_tan;
262 }
263
264 class Radial : public Fields<1, 2> {
265 double dV(const Eigen::Array<int, 1, 1> &index);
266 void doSample(double weight);
267
268 public:
269 inline static std::unordered_map<std::string, double (*)(System *s, Particle &p0, Particle &p1,
270 const Vec &e_r)>
271 functions = {
272 {"1", [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double { return 1.; }},
273 {"F0_rad",
274 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
275 ✗ return p0.F().dot(e_r);
276 }},
277 {"F1_rad",
278 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
279 ✗ return p1.F().dot(e_r);
280 }},
281 {"F0_rad_F1_rad",
282 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
283 ✗ return p0.F().dot(e_r) * p1.F().dot(e_r);
284 }},
285 {"F0_tan_dot_F1_tan",
286 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
287 ✗ Vec F0 = p0.F();
288 ✗ Vec F1 = p1.F();
289 ✗ Vec F0tan = F0 - F0.dot(e_r) * e_r;
290 ✗ Vec F1tan = F1 - F1.dot(e_r) * e_r;
291 ✗ return F0tan.dot(F1tan);
292 }},
293 {"(F0_tan_cross_F1_tan)_rad",
294 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
295 ✗ if constexpr (DIM == 3) {
296 ✗ Vec F0 = p0.F();
297 ✗ Vec F1 = p1.F();
298 ✗ Vec F0tan = F0 - F0.dot(e_r) * e_r;
299 ✗ Vec F1tan = F1 - F1.dot(e_r) * e_r;
300 ✗ return F0tan.cross(F1tan).dot(e_r);
301 } else {
302 return 0.;
303 }
304 }},
305 {"F0_dot_F1",
306 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
307 ✗ return p0.F().dot(p1.F());
308 }},
309 {"grad0(F0)_rad",
310 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
311 ✗ return deltaF_deltar(s, p0, e_r).dot(e_r);
312 }},
313 {"grad0(F0)_tan",
314 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
315 ✗ Vec e_tan = e_tan_safe(e_r, p0.F());
316 ✗ return deltaF_deltar(s, p0, e_tan).dot(e_tan);
317 }},
318 {"grad0(F1)_rad",
319 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
320 ✗ return deltaF1_deltar0(s, p0, p1, e_r).dot(e_r);
321 }},
322 {"grad0(F1)_tan",
323 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
324 ✗ Vec e_tan = e_tan_safe(e_r, p0.F(), p1.F());
325 ✗ return deltaF1_deltar0(s, p0, p1, e_tan).dot(e_tan);
326 }},
327 {"grad0(F1)_tan_u",
328 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
329 ✗ Vec e_tan = e_tan_safe(e_r, p0.u, p1.u);
330 ✗ return deltaF1_deltar0(s, p0, p1, e_tan).dot(e_tan);
331 }},
332 {"grad0(F1)_tan_random",
333 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
334 ✗ Vec e_tan = e_tan_random(e_r);
335 ✗ return deltaF1_deltar0(s, p0, p1, e_tan).dot(e_tan);
336 }},
337 {"grad0(F1)_rad_fullUpdate",
338 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
339 ✗ return deltaF1_deltar0_fullUpdate(s, p0, p1, e_r).dot(e_r);
340 }},
341 {"grad0(F1)_tan_fullUpdate",
342 ✗ [](System *s, Particle &p0, Particle &p1, const Vec &e_r) -> double {
343 ✗ Vec e_tan = e_tan_safe(e_r, p0.F(), p1.F());
344 ✗ return deltaF1_deltar0_fullUpdate(s, p0, p1, e_tan).dot(e_tan);
345 }},
346 };
347
348 Radial(System *system, const YAML::Node &config);
349
350 ✗ std::string type() { return "radial"; }
351 ✗ std::vector<std::string> axisLabels() { return {"r"}; }
352 };
353
354 class TwobodyPlanar : public Fields<3, 2> {
355 double dV(const Eigen::Array<int, 3, 1> &index);
356 void doSample(double weight);
357
358 public:
359 inline static std::unordered_map<std::string,
360 double (*)(System *s, const Particle &p0, const Particle &p1)>
361 functions = {
362 {"1", [](System *s, const Particle &p0, const Particle &p1) -> double { return 1.; }},
363 };
364
365 TwobodyPlanar(System *system, const YAML::Node &config);
366
367 ✗ std::string type() { return "twobodyplanar"; }
368 ✗ std::vector<std::string> axisLabels() { return {"z_1", "z_2", "r^*"}; }
369 };
370
371 class Clusters : public Observables {
372 void doSample(double weight);
373
374 public:
375 inline static std::unordered_map<std::string,
376 double (*)(PercolationAnalyzer *percolation, int i, int n)>
377 functions = {
378 {"clusters",
379 [](PercolationAnalyzer *percolation, int i, int n) -> double {
380 return (double)percolation->clusterSizeCount[i] / n;
381 }},
382 {"coordinations",
383 [](PercolationAnalyzer *percolation, int i, int n) -> double {
384 return (double)percolation->coordinationCount[i] / n;
385 }},
386 {"P0",
387 [](PercolationAnalyzer *percolation, int i, int n) -> double {
388 return percolation->coordinationCount[0] == i ? 1. : 0.;
389 }},
390 {"P1",
391 [](PercolationAnalyzer *percolation, int i, int n) -> double {
392 return percolation->coordinationCount[1] == i ? 1. : 0.;
393 }},
394 {"P2",
395 [](PercolationAnalyzer *percolation, int i, int n) -> double {
396 return percolation->coordinationCount[2] == i ? 1. : 0.;
397 }},
398 {"P3",
399 [](PercolationAnalyzer *percolation, int i, int n) -> double {
400 return percolation->coordinationCount[3] == i ? 1. : 0.;
401 }},
402 {"P4",
403 [](PercolationAnalyzer *percolation, int i, int n) -> double {
404 return percolation->coordinationCount[4] == i ? 1. : 0.;
405 }},
406 };
407
408 int nMax;
409 std::unordered_map<std::string, std::vector<double>> distributions;
410 std::unordered_map<std::string, std::vector<double>> distributionsNormed;
411
412 Clusters(System *system, const YAML::Node &config);
413
414 ✗ std::string type() { return "clusters"; }
415 ✗ std::vector<std::string> axisLabels() { return {"n"}; }
416 ✗ std::vector<std::string> analyzers() { return {"percolation"}; };
417 std::vector<int> shape();
418 std::vector<double> coordinates(int axis);
419 std::vector<double> &data(const std::string &key = "clusters",
420 Eigen::Array<int, Eigen::Dynamic, 1> species = {});
421 };
422
423 class ObservablesFactory {
424 public:
425 static Observables *create(System *system, const std::string &category, const YAML::Node &config);
426 };
427