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 |