GCC Code Coverage Report

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