franzblau.cpp
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------------
2 // d-SEAMS is free software: you can redistribute it and/or modify
3 // it under the terms of the GNU General Public License as published by
4 // the Free Software Foundation, either version 3 of the License, or
5 // (at your option) any later version.
6 //
7 // A copy of the GNU General Public License is available at
8 // http://www.gnu.org/licenses/
9 //-----------------------------------------------------------------------------------
10 
11 #include <franzblau.hpp>
12 
29  //
30  primitive::Graph fullGraph; // Graph object, contains the connectivity
31  // information from the neighbourlist
32 
33  // Find all possible rings, using backtracking. This may contain non-primitive
34  // rings as well.
35  fullGraph = primitive::countAllRingsFromIndex(nList, maxDepth);
36 
37  // Remove all non-SP rings using the Franzblau algorithm.
38  fullGraph = primitive::removeNonSPrings(&fullGraph);
39 
40  // The rings vector of vectors inside the fullGraph graph object is the ring
41  // network information we want
42  return fullGraph.rings;
43 }
44 
61  int maxDepth) {
62  //
63  primitive::Graph fullGraph; // Graph object
65  visited; // Contains the indices (NOT IDs) visited for a particular node
66  int depth; // Current depth
67 
68  // ------------------------------
69  // Init
70  // Initialize the graph object with all the information from the neighbour
71  // list (of indices only)
72  fullGraph = primitive::populateGraphFromIndices(neighHbondList);
73  // ------------------------------
74  // Loop through every vertex
75  for (int iatom = 0; iatom < neighHbondList.size(); iatom++) {
76  visited.clear();
77  depth = 0;
78  // Call the function for getting rings
79  primitive::findRings(&fullGraph, iatom, &visited, maxDepth, depth);
80  } // loop through every vertex
81  // ------------------------------
82  // Restore back-up of the edges (some may have been removed)
83  fullGraph = primitive::restoreEdgesFromIndices(&fullGraph, neighHbondList);
84  // ------------------------------
85 
86  return fullGraph;
87 }
88 
112 int primitive::findRings(Graph *fullGraph, int v, std::vector<int> *visited,
113  int maxDepth, int depth, int root) {
114  //
115  int nnumNeighbours; // Number of nearest neighbours for iNode
116  int n; // Node of a neighbour (index)
117  // -------------
118  // For the first call:
119  if (root == -1) {
120  root = v; // Init the root as the current node
121  fullGraph->pts[v].inGraph = false; // Remove the root from the graph
122  } // first call
123  //
124  // Checks
125  if (depth >= maxDepth) {
126  return 1; // false?
127  } // searched until the maximum length specified
128  //
129  // Add the current node to the visited vector
130  visited->push_back(v);
131  // -------------
132  // Depth-first search
133  // 1. Pick a root node (above)
134  // 2. Search all the neighbouring nodes
135  // 3. Start a recursion call, from the second nearest neighbours
136  // 3(i) If the neighbour equals the root (selected in 1), ring has been found!
137  // 3(ii) Or else search for all the unvisited neighbours
138  // 4. Remove the root node from the graph
139  // 5. Update the vector of vectors containing all rings
140  // -------------
141  // Start the search
142  depth += 1; // Go to the next layer
143  nnumNeighbours = fullGraph->pts[v].neighListIndex.size();
144  // Loop through the neighbours of iNode
145  for (int j = 0; j < nnumNeighbours; j++) {
146  n = fullGraph->pts[v].neighListIndex[j]; // Neighbour index
147  // Has a ring been found?!
148  if (depth > 2 && n == root) {
149  // Add the visited vector to the rings vector of vector
150  fullGraph->rings.push_back(*visited);
151  } // A ring has been found!
152  // Otherwise search all the neighbours which have not been searched already
153  else if (fullGraph->pts[n].inGraph) {
154  fullGraph->pts[n].inGraph = false; // Set to false now
155  // Recursive call
156  primitive::findRings(fullGraph, n, visited, maxDepth, depth, root);
157  fullGraph->pts[n].inGraph = true; // Put n back in the graph
158  } // Search the other unsearched neighbours
159  } // end of search through all neighbours
160  // -------------
161  // When the depth is 2, remove just the root from the neighbours of iNode
162  if (depth == 2) {
163  // Search for root in the neighbour list of v
164  for (int j = 0; j < fullGraph->pts[v].neighListIndex.size(); j++) {
165  if (root == fullGraph->pts[v].neighListIndex[j]) {
166  fullGraph->pts[v].neighListIndex.erase(
167  fullGraph->pts[v].neighListIndex.begin() + j);
168  } // end of erase
169  } // end of search
170  } // remove root not edges
171 
172  //
173  (*visited).pop_back();
174  //
175  return 0;
176 }
177 
189  molSys::PointCloud<molSys::Point<double>, double> *yCloud,
190  std::vector<std::vector<int>> neighHbondList) {
191  //
192  primitive::Graph fullGraph; // Contains all the information of the pointCloud
193  primitive::Vertex iVertex; // The vertex corresponding to a particular point
194  int nnumNeighbours; // Number of nearest neighbours for iatom
195  std::vector<int> iNeigh; // Neighbours of the current iatom
196  int jatomID; // Atom ID of the nearest neighbour
197  int jatomIndex; // Atom index of the nearest neighbour
198  // ------------------------------
199  // Loop through every point in yCloud
200  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
201  // -----
202  // Update the neighbours of iatom (only the indices!)
203  iNeigh.clear(); // init
204  nnumNeighbours = neighHbondList[iatom].size() -
205  1; // The first element is atomID of iatom
206  for (int j = 1; j <= nnumNeighbours; j++) {
207  jatomID = neighHbondList[iatom][j]; // Atom ID
208  // Get the atom index for the vector nearest neighbour list
209  auto it = yCloud->idIndexMap.find(jatomID);
210  if (it != yCloud->idIndexMap.end()) {
211  jatomIndex = it->second;
212  } // found jatomIndex
213  else {
214  std::cerr << "Panic induced!\n";
215  return fullGraph;
216  } // jatomIndex not found
217  // Update iNeigh
218  iNeigh.push_back(jatomIndex);
219  } // end of loop through nearest neighbours
220  // -----
221  // Update iVertex
222  iVertex.atomIndex = iatom;
223  iVertex.neighListIndex = iNeigh;
224  // Add to the Graph object
225  fullGraph.pts.push_back(iVertex);
226  } // end of loop through every iatom
227  // ------------------------------
228 
229  return fullGraph;
230 }
231 
245  //
246  primitive::Graph fullGraph; // Contains all the information of the pointCloud
247  primitive::Vertex iVertex; // The vertex corresponding to a particular point
248  int nnumNeighbours; // Number of nearest neighbours for iatom
249  int iatom, jatom; // Atom index being saved
250  // ------------------------------
251  // Loop through every point in nList
252  for (int i = 0; i < nList.size(); i++) {
253  iatom = nList[i][0]; // Atom index of i
254  // neighListIndex is simply the i^th row of nList
255  //
256  // Update iVertex
257  iVertex.atomIndex = iatom;
258  iVertex.neighListIndex = nList[i];
259  // Add to the Graph object
260  fullGraph.pts.push_back(iVertex);
261  } // end of loop through iatom
262 
263  return fullGraph;
264 }
265 
280  std::vector<std::vector<int>> nList) {
281  //
282  // ------------------------------
283  // Loop through every point in nList
284  for (int i = 0; i < nList.size(); i++) {
285  // neighListIndex is simply the i^th row of nList
286  // Update the neighListIndex list in the graph object
287  fullGraph->pts[i].neighListIndex = nList[i];
288  } // end of loop through iatom
289 
290  return *fullGraph;
291 }
292 
304  //
305  int nVertices = fullGraph->pts.size(); // Number of vertices in the graph
306  int nRings = fullGraph->rings.size(); // Number of rings
307  std::vector<bool> ringsToRemove; // Vector containing the logical values for
308  // removal of the current ring index
309  std::vector<int> currentRing; // Current ring being evaluated
310  int ringSize; // Length of the current ring
311  bool removeRing; // Logical for removing the current ring (true) or not
313  emptyTempRings; // Empty vector of vectors to swap
315  primitiveRings; // Vector of vectors of rings after removing non SP rings
316  int currentV; // Current vertex
317  int currentN; // Current neighbour
318  int dist_r; // Distance over ring
319  int dist_g; // Distance over the entire graph
320  int d_jk; // Absolute difference between j and k
321  std::vector<int> path; // Vector containing a path
322  std::vector<int> visited; // Vector containing the visited points
323  // -------------------
324  // Make sure all the vertices are in the graph before removing non-SP rings
325  for (int iVer = 0; iVer < nVertices; iVer++) {
326  fullGraph->pts[iVer].inGraph = true;
327  } // end of loop through every vertex
328  // -------------------
329  // Loop through every ring
330  for (int iRing = 0; iRing < nRings; iRing++) {
331  currentRing = fullGraph->rings[iRing]; // Current ring
332  ringSize = currentRing.size(); // Length of the current ring
333  removeRing = false; // init
334  // Loop through every j^th vertex
335  for (int jVer = 0; jVer < ringSize; jVer++) {
336  // connect with all other, skip j-j (distance=0) and j-(j+1) (nearest
337  // neighbors)
338  for (int kVer = jVer + 2; kVer < ringSize; kVer++) {
339  // If not remove
340  if (!removeRing) {
341  currentV = currentRing[jVer]; // current 'vertex'
342  currentN = currentRing[kVer]; // Current 'neighbour'
343  d_jk = std::abs(jVer - kVer);
344  dist_r = std::min(d_jk, std::abs(d_jk - ringSize)) +
345  1; // Distance over the ring
346  path.clear(); // init
347  visited.clear();
348  // Call shortest path function
349  primitive::shortestPath(fullGraph, currentV, currentN, &path,
350  &visited, dist_r, 0);
351  dist_g = path.size(); // Length of the path over the graph
352  if (dist_g < dist_r) {
353  removeRing = true;
354  } // Decide whether to keep or remove the ring
355  } // If ring is not to be removed
356  } // end of loop through k^th vertex
357  } // end of loop through j^th vertex
358  // Update bool value for removal of currentRing
359  ringsToRemove.push_back(removeRing);
360  } // end of loop through rings
361  // -------------------
362  // Remove all the rings whose indices are given in the ringsToRemove vector
363  for (int i = 0; i < ringsToRemove.size(); i++) {
364  if (!ringsToRemove[i]) {
365  primitiveRings.push_back(fullGraph->rings[i]);
366  } // updates new copy
367  } // end of loop through ringsToRemove
368  // -------------------
369  // Update the graph rings with the primitiveRings
370  emptyTempRings.swap(fullGraph->rings);
371  fullGraph->rings = primitiveRings;
372  // -------------------
373  return *fullGraph;
374 }
375 
393 int primitive::shortestPath(Graph *fullGraph, int v, int goal,
394  std::vector<int> *path, std::vector<int> *visited,
395  int maxDepth, int depth) {
396  int len_path = 0; // Length of the path
397  int nnumNeighbours; // Number of neighbours for a particular v
398  int n; // Index of the nearest neighbour of v
399  // Start the search for the shortest path
400  if (depth < maxDepth) {
401  depth += 1; // One layer below
402  (*visited).push_back(
403  v); // Add the current vertex to the path (visited points)
404  //
405  if (v == goal) {
406  len_path = (*path).size(); // Path of the length of visited points
407  // If the current path is shorter OR this is the first path found
408  if (depth < len_path || len_path == 0) {
409  (*path) = (*visited);
410  maxDepth = depth;
411  } // Current path is the shortest
412  } // Goal reached
413  // Recursive calls to function
414  else {
415  nnumNeighbours = fullGraph->pts[v].neighListIndex.size();
416  // Search all the neighbours
417  for (int j = 0; j < nnumNeighbours; j++) {
418  n = fullGraph->pts[v].neighListIndex[j]; // Index of nearest neighbour
419  // If n has not already been searched:
420  if (fullGraph->pts[n].inGraph == true) {
421  fullGraph->pts[n].inGraph = false; // Set to false
422  primitive::shortestPath(fullGraph, n, goal, path, visited, maxDepth,
423  depth); // Recursive call
424  fullGraph->pts[n].inGraph = true; // Put back in the graph
425  } // If n is in the graph, call recursively
426  } // End of loop over all neighbours
427  } // Goal not reached
428  //
429  // Pop the vector
430  (*visited).pop_back();
431  } // for depth less than maxDepth
432  return 0;
433 }
434 
441  //
443  std::vector<std::vector<int>> tempRings;
444  tempPts.swap(currentGraph->pts);
445  tempRings.swap(currentGraph->rings);
446  return *currentGraph;
447 }
T clear(T... args)
File for generating shortest-path rings according to the Franzblau algorithm.
Graph restoreEdgesFromIndices(Graph *fullGraph, std::vector< std::vector< int >> nList)
Definition: franzblau.cpp:279
Graph populateGraphFromNListID(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> neighHbondList)
Definition: franzblau.cpp:188
Graph countAllRingsFromIndex(std::vector< std::vector< int >> neighHbondList, int maxDepth)
Creates a vector of vectors of all possible rings.
Definition: franzblau.cpp:60
Graph populateGraphFromIndices(std::vector< std::vector< int >> nList)
Definition: franzblau.cpp:244
Graph removeNonSPrings(Graph *fullGraph)
Removes the non-SP rings, using the Franzblau shortest path criterion.
Definition: franzblau.cpp:303
std::vector< std::vector< int > > ringNetwork(std::vector< std::vector< int >> nList, int maxDepth)
Definition: franzblau.cpp:28
std::vector< int > neighListIndex
This is the index according to pointCloud.
Definition: franzblau.hpp:88
std::vector< std::vector< int > > rings
Definition: franzblau.hpp:108
Graph clearGraph(Graph *currentGraph)
Function for clearing vectors in Graph after multiple usage.
Definition: franzblau.cpp:440
int shortestPath(Graph *fullGraph, int v, int goal, std::vector< int > *path, std::vector< int > *visited, int maxDepth, int depth=1)
Calculates the shortest path.
Definition: franzblau.cpp:393
int findRings(Graph *fullGraph, int v, std::vector< int > *visited, int maxDepth, int depth, int root=-1)
Main function that searches for all rings.
Definition: franzblau.cpp:112
std::vector< Vertex > pts
Definition: franzblau.hpp:105
T min(T... args)
T push_back(T... args)
T size(T... args)
This contains a collection of points; contains information for a particular frame.
Definition: mol_sys.hpp:166
This contains per-particle information.
Definition: mol_sys.hpp:145
This is a per-frame object, containing all the vertices for the particular frame, along with the vect...
Definition: franzblau.hpp:104
This is a collection of elements, for each point, required for graph traversal.
Definition: franzblau.hpp:86
T swap(T... args)