GCC Code Coverage Report

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