GCC Code Coverage Report

Directory: ./
File: src/observables.cpp
Date: 2025-11-05 16:49:12
Exec Total Coverage
Lines: 78 164 47.6%
Functions: 13 25 52.0%
Branches: 57 115 49.6%
Line Branch Exec Source
1 #include <H5Cpp.h>
2 #include <numeric>
3 #include <string>
4 #include <vector>
5
6 #include "histogram.hpp"
7 #include "observables.hpp"
8 #include "particle.hpp"
9 #include "analyzers.hpp"
10 #include "system.hpp"
11
12 12 template <class T> static void checkFunctionsExist(T *t, bool withUnweighted = false) {
13
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 6 times.
 42 for (std::string &key : t->keys) {
14 30 std::string functionKey = key;
15
7/10
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 3 times.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 9 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 15 times.
 48 if (withUnweighted && key.length() > 3 && key.substr(key.length() - 3) == "_dt") {
16 ✗ functionKey = key.substr(0, key.length() - 3);
17 }
18
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 15 times.
 30 if (t->functions.find(functionKey) == t->functions.end()) {
19 ✗ throw std::invalid_argument("There is no observable " + key + " for type " + t->type());
20 }
21 }
22 }
23
24 9 Observables::Observables(System *system, int numSpeciesIndices, const YAML::Node &config)
25 18 : system(system), numSpeciesIndices(numSpeciesIndices), numSpecies(system->getNumSpecies()),
26
2/4
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
 18 totalWeight(0.), tLast(-INFINITY), every(YAML::parse<double>(config, "every", 0.)) {}
27
28 3 void Observables::sample(double weight) {
29
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if (system->t > tLast + every) {
30 3 doSample(weight);
31 3 totalWeight += weight;
32 3 tLast = system->t;
33 }
34 3 }
35
36 3 Scalar::Scalar(System *system, const YAML::Node &config) : Fields(system, {}, {}, config) {}
37
38 ✗ std::vector<std::string> Scalar::analyzers() {
39 ✗ for (const std::string &key : keys) {
40 ✗ if (key == "clusterSizeMax") {
41 ✗ return {"percolation"};
42 }
43 }
44 ✗ return {};
45 }
46
47 1 double Scalar::dV(const Eigen::Array<int, 0, 1> &index) { return 1.; }
48
49 1 void Scalar::doSample(double weight) {
50
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for (const std::string &key : keys) {
51 1 const auto &function = functions.at(key);
52 1 auto &histogram = histograms.at(key)[0];
53 1 histogram[0] += function(system) * weight;
54 }
55 1 }
56
57 3 Onebody::Onebody(System *system, const YAML::Node &config)
58
4/8
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
 6 : Fields(system, YAML::parse<ArrayI>(config, "bins"), system->L, config), dV_(dr.prod()) {
59
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 checkFunctionsExist(this, true);
60 3 }
61
62 ✗ std::vector<std::string> Onebody::analyzers() {
63 ✗ for (const std::string &key : keys) {
64 ✗ if (key == "clusterSizeMax") {
65 ✗ return {"percolation"};
66 }
67 }
68 ✗ return {};
69 }
70
71 720 double Onebody::dV(const ArrayI &index) { return dV_; }
72
73 1 void Onebody::doSample(double weight) {
74 1 const auto &particles = system->particles;
75 1 Eigen::Array<int, 1, 1> species;
76
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (const std::string &key : keys) {
77 4 double maybeWeight = weight;
78
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::string functionKey = key;
79
5/8
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4 times.
 7 if (key.length() > 3 && key.substr(key.length() - 3) == "_dt") {
80 ✗ maybeWeight = 1.;
81 ✗ functionKey = key.substr(0, key.length() - 3);
82 }
83
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto &function = functions.at(functionKey);
84
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 auto &hKey = histograms.at(key);
85
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 4 times.
 12 for (const Particle &p : particles) {
86
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 species[0] = p.species;
87 8 hKey[speciesMultiindex.multi2lin(species)][real2index(p.r)] +=
88
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 function(system, p) * maybeWeight;
89 }
90 4 }
91 1 }
92
93 3 Radial::Radial(System *system, const YAML::Node &config)
94
2/4
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
 6 : Fields(system, YAML::parse<Eigen::Array<int, 1, 1>>(config, "bins"),
95
1/2
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 6 Eigen::Array<double, 1, 1>{0.5 * system->L.minCoeff()}, config) {
96
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 checkFunctionsExist(this);
97
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
 6 customExtent(config, "rMax", 0);
98 3 }
99
100 25 double Radial::dV(const Eigen::Array<int, 1, 1> &index) {
101 // Calculate the volume of the spherical shell from r to r + dr;
102 25 double r = index2real(index)[0];
103 25 double R = r + dr[0];
104 25 return 4. / 3. * M_PI * (R * R * R - r * r * r);
105 }
106
107 1 void Radial::doSample(double weight) {
108 1 Eigen::Array<double, 1, 1> real;
109 1 Eigen::Array<int, 2, 1> species;
110 1 Vec rVec;
111 1 Vec e_r;
112 1 double V = system->L.prod();
113 1 auto &particles = system->particles;
114 1 std::vector<size_t> Ns(system->getNumSpecies(), 0);
115
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (const Particle &p : particles) {
116 2 ++Ns[p.species];
117 }
118
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for (const std::string &key : keys) {
119
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto &histogram = histograms.at(key);
120
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto &function = functions.at(key);
121
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 for (Particle &p0 : particles) {
122
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 for (Particle &p1 : particles) {
123
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 if (&p0 == &p1) {
124 2 continue;
125 }
126
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 rVec = r_ij(system->L, p0.r, p1.r);
127
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 real[0] = rVec.norm();
128 2 e_r = rVec / real[0];
129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (real[0] >= extent[0]) {
130 ✗ continue;
131 }
132
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 species[0] = p0.species;
133 2 species[1] = p1.species;
134 2 histogram[speciesMultiindex.multi2lin(species)][real2index(real)] +=
135
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 function(system, p0, p1, e_r) * V / (Ns[p0.species] * Ns[p1.species]) * weight;
136 }
137 }
138 }
139 1 }
140
141 ✗ TwobodyPlanar::TwobodyPlanar(System *system, const YAML::Node &config)
142 ✗ : Fields(system, YAML::parse<Eigen::Array<int, 3, 1>>(config, "bins"),
143 ✗ Eigen::Array<double, 3, 1>{system->L[DIM - 1], system->L[DIM - 1],
144 ✗ 0.5 * system->L.head<DIM - 1>().minCoeff()},
145 ✗ config) {
146 ✗ checkFunctionsExist(this);
147 ✗ customExtent(config, "rMax", 2);
148 }
149
150 ✗ double TwobodyPlanar::dV(const Eigen::Array<int, 3, 1> &index) {
151 ✗ double r = index2real(index)[2];
152 ✗ double R = r + dr[2];
153 ✗ return dr[0] * dr[1] * M_PI * (R * R - r * r);
154 }
155
156 ✗ void TwobodyPlanar::doSample(double weight) {
157 ✗ Eigen::Array<double, 3, 1> real;
158 ✗ Eigen::Array<int, 2, 1> species;
159 ✗ Vec rij;
160 ✗ const auto &particles = system->particles;
161 ✗ for (const std::string &key : keys) {
162 ✗ auto &histogram = histograms.at(key);
163 ✗ const auto &function = functions.at(key);
164 ✗ for (auto &p0 : particles) {
165 ✗ for (auto &p1 : particles) {
166 ✗ if (&p0 == &p1) {
167 ✗ continue;
168 }
169 ✗ rij = r_ij(system->L, p0.r, p1.r);
170 ✗ rij[DIM - 1] = 0.;
171 ✗ real = {p0.r[DIM - 1], p1.r[DIM - 1], rij.norm()};
172 ✗ if (real[2] >= extent[2]) {
173 ✗ continue;
174 }
175 ✗ species[0] = p0.species;
176 ✗ species[1] = p0.species;
177 ✗ histogram[speciesMultiindex.multi2lin(species)][real2index(real)] +=
178 ✗ function(system, p0, p1) * weight;
179 }
180 }
181 }
182 }
183
184 ✗ Clusters::Clusters(System *system, const YAML::Node &config) : Observables(system, 0, config) {
185 ✗ this->keys = {"clusters", "coordinations", "P0", "P1", "P2", "P3", "P4"};
186 ✗ for (const std::string &key : this->keys) {
187 ✗ distributions[key] = std::vector<double>(system->particles.size());
188 ✗ distributionsNormed[key] = std::vector<double>(system->particles.size());
189 }
190 }
191
192 ✗ void Clusters::doSample(double weight) {
193 ✗ size_t n = system->particles.size();
194 ✗ PercolationAnalyzer *percolation = static_cast<PercolationAnalyzer *>(system->analyzers["percolation"]);
195 ✗ for (const std::string &key : keys) {
196 ✗ auto &dist = distributions.at(key);
197 ✗ const auto &function = functions.at(key);
198 ✗ nMax = std::max(dist.size() - 1, n);
199 ✗ dist.resize(nMax + 1, 0);
200 ✗ for (size_t i = 0; i <= n; ++i) {
201 ✗ dist[i] += function(percolation, i, n) * weight;
202 }
203 }
204 }
205
206 ✗ std::vector<int> Clusters::shape() { return {nMax + 1}; }
207
208 ✗ std::vector<double> Clusters::coordinates(int axis) {
209 ✗ std::vector<double> result(nMax + 1);
210 ✗ std::iota(result.begin(), result.end(), 0);
211 ✗ return result;
212 }
213
214 ✗ std::vector<double> &Clusters::data(const std::string &key,
215 Eigen::Array<int, Eigen::Dynamic, 1> species) {
216 ✗ auto &dist = distributions.at(key);
217 ✗ auto &distNormed = distributionsNormed.at(key);
218 ✗ distNormed.resize(dist.size());
219 ✗ for (size_t i = 0; i < dist.size(); ++i) {
220 ✗ distNormed[i] = dist[i] / totalWeight;
221 }
222 ✗ return distNormed;
223 }
224
225 ✗ Observables *ObservablesFactory::create(System *system, const std::string &category,
226 const YAML::Node &config) {
227 ✗ if (category == "Scalar") {
228 ✗ return new Scalar(system, config);
229 ✗ } else if (category == "Onebody") {
230 ✗ return new Onebody(system, config);
231 ✗ } else if (category == "Radial") {
232 ✗ return new Radial(system, config);
233 ✗ } else if (category == "TwobodyPlanar") {
234 ✗ return new TwobodyPlanar(system, config);
235 ✗ } else if (category == "Clusters") {
236 ✗ return new Clusters(system, config);
237 } else {
238 ✗ throw std::invalid_argument("No observable category of type " + category + " found");
239 }
240 }
241