Changeset 8d9c38 for src

Timestamp:

May 7, 2008, 2:46:33 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

fe5589

Parents:

6d35e4

Message:

rewritten CyclicStructureAnalysis(): uses DFS spanning tree and limited BFS from all back edges

The DFS tree is given via a StackClass and from each BackEdge a BFS is started to find the leftatom from rightatom but not by taking BackEdge itself. Hence, it finds the smallest circle and returns it.
Also, it limits itself by constantly redefining MinimumRingSize as the minimum of this and the current RingSize.
Tested to good results on buckyball. And may easily be enhanced to scanning and returning all rings (but so far we only need the minimum ring size).

Location:

src

Files:

• molecules.cpp (modified) (3 diffs)
• molecules.hpp (modified) (1 diff)

src/molecules.cpp

@@ -1572 +1572 @@
 
   // analysis of the cycles (print rings, get minimum cycle length)
-  CyclicStructureAnalysis(out, MinimumRingSize);
+  CyclicStructureAnalysis(out, BackEdgeStack, MinimumRingSize);
 
   *out << Verbose(1) << "Final graph info for each atom is:" << endl;
@@ -1606 +1606 @@
 /** Analyses the cycles found and returns minimum of all cycle lengths.
  * \param *out output stream for debugging
- * \param MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found
+ * \param *BackEdgeStack stack with all back edges found during DFS scan
+ * \param MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
  * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond
  */
-void molecule::CyclicStructureAnalysis(ofstream *out, int &MinimumRingSize)
-{
-  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
-  int LP;
-  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Runner = NULL;
-  bond *Binder = NULL;
+void molecule::CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *  BackEdgeStack, int &MinimumRingSize)
+{
+  atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CyclicStructureAnalysis: **PredecessorList");
+  int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CyclicStructureAnalysis: *ShortestPathList");
+  enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
+  class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring 
+  class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms 
+  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
+  bond *Binder = NULL, *BackEdge = NULL;
   int RingSize, NumCycles;
 
-  // count number of cyclic bonds per atom and push all with 3 cyclic bonds onto stack
-  AtomStack->ClearStack();
-  int *NoCyclicBondsPerAtom = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CyclicStructureAnalysis: *NoCyclicBondsPerAtom");
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    NoCyclicBondsPerAtom[Walker->nr] = 0;
-    // check whether it's connected to cyclic bonds and count per atom
-    // 0 means not part of a cycle, 2 means in a cycle, 3 or more means interconnection site of cycles
-    for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
-      Binder = ListOfBondsPerAtom[Walker->nr][i];
-      NoCyclicBondsPerAtom[Walker->nr] += (int) Binder->Cyclic;
-      if (NoCyclicBondsPerAtom[Walker->nr] == 3)  //push all intersections
-        AtomStack->Push(Walker);
-    }
-  }
-  *out << Verbose(1) << "NoCyclicBondsPerAtom: ";
-  for(int i=0;i<AtomCount;i++) {
-    *out << NoCyclicBondsPerAtom[i] << " ";
-  }
-  *out << endl;
+  // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
+  for (int i=0;i<AtomCount;i++) {
+    PredecessorList[i] = NULL;
+    ShortestPathList[i] = -1;
+    ColorList[i] = white;
+  }
+  
+  *out << Verbose(1) << "Back edge list - ";
+  BackEdgeStack->Output(out);
+  
   *out << Verbose(1) << "Analysing cycles ... " << endl;
-  MinimumRingSize = -1;
+  if ((MinimumRingSize <= 2) || (MinimumRingSize > AtomCount))
+    MinimumRingSize = AtomCount;
   NumCycles = 0;
-  while (!AtomStack->IsEmpty()) {
-    Walker = AtomStack->PopFirst();
-    if (NoCyclicBondsPerAtom[Walker->nr] > 1) {
-      NoCyclicBondsPerAtom[Walker->nr]--; // remove one for being intersection
-      RingSize = 0;
-      *out << Verbose(2) << "Current intersection is " << *Walker << ", expect to find " << NoCyclicBondsPerAtom[Walker->nr] << " cycles." << endl;
+  while (!BackEdgeStack->IsEmpty()) {
+    BackEdge = BackEdgeStack->PopFirst();
+    // this is the target
+    Root = BackEdge->leftatom;  
+    // this is the source point
+    Walker = BackEdge->rightatom; 
+    ShortestPathList[Walker->nr] = 0;
+    BFSStack->ClearStack();  // start with empty BFS stack
+    BFSStack->Push(Walker);
+    TouchedStack->Push(Walker);
+    //*out << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
+    OtherAtom = NULL;
+    while ((Walker != Root) && ((OtherAtom == NULL) || (ShortestPathList[OtherAtom->nr] < MinimumRingSize))) {  // look for Root
+      Walker = BFSStack->PopFirst();
+      //*out << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
       for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
         Binder = ListOfBondsPerAtom[Walker->nr][i];
-        OtherAtom = Binder->GetOtherAtom(Walker);
-        // note down the LowPoint number of this cycle
-        if (NoCyclicBondsPerAtom[OtherAtom->nr] > 1) {
-          LP = OtherAtom->LowpointNr;
-          NoCyclicBondsPerAtom[Walker->nr]--; // walker is start of cycle
-          if (LP != Walker->LowpointNr)
-            *out << Verbose(2) << "Tributary cycle: ... <-> " << Walker->Name;
-          else 
-            *out << Verbose(2) << "Main cycle: ... <-> " << Walker->Name;
-          Root = Walker;  // root acts as predecessor marker so that we don't step back accidentally
-          RingSize = 1;
-          do {
-            for(int j=0;j<NumberOfBondsPerAtom[OtherAtom->nr];j++) {  // search among its bonds for next in cycle (same lowpoint nr)
-              Runner = ListOfBondsPerAtom[OtherAtom->nr][j]->GetOtherAtom(OtherAtom);
-              if (((Runner->LowpointNr == LP) || (Runner->LowpointNr == Walker->LowpointNr)) && (Runner != Root)) {
-                // first check is to stay in the cycle
-                // middle check is allow for getting back into main cycle briefly from tributary cycle (just one step, then while further down stops)
-                // last check is not step back
-                *out << " <-> " << OtherAtom->Name;
-                NoCyclicBondsPerAtom[OtherAtom->nr]--;
-                Root = OtherAtom;
-                OtherAtom = Runner;
-                NoCyclicBondsPerAtom[Root->nr]--;
-                RingSize++;
-                break;
-              }
+        if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
+          OtherAtom = Binder->GetOtherAtom(Walker);
+          //*out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
+          if (ColorList[OtherAtom->nr] == white) {
+            TouchedStack->Push(OtherAtom);
+            ColorList[OtherAtom->nr] = lightgray;
+            PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
+            ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
+            //*out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
+            if (ShortestPathList[OtherAtom->nr] < MinimumRingSize) { // Check for maximum distance
+              //*out << Verbose(3) << "Putting OtherAtom into queue." << endl;
+              BFSStack->Push(OtherAtom);
             }
-          } while ((OtherAtom->LowpointNr == LP) && (Walker != OtherAtom) && (Root->LowpointNr == OtherAtom->LowpointNr));
-          // now check if the LP is different from Walker's, as then there is one more bond to follow
-          if (LP != Walker->LowpointNr) {
-            // find last bond to home base
-            for(int j=0;j<NumberOfBondsPerAtom[OtherAtom->nr];j++)
-              if (ListOfBondsPerAtom[OtherAtom->nr][j]->GetOtherAtom(OtherAtom) == Root) {
-                *out << " <-> " << OtherAtom->Name;
-                RingSize++;
-                NoCyclicBondsPerAtom[OtherAtom->nr]--;
-              }
           } else {
-            // we have made the complete cycle
+            //*out << Verbose(3) << "Not Adding, has already been visited." << endl;
           }
-          *out << " <-> ... with cycle length of " << RingSize << "." << endl;
-          NumCycles++;
-          if ((RingSize < MinimumRingSize) || (MinimumRingSize == -1))
-            MinimumRingSize = RingSize;
+        } else {
+          //*out << Verbose(3) << "Not Visiting, is a back edge." << endl;
         }
       }
-    }
-  }
-
-  // print NoCyclicBondsPerAtom to visually check of all are zero
-  *out << Verbose(1) << "NoCyclicBondsPerAtom: ";
-  for(int i=0;i<AtomCount;i++) {
-    if (NoCyclicBondsPerAtom[i] > 0)
-      cerr << "There was an error in  molecule::CyclicStructureAnalysis!" << endl;
-    *out << NoCyclicBondsPerAtom[i] << " ";
-  }
-  *out << endl;
+      ColorList[Walker->nr] = black;
+      //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
+    }
+    
+    if (Walker == Root) {
+      // now climb back the predecessor list and thus find the cycle members
+      NumCycles++;
+      RingSize = 1;
+      Walker = Root;
+      *out << Verbose(1) << "Found ring contains: ";
+      while (Walker != BackEdge->rightatom) {
+        *out << Walker->Name << " <-> ";
+        Walker = PredecessorList[Walker->nr];
+        RingSize++;
+      }
+      *out << Walker->Name << "  with a length of " << RingSize << "." << endl << endl;
+      if (RingSize < MinimumRingSize)
+        MinimumRingSize = RingSize;
+    } else {
+      *out << Verbose(1) << "No ring with length equal to or smaller than " << MinimumRingSize << " found." << endl;
+    }
+    
+    // now clean the lists
+    while (!TouchedStack->IsEmpty()){
+      Walker = TouchedStack->PopFirst();
+      PredecessorList[Walker->nr] = NULL;
+      ShortestPathList[Walker->nr] = -1;
+      ColorList[Walker->nr] = white;
+    }
+  }
+        
 
   if (MinimumRingSize != -1)
@@ -1712 +1709 @@
     *out << Verbose(1) << "No rings were detected in the molecular structure." << endl;
 
-  Free((void **)&NoCyclicBondsPerAtom, "molecule::CyclicStructureAnalysis: *NoCyclicBondsPerAtom");
-  delete(AtomStack);
+  Free((void **)&PredecessorList, "molecule::CyclicStructureAnalysis: **PredecessorList");
+  Free((void **)&ShortestPathList, "molecule::CyclicStructureAnalysis: **ShortestPathList");
+  Free((void **)&ColorList, "molecule::CyclicStructureAnalysis: **ColorList");
+  delete(BFSStack);
 };
 

src/molecules.hpp

@@ -326 +326 @@
   // Graph analysis
   MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, bool ReturnStack, int &MinimumRingSize);
-  void CyclicStructureAnalysis(ofstream *out, int &MinimumRingSize);
+  void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int &MinimumRingSize);
   bond * FindNextUnused(atom *vertex);
   void SetNextComponentNumber(atom *vertex, int nr);