/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /* * CyclicStructureAnalysis.cpp * * Created on: Feb 16, 2011 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "CyclicStructureAnalysis.hpp" #include "atom.hpp" #include "Bond/bond.hpp" #include "CodePatterns/Assert.hpp" #include "CodePatterns/Info.hpp" #include "CodePatterns/Log.hpp" #include "CodePatterns/Verbose.hpp" #include "Element/element.hpp" #include "molecule.hpp" CyclicStructureAnalysis::CyclicStructureAnalysis() {} CyclicStructureAnalysis::~CyclicStructureAnalysis() {} /** Initialise vertex as white with no predecessor, no shortest path(-1), color white. * \param atom_id id of atom whose node we address */ void CyclicStructureAnalysis::InitNode(atomId_t atom_id) { ShortestPathList[atom_id] = -1; PredecessorList[atom_id] = 0; ColorList[atom_id] = GraphEdge::white; } void CyclicStructureAnalysis::Reset() { // clear what's present ShortestPathList.clear(); PredecessorList.clear(); ColorList.clear(); BFSStack.clear(); TouchedStack.clear(); } /** Clean the accounting structure for all nodes touched so far. */ void CyclicStructureAnalysis::CleanAllTouched() { atom *Walker = NULL; while (!TouchedStack.empty()) { Walker = TouchedStack.front(); TouchedStack.pop_front(); PredecessorList[Walker->getNr()] = NULL; ShortestPathList[Walker->getNr()] = -1; ColorList[Walker->getNr()] = GraphEdge::white; } } /** Resets shortest path list and BFSStack. * \param *&Walker current node, pushed onto BFSStack and TouchedStack */ void CyclicStructureAnalysis::InitializeToRoot(atom *&Root) { ShortestPathList[Root->getNr()] = 0; BFSStack.clear(); // start with empty BFS stack BFSStack.push_front(Root); TouchedStack.push_front(Root); } /** Performs a BFS from \a *Root, trying to find the same node and hence a cycle. * \param *&BackEdge the edge from root that we don't want to move along * \param &BFS accounting structure */ void CyclicStructureAnalysis::CyclicBFSFromRootToRoot(bond *&BackEdge) { atom *Walker = NULL; atom *OtherAtom = NULL; do { // look for Root ASSERT(!BFSStack.empty(), "CyclicStructureAnalysis_CyclicBFSFromRootToRoot() - BFSStack is empty!"); Walker = BFSStack.front(); BFSStack.pop_front(); LOG(2, "INFO: Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl); const BondList& ListOfBonds = Walker->getListOfBonds(); for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) { if ((*Runner) != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder) OtherAtom = (*Runner)->GetOtherAtom(Walker); #ifdef ADDHYDROGEN if (OtherAtom->getType()->getAtomicNumber() != 1) { #endif LOG(2, "INFO: Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << "." << endl); if (ColorList[OtherAtom->getNr()] == GraphEdge::white) { TouchedStack.push_front(OtherAtom); ColorList[OtherAtom->getNr()] = GraphEdge::lightgray; PredecessorList[OtherAtom->getNr()] = Walker; // Walker is the predecessor ShortestPathList[OtherAtom->getNr()] = ShortestPathList[Walker->getNr()] + 1; LOG(2, "INFO: Coloring OtherAtom " << OtherAtom->getName() << " lightgray, its predecessor is " << Walker->getName() << " and its Shortest Path is " << ShortestPathList[OtherAtom->getNr()] << " egde(s) long." << endl); //if (ShortestPathList[OtherAtom->getNr()] < MinimumRingSize[Walker->GetTrueFather()->getNr()]) { // Check for maximum distance LOG(3, "ACCEPT: Putting OtherAtom into queue." << endl); BFSStack.push_front(OtherAtom); //} } else { LOG(3, "REJECT: Not Adding, has already been visited." << endl); } if (OtherAtom == Root) break; #ifdef ADDHYDROGEN } else { LOG(2, "INFO: Skipping hydrogen atom " << *OtherAtom << "." << endl); ColorList[OtherAtom->getNr()] = GraphEdge::black; } #endif } else { LOG(2, "REJECT: Bond " << *(*Runner) << " not Visiting, is the back edge." << endl); } } ColorList[Walker->getNr()] = GraphEdge::black; LOG(1, "INFO: Coloring Walker " << Walker->getName() << " " << GraphEdge::getColorName(ColorList[Walker->getNr()]) << "." << endl); if (OtherAtom == Root) { // if we have found the root, check whether this cycle wasn't already found beforehand // step through predecessor list while (OtherAtom != BackEdge->rightatom) { if (!OtherAtom->GetTrueFather()->IsCyclic) // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet break; else OtherAtom = PredecessorList[OtherAtom->getNr()]; } if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already LOG(3, "INFO This cycle was already found before, skipping and removing seeker from search." << endl); do { ASSERT(!TouchedStack.empty(), "CyclicStructureAnalysis_CyclicBFSFromRootToRoot() - TouchedStack is empty!"); OtherAtom = TouchedStack.front(); TouchedStack.pop_front(); if (PredecessorList[OtherAtom->getNr()] == Walker) { LOG(4, "INFO: Removing " << *OtherAtom << " from lists and stacks." << endl); PredecessorList[OtherAtom->getNr()] = NULL; ShortestPathList[OtherAtom->getNr()] = -1; ColorList[OtherAtom->getNr()] = GraphEdge::white; // rats ... deque has no find() std::deque::iterator iter = find( BFSStack.begin(), BFSStack.end(), OtherAtom); ASSERT(iter != BFSStack.end(), "CyclicStructureAnalysis_CyclicBFSFromRootToRoot() - can't find "+toString(*OtherAtom)+" on stack!"); BFSStack.erase(iter); } } while ((!TouchedStack.empty()) && (PredecessorList[OtherAtom->getNr()] == NULL)); TouchedStack.push_front(OtherAtom); // last was wrongly popped OtherAtom = BackEdge->rightatom; // set to not Root } else OtherAtom = Root; } } while ((!BFSStack.empty()) && (OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->getNr()] < MinimumRingSize[Walker->GetTrueFather()->getNr()]))); } /** Climb back the BFSAccounting::PredecessorList and find cycle members. * \param *&OtherAtom * \param *&BackEdge denotes the edge we did not want to travel along when doing CyclicBFSFromRootToRoot() * \param &BFS accounting structure * \param &MinRingSize global minimum distance from one node without encountering oneself, set on return */ void CyclicStructureAnalysis::RetrieveCycleMembers(atom *&OtherAtom, bond *&BackEdge, int &MinRingSize) { atom *Walker = NULL; int NumCycles = 0; int RingSize = -1; if (OtherAtom == Root) { // now climb back the predecessor list and thus find the cycle members NumCycles++; RingSize = 1; Root->GetTrueFather()->IsCyclic = true; std::stringstream output; output << "Found ring contains: "; Walker = Root; while (Walker != BackEdge->rightatom) { output << Walker->getName() << " <-> "; Walker = PredecessorList[Walker->getNr()]; Walker->GetTrueFather()->IsCyclic = true; RingSize++; } output << Walker->getName() << " with a length of " << RingSize << "."; LOG(0, "INFO: " << output.str()); // walk through all and set MinimumRingSize Walker = Root; ASSERT(!MinimumRingSize.count(Walker->GetTrueFather()->getNr()), "CyclicStructureAnalysis::RetrieveCycleMembers() - setting MinimumRingSize of " +toString(*(Walker->GetTrueFather()))+" to " +toString(RingSize)+" which is already set to " +toString(MinimumRingSize[Walker->GetTrueFather()->getNr()])+"."); MinimumRingSize[Walker->GetTrueFather()->getNr()] = RingSize; while (Walker != BackEdge->rightatom) { Walker = PredecessorList[Walker->getNr()]; if (RingSize < MinimumRingSize[Walker->GetTrueFather()->getNr()]) MinimumRingSize[Walker->GetTrueFather()->getNr()] = RingSize; } if ((RingSize < MinRingSize) || (MinRingSize == -1)) MinRingSize = RingSize; } else { LOG(1, "INFO: No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Root->GetTrueFather()->getNr()] << " found." << endl); } } /** From a given node performs a BFS to touch the next cycle, for whose nodes \a MinimumRingSize is set and set it accordingly. * \param *&Root node to look for closest cycle from, i.e. \a MinimumRingSize is set for this node * \param AtomCount number of nodes in graph */ void CyclicStructureAnalysis::BFSToNextCycle(atom *&Root, atom *&Walker) { atom *OtherAtom = Walker; Reset(); InitializeToRoot(Walker); while (OtherAtom != NULL) { // look for Root ASSERT(!BFSStack.empty(), "CyclicStructureAnalysis_BFSToNextCycle() - BFSStack is empty!"); Walker = BFSStack.front(); BFSStack.pop_front(); LOG(2, "INFO: Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "."); const BondList& ListOfBonds = Walker->getListOfBonds(); for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) { // "removed (*Runner) != BackEdge) || " from next if, is u if ((ListOfBonds.size() == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check OtherAtom = (*Runner)->GetOtherAtom(Walker); LOG(2, "INFO: Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << "."); if (ColorList[OtherAtom->getNr()] == GraphEdge::white) { TouchedStack.push_front(OtherAtom); ColorList[OtherAtom->getNr()] = GraphEdge::lightgray; PredecessorList[OtherAtom->getNr()] = Walker; // Walker is the predecessor ShortestPathList[OtherAtom->getNr()] = ShortestPathList[Walker->getNr()] + 1; LOG(2, "ACCEPT: Coloring OtherAtom " << OtherAtom->getName() << " lightgray, its predecessor is " << Walker->getName() << " and its Shortest Path is " << ShortestPathList[OtherAtom->getNr()] << " egde(s) long."); if (OtherAtom->GetTrueFather()->IsCyclic) { // if the other atom is connected to a ring ASSERT(!MinimumRingSize.count(Root->GetTrueFather()->getNr()), "CyclicStructureAnalysis::BFSToNextCycle() - setting MinimumRingSize of " +toString(*(Root->GetTrueFather()))+" to "+ toString(ShortestPathList[OtherAtom->getNr()] + MinimumRingSize[OtherAtom->GetTrueFather()->getNr()]) +" which is already set to " +toString(MinimumRingSize[Root->GetTrueFather()->getNr()])+"."); MinimumRingSize[Root->GetTrueFather()->getNr()] = ShortestPathList[OtherAtom->getNr()] + MinimumRingSize[OtherAtom->GetTrueFather()->getNr()]; OtherAtom = NULL; //break; break; } else BFSStack.push_front(OtherAtom); } else { LOG(3, "REJECT: Not Adding, has already been visited."); } } else { LOG(3, "REJECT: Not Visiting, is a back edge."); } } ColorList[Walker->getNr()] = GraphEdge::black; LOG(1, "INFO: Coloring Walker " << Walker->getName() << " " << GraphEdge::getColorName(ColorList[Walker->getNr()]) << "."); } } /** All nodes that are not in cycles get assigned a \a *&MinimumRingSizeby BFS to next cycle. * \param *&MinimumRingSize array with minimum distance without encountering onself for each atom * \param &MinRingSize global minium distance * \param &NumCyles number of cycles in graph */ void CyclicStructureAnalysis::AssignRingSizetoNonCycleMembers(int &MinRingSize, int &NumCycles) { atom *Root = NULL; atom *Walker = NULL; if (MinRingSize != -1) { // if rings are present // go over all atoms World::AtomComposite allatoms = World::getInstance().getAllAtoms(); for (World::AtomComposite::const_iterator iter = allatoms.begin(); iter != allatoms.end(); ++iter) { Root = *iter; if (MinimumRingSize.find(Root->GetTrueFather()->getNr()) != MinimumRingSize.end()) { // check whether MinimumRingSize is set, if not BFS to next where it is Walker = Root; LOG(1, "---------------------------------------------------------------------------------------------------------"); BFSToNextCycle(Root, Walker); } ASSERT(MinimumRingSize.find(Root->GetTrueFather()->getNr()) != MinimumRingSize.end(), "CyclicStructureAnalysis::AssignRingSizetoNonCycleMembers() - BFSToNextCycle did not set MinimumRingSize of " +toString(*(Root->GetTrueFather()))+"."); LOG(1, "INFO: Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->getNr()] << "." << endl); } LOG(1, "INFO: Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl); } else LOG(1, "INFO: No rings were detected in the molecular structure." << endl); } /** Analyses the cycles found and returns minimum of all cycle lengths. * We begin with a list of Back edges found during DepthFirstSearchAnalysis(). We go through this list - one end is the Root, * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds * as cyclic and print out the cycles. * \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! * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond */ void CyclicStructureAnalysis::operator()(std::deque * BackEdgeStack) { Info FunctionInfo("CyclicStructureAnalysis"); atom *Walker = NULL; atom *OtherAtom = NULL; bond *BackEdge = NULL; int NumCycles = 0; int MinRingSize = -1; //Log() << Verbose(1) << "Back edge list - "; //BackEdgeStack->Output(out); LOG(1, "STATUS: Analysing cycles ... " << endl); NumCycles = 0; while (!BackEdgeStack->empty()) { BackEdge = BackEdgeStack->front(); BackEdgeStack->pop_front(); // this is the target Root = BackEdge->leftatom; // this is the source point Walker = BackEdge->rightatom; InitializeToRoot(Walker); LOG(1, "---------------------------------------------------------------------------------------------------------" << endl); OtherAtom = NULL; // go to next cycle via BFS CyclicBFSFromRootToRoot(BackEdge); // get all member nodes of this cycle RetrieveCycleMembers(OtherAtom, BackEdge, MinRingSize); CleanAllTouched(); } AssignRingSizetoNonCycleMembers(MinRingSize, NumCycles); } /** Output a list of flags, stating whether the bond was visited or not. * \param *list list to print */ void CyclicStructureAnalysis::OutputAlreadyVisited(int *list) { std::stringstream output; output << "Already Visited Bonds:\t"; for (int i = 1; i <= list[0]; i++) output << list[i] << " "; LOG(0, output.str()); } const std::map& CyclicStructureAnalysis::getMinimumRingSize() const { return MinimumRingSize; }