GCC Code Coverage Report


Directory: ./
File: src/analyzers.cpp
Date: 2025-02-03 10:58:24
Exec Total Coverage
Lines: 0 85 0.0%
Functions: 0 7 0.0%
Branches: 0 33 0.0%

Line Branch Exec Source
1 #include "analyzers.hpp"
2 #include "interaction.hpp"
3 #include "particle.hpp"
4 #include "potentials.hpp"
5 #include "types.hpp"
6 #include <algorithm>
7 #include <iterator>
8 #include <map>
9 #include <numeric>
10 #include <iostream>
11
12 Clusterer::Clusterer(std::vector<Particle> &particles, const Vec &L, double cutoff)
13 : particles(particles), L(L), cutoff(cutoff), cutoffsq(cutoff * cutoff) {
14 }
15
16 void Clusterer::reset() {
17 int n = particles.size();
18 particle2Coordination.resize(n);
19 std::fill(particle2Coordination.begin(), particle2Coordination.end(), 0);
20 particle2Cluster.resize(n);
21 cluster2Particles.clear();
22 for (int i = 0; i < n; ++i) {
23 particle2Cluster[i] = i;
24 cluster2Particles[i] = {i};
25 }
26 }
27
28 double Clusterer::operator()(Particle &p1, Particle &p2) {
29 int i = std::distance(&particles.front(), &p1);
30 int j = std::distance(&particles.front(), &p2);
31 double distsq = r_ij(L, particles[i].r, particles[j].r).squaredNorm();
32 if (distsq < cutoffsq && particle2Cluster[i] != particle2Cluster[j]) {
33 particle2Coordination[i]++;
34 particle2Coordination[j]++;
35 int oldCluster = std::max(particle2Cluster[i], particle2Cluster[j]);
36 int newCluster = std::min(particle2Cluster[i], particle2Cluster[j]);
37 for (int i : cluster2Particles[oldCluster]) {
38 particle2Cluster[i] = newCluster;
39 }
40 cluster2Particles[newCluster].splice(cluster2Particles[newCluster].end(),
41 cluster2Particles[oldCluster]);
42 }
43 return 0.;
44 }
45
46 void Clusterer::sortClusters() {
47 // At this point, particles of the same cluster have the same id in particle2Cluster, which can
48 // still be in the range [0, total number of particles)
49 // Now change data structures such that cluster ids go from 0 to numClusters - 1 and are ordered
50 // by cluster size
51 int n = particles.size();
52 std::set<int> clusterIdSet(particle2Cluster.begin(),
53 particle2Cluster.end()); // get all unique ids
54 std::vector<int> clusterIdVector(clusterIdSet.begin(),
55 clusterIdSet.end()); // transform to vector for accessible index
56 int numClusters = clusterIdVector.size();
57 // Construct map oldId -> newId where oldId in [0, numParticles) and newId in [0, numClusters)
58 std::map<int, int> clusterIdMap;
59 for (int newId = 0; newId < numClusters; ++newId) {
60 int oldId = clusterIdVector[newId];
61 clusterIdMap[oldId] = newId;
62 }
63 // Change cluster ids of data structures to new restricted range
64 for (int &clusterId : particle2Cluster) {
65 clusterId = clusterIdMap[clusterId];
66 }
67 std::map<int, std::list<int>> cluster2ParticlesRestricted;
68 for (int i = 0; i < n; ++i) {
69 int newId = clusterIdMap[i];
70 cluster2ParticlesRestricted[newId].splice(cluster2ParticlesRestricted[newId].end(),
71 cluster2Particles[i]);
72 }
73 // At this point, cluster ids are in [0, numClusters), but they are not yet sorted by the cluster
74 // size
75 std::vector<int> clusterSizes(numClusters, 0);
76 for (int clusterId : particle2Cluster) {
77 clusterSizes[clusterId]++;
78 }
79 std::vector<int> permutation(numClusters);
80 std::iota(permutation.begin(), permutation.end(), 0);
81 std::sort(permutation.begin(), permutation.end(),
82 [&clusterSizes](int a, int b) { return clusterSizes[a] > clusterSizes[b]; });
83 std::vector<int> inversePermutation(numClusters);
84 for (int i = 0; i < numClusters; ++i) {
85 inversePermutation[permutation[i]] = i;
86 }
87 for (int &clusterId : particle2Cluster) {
88 clusterId = inversePermutation[clusterId];
89 }
90 cluster2Particles.clear();
91 for (int i = 0; i < numClusters; ++i) {
92 cluster2Particles[i].splice(cluster2Particles[i].end(), cluster2ParticlesRestricted[permutation[i]]);
93 }
94 }
95
96 PercolationAnalyzer::PercolationAnalyzer(System *system)
97 : system(system), tLast(std::numeric_limits<double>::quiet_NaN()),
98 clusterer(system->particles, system->L, system->interaction->internal->range) {}
99
100 void PercolationAnalyzer::update() {
101 if (system->t == tLast) {
102 return;
103 }
104 this->tLast = system->t;
105 int n = system->particles.size();
106 // Do the clustering algorithm
107 clusterer.reset();
108 system->interaction->iteratePairs(clusterer, true, false);
109 clusterer.sortClusters();
110 // Do the coordination statistics
111 // Just go through all particles and make a histogram of the coordinations
112 coordinationCount.clear();
113 coordinationCount.resize(n + 1);
114 clusterSizeCount.clear();
115 clusterSizeCount.resize(n + 1);
116 for (int i = 0; i < n; ++i) {
117 coordinationCount[i] = 0;
118 clusterSizeCount[i] = 0;
119 }
120 for (int coord : clusterer.particle2Coordination) {
121 coordinationCount[coord]++;
122 }
123 numClusters = clusterer.cluster2Particles.size();
124 clusterSizeMax = numClusters == 0 ? 0 : clusterer.cluster2Particles[0].size();
125 for (int i = 0; i < numClusters; ++i) {
126 clusterSizeCount[i] = clusterer.cluster2Particles[i].size();
127 }
128 }
129
130 Analyzer *AnalyzerFactory::create(System *system, const std::string &key) {
131 if (key == "percolation") {
132 return new PercolationAnalyzer(system);
133 } else {
134 throw std::invalid_argument("No " + key + " analyzer found");
135 }
136 }
137