Changeset b86de7 for molecuilder/src/molecules.cpp

Timestamp:

Oct 9, 2008, 6:27:56 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

d01f4d

Parents:

a7b3b8

Message:

BUGFIXES: CyclicStructureAnalysis() now compatible to disconnected subgraphs, AssignKeySetsToFragment() and FillBondStructureFromReference() memory cleanup corrected

+ molecule::DepthFirstSearchAnalysis() now just returns BackEdgeStack, not MinimumRingSize. CyclicStructureAnalysis() is called during FragmentMolecule(), after subgraphs bonds list have been filled by FillBondStructureFromReference().
+ new function molecule::PickLocalBackEdges(), as the BackEdgeStack returned by DepthFirstSearchAnalysis() co
ntains only global bonds, not the local ones for the subgraph, we have to step through it and pick the right
ones out.
+ molecule::FragmentMolecule() now calls molecule::CyclicStructureAnalysis() separately for each subgraph, along with a BackEdgeStack filled by PickLocalBackEdges(), and allocates&initialises MinimumRingSize array. Als
o AssignKeySetsToFragment() frees the LocalListOfAtoms now (FreeList=true), now longer after the following wh
ile
+ molecule::CyclicStructureAnalysis() takes a local BackEdgeStack and analysis the subgraphs cycles, returnin
g minimum ring size
+ MoleculeLeafClass::AssignKeySetsToFragment() now frees memory for ListOfLocalAtoms when FreeList is set. BUGFIX: test of first key was testing against ..->nr != -1. However, LocalListOfAtoms was not even initialised correctly to NULL, hence ...->nr pointed in some cases to nowhere. Now it test atom* against NULL.
+ MoleculeLeafClass::FillBondStructureFromReference() now frees memory for ListOfLocalAtoms when FreeList is set correctly (only free initial pointer when FragmentCounter == 0, also it was decreased not before but after freeing, hence we free'd the wrong list). Also, father replaced by GetTrueFather() (makes the function moregenerally useable, was not a bug).
+ ParseCommandLineOptions() option 'D': adapted to changes in DepthFirstSearchAnalysis() in a similar manner
to FragmentMolecule()
+ molecule::CountCyclicBonds() adapted but does not perform CyclicStructureAnalysis()
+ molecule::CreateAdjacencyList() counts the bonds that could not be brought to covalently corrected degree (i.e. the remaining ionic atoms)
+ molecule::CreateListOfBondsPerAtom() prints atom names and number, which is helpful as name contains global

and number contains local number (helped in finding above bugs)

+ CreateFatherLookupTable(): BUGFIX: LookupTable was not initialised to NULL (see above)

File:

• molecuilder/src/molecules.cpp (modified) (20 diffs)

molecuilder/src/molecules.cpp

@@ -2244 +2244 @@
   int *MinimumRingSize = NULL;
   MoleculeLeafClass *Subgraphs = NULL;
+  class StackClass<bond *> *BackEdgeStack = NULL;
   bond *Binder = first;
   if ((Binder->next != last) && (Binder->next->Type == Undetermined)) {
     *out << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl;
-    Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
+    Subgraphs = DepthFirstSearchAnalysis(out, BackEdgeStack);
     while (Subgraphs->next != NULL) {
       Subgraphs = Subgraphs->next;
@@ -2260 +2261 @@
       No++;    
   }
+  delete(BackEdgeStack);
   return No;
 };
@@ -2341 +2343 @@
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
   Vector x;
+  int FalseBondDegree = 0;
  
   BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
@@ -2359 +2362 @@
       j += i+1;
       divisor[i] = (int)floor(cell_size[j]/bonddistance); // take smaller value such that size of linked cell is at least bonddistance
-      *out << Verbose(1) << "divisor[" << i << "]  = " << divisor[i] << "." << endl;
+      //*out << Verbose(1) << "divisor[" << i << "]  = " << divisor[i] << "." << endl;
     }
     // 2a. allocate memory for the cell list
@@ -2384 +2387 @@
       }
       index = n[2] + (n[1] + n[0] * divisor[1]) * divisor[2];
-      *out << Verbose(1) << "Atom " << *Walker << " goes into cell number [" << n[0] << "," << n[1] << "," << n[2] << "] = " << index << "." << endl;
+      //*out << Verbose(1) << "Atom " << *Walker << " goes into cell number [" << n[0] << "," << n[1] << "," << n[2] << "] = " << index << "." << endl;
       // add copy atom to this cell
       if (CellList[index] == NULL)  // allocate molecule if not done
@@ -2426 +2429 @@
                         distance =   OtherWalker->x.PeriodicDistance(&(Walker->x), cell_size);
                         if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
-                          *out << Verbose(0) << "Adding Bond between " << *Walker << " and " << *OtherWalker << "." << endl;
+                          //*out << Verbose(0) << "Adding Bond between " << *Walker << " and " << *OtherWalker << "." << endl;
                           AddBond(Walker->father, OtherWalker->father, 1);  // also increases molecule::BondCount
                           BondCount++;
@@ -2485 +2488 @@
               ListOfBondsPerAtom[Walker->nr][CandidateBondNo]->BondDegree++;
               *out << Verbose(2) << "Increased bond degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
-            }
+            } else
+              *out << Verbose(2) << "Could not find correct degree for bond " << *ListOfBondsPerAtom[Walker->nr][CandidateBondNo] << "." << endl;
+              FalseBondDegree++;
           }
                     }
@@ -2492 +2497 @@
                 } else
                         *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
-          *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << "." << endl;
+          *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
                 
           // output bonds for debugging (if bond chain list was correctly installed)
@@ -2513 +2518 @@
  * We use the algorithm from [Even, Graph Algorithms, p.62].
  * \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 NULL pointer to StackClass with all the found back edges, allocated and filled on return
  * \return list of each disconnected subgraph as an individual molecule class structure
  */
-MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, int *&MinimumRingSize)
-{
-  class StackClass<atom *> *AtomStack;
-  AtomStack = new StackClass<atom *>(AtomCount);
-  class StackClass<bond *> *BackEdgeStack = new StackClass<bond *> (BondCount);
+MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack)
+{
+  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
+  BackEdgeStack = new StackClass<bond *> (BondCount);
   MoleculeLeafClass *SubGraphs = new MoleculeLeafClass(NULL);
   MoleculeLeafClass *LeafWalker = SubGraphs;
@@ -2652 +2656 @@
     LeafWalker->Leaf->Output(out);
     *out << endl;
-    
+
     // step on to next root
     while ((Root != end) && (Root->GraphNr != -1)) {
@@ -2670 +2674 @@
   }
 
-  // analysis of the cycles (print rings, get minimum cycle length)
-  CyclicStructureAnalysis(out, BackEdgeStack, MinimumRingSize);
 
   *out << Verbose(1) << "Final graph info for each atom is:" << endl;
@@ -2709 +2711 @@
  * as cyclic and print out the cycles. 
  * \param *out output stream for debugging
- * \param *BackEdgeStack stack with all back edges found during DFS scan
+ * \param *BackEdgeStack stack with all back edges found during DFS scan. Beware: This stack contains the bonds from the total molecule, not from the subgraph!
  * \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
@@ -2730 +2732 @@
     ColorList[i] = white;
   }
-  MinimumRingSize = new int[AtomCount];
-  for(int i=AtomCount;i--;)
-    MinimumRingSize[i] = AtomCount;
-
   
   *out << Verbose(1) << "Back edge list - ";
@@ -3535 +3533 @@
         MoleculeListClass *BondFragments = NULL;
   int *SortIndex = NULL;
-  int *MinimumRingSize = NULL;
+  int *MinimumRingSize = new int[AtomCount];
   int FragmentCounter;
   MoleculeLeafClass *MolecularWalker = NULL;
@@ -3541 +3539 @@
   fstream File;
   bool FragmentationToDo = true;
+  class StackClass<bond *> *BackEdgeStack = NULL, *LocalBackEdgeStack = NULL;
   bool CheckOrder = false;
   Graph **FragmentList = NULL;
@@ -3572 +3571 @@
 
   // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
-  Subgraphs = DepthFirstSearchAnalysis(out, MinimumRingSize);
+  Subgraphs = DepthFirstSearchAnalysis(out, BackEdgeStack);
   // fill the bond structure of the individually stored subgraphs
   Subgraphs->next->FillBondStructureFromReference(out, this, (FragmentCounter = 0), ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
-
+  // analysis of the cycles (print rings, get minimum cycle length) for each subgraph
+  for(int i=AtomCount;i--;)
+    MinimumRingSize[i] = AtomCount;
+  MolecularWalker = Subgraphs;
+  FragmentCounter = 0;
+  while (MolecularWalker->next != NULL) {
+    MolecularWalker = MolecularWalker->next;
+    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
+    LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
+//    // check the list of local atoms for debugging
+//    *out << Verbose(0) << "ListOfLocalAtoms for this subgraph is:" << endl;
+//    for (int i=0;i<AtomCount;i++)
+//      if (ListOfLocalAtoms[FragmentCounter][i] == NULL)
+//        *out << "\tNULL";
+//      else
+//        *out << "\t" << ListOfLocalAtoms[FragmentCounter][i]->Name;
+    MolecularWalker->Leaf->PickLocalBackEdges(out, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
+    MolecularWalker->Leaf->CyclicStructureAnalysis(out, LocalBackEdgeStack, MinimumRingSize);
+    delete(LocalBackEdgeStack);
+  }
+  
   // ===== 3. if structure still valid, parse key set file and others =====
   FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
@@ -3584 +3603 @@
   // =================================== Begin of FRAGMENTATION =============================== 
   // ===== 6a. assign each keyset to its respective subgraph ===== 
-  Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), false);
+  Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
 
   // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
@@ -3621 +3640 @@
   delete[](MinimumRingSize);
   
-  // free the index lookup list
-  for (int i=FragmentCounter;i--;)
-    Free((void **)&ListOfLocalAtoms[i], "molecule::FragmentMolecule - *ListOfLocalAtoms[]");
-  Free((void **)&ListOfLocalAtoms, "molecule::FragmentMolecule - **ListOfLocalAtoms");
 
   // ==================================== End of FRAGMENTATION ============================================
@@ -3698 +3713 @@
 };
 
+
+/** Picks from a global stack with all back edges the ones in the fragment.
+ * \param *out output stream for debugging
+ * \param **ListOfLocalAtoms array of father atom::nr to local atom::nr (reverse of atom::father)
+ * \param *ReferenceStack stack with all the back egdes
+ * \param *LocalStack stack to be filled
+ * \return true - everything ok, false - ReferenceStack was empty
+ */
+bool molecule::PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack)
+{
+  bool status = true;
+  if (ReferenceStack->IsEmpty()) {
+    cerr << "ReferenceStack is empty!" << endl;
+    return false;
+  }
+  bond *Binder = ReferenceStack->PopFirst();
+  bond *FirstBond = Binder;   // mark the first bond, so that we don't loop through the stack indefinitely
+  atom *Walker = NULL, *OtherAtom = NULL;
+  ReferenceStack->Push(Binder);
+  
+  do {  // go through all bonds and push local ones
+    Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
+    if (Walker == NULL) // if this Walker exists in the subgraph ... 
+      continue;
+    for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
+      OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+      if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
+        LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
+        break;
+      }
+    }
+    Binder = ReferenceStack->PopFirst();  // loop the stack for next item
+    ReferenceStack->Push(Binder);
+  } while (FirstBond != Binder);
+  
+  return status;
+};
+
 /** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
  * Atoms not present in the file get "-1".
@@ -3840 +3893 @@
   while (Walker->next != end) {
     Walker = Walker->next; 
-    *out << Verbose(4) << "Atom " << Walker->Name << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
+    *out << Verbose(4) << "Atom " << Walker->Name << "/" << Walker->nr << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
     TotalDegree = 0;
     for (int j=0;j<NumberOfBondsPerAtom[Walker->nr];j++) {