GCC Code Coverage Report

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