/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010-2012 University of Bonn. All rights reserved. * Copyright (C) 2013 Frederik Heber. All rights reserved. * * * This file is part of MoleCuilder. * * MoleCuilder is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 2 of the License, or * (at your option) any later version. * * MoleCuilder is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MoleCuilder. If not, see . */ /* * DepthFirstSearchAnalysis.cpp * * Created on: Feb 16, 2011 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "DepthFirstSearchAnalysis.hpp" #include #include #include "Atom/atom.hpp" #include "Bond/bond.hpp" #include "CodePatterns/Assert.hpp" #include "CodePatterns/Info.hpp" #include "CodePatterns/Log.hpp" #include "CodePatterns/Verbose.hpp" #include "Descriptors/AtomDescriptor.hpp" #include "Descriptors/MoleculeDescriptor.hpp" #include "Graph/ListOfLocalAtoms.hpp" #include "molecule.hpp" #include "MoleculeLeafClass.hpp" #include "World.hpp" DepthFirstSearchAnalysis::DepthFirstSearchAnalysis() : CurrentGraphNr(0), ComponentNumber(0), BackStepping(false) { ResetAllBondsToUnused(); } DepthFirstSearchAnalysis::~DepthFirstSearchAnalysis() {} void DepthFirstSearchAnalysis::Init() { CurrentGraphNr = 0; ComponentNumber = 0; BackStepping = false; std::for_each(World::getInstance().getAtomIter(),World::getInstance().atomEnd(), std::mem_fun(&atom::resetGraphNr)); std::for_each(World::getInstance().getAtomIter(),World::getInstance().atomEnd(), std::mem_fun(&atom::InitComponentNr)); } bond::ptr DepthFirstSearchAnalysis::FindNextUnused(atom *vertex) const { const BondList& ListOfBonds = vertex->getListOfBonds(); for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) if ((*Runner)->IsUsed() == GraphEdge::white) return ((*Runner)); return bond::ptr(); } void DepthFirstSearchAnalysis::ResetAllBondsToUnused() const { World::AtomComposite allatoms = World::getInstance().getAllAtoms(); for(World::AtomComposite::const_iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) { const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds(); for(BondList::const_iterator BondRunner = ListOfBonds.begin(); BondRunner != ListOfBonds.end(); ++BondRunner) if ((*BondRunner)->leftatom == *AtomRunner) (*BondRunner)->ResetUsed(); } } void DepthFirstSearchAnalysis::SetNextComponentNumber(atom *vertex, int nr) const { size_t i = 0; ASSERT(vertex != NULL, "DepthFirstSearchAnalysis::SetNextComponentNumber() - Given vertex is NULL!"); const BondList& ListOfBonds = vertex->getListOfBonds(); for (; i < ListOfBonds.size(); i++) { if (vertex->ComponentNr[i] == -1) { // check if not yet used vertex->ComponentNr[i] = nr; break; } else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time break; // breaking here will not cause error! } ASSERT(i < ListOfBonds.size(), "DepthFirstSearchAnalysis::SetNextComponentNumber() - All Component entries are already occupied!"); } bool DepthFirstSearchAnalysis::PickLocalBackEdges(const ListOfLocalAtoms_t &ListOfLocalAtoms, std::deque *&LocalStack) const { bool status = true; if (BackEdgeStack.empty()) { ELOG(1, "Reference BackEdgeStack is empty!"); return false; } std::deque MyBackEdgeStack = BackEdgeStack; do { // go through all bonds and push local ones const bond::ptr &Binder = MyBackEdgeStack.front(); // loop the stack for next item MyBackEdgeStack.pop_front(); LOG(3, "INFO: Current candidate edge " << *Binder << "."); const ListOfLocalAtoms_t::const_iterator leftiter = ListOfLocalAtoms.find(Binder->leftatom->getNr()); ASSERT( leftiter != ListOfLocalAtoms.end(), "DepthFirstSearchAnalysis::PickLocalBackEdges() - could not find atom id " +toString(Binder->leftatom->getNr())+" in ListOfLocalAtoms."); atom * const Walker = leftiter->second; // get one atom in the reference molecule if (Walker != NULL) { // if this Walker exists in the subgraph ... const BondList& ListOfBonds = Walker->getListOfBonds(); for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) { atom * const OtherAtom = (*Runner)->GetOtherAtom(Walker); const ListOfLocalAtoms_t::const_iterator rightiter = ListOfLocalAtoms.find((*Runner)->rightatom->getNr()); if (OtherAtom == rightiter->second) { // found the bond LocalStack->push_front((*Runner)); LOG(3, "INFO: Found local edge " << *(*Runner) << "."); break; } } } } while (!MyBackEdgeStack.empty()); return status; } void DepthFirstSearchAnalysis::OutputGraphInfoPerAtom() const { LOG(1, "Final graph info for each atom is:"); World::ConstAtomComposite allatoms = const_cast(World::getInstance()). getAllAtoms(); for_each(allatoms.begin(),allatoms.end(),mem_fun(&atom::OutputGraphInfo)); } void DepthFirstSearchAnalysis::OutputGraphInfoPerBond() const { LOG(1, "Final graph info for each bond is:"); World::ConstAtomComposite allatoms = const_cast(World::getInstance()). getAllAtoms(); for(World::ConstAtomComposite::const_iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) { const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds(); for(BondList::const_iterator BondRunner = ListOfBonds.begin(); BondRunner != ListOfBonds.end(); ++BondRunner) if ((*BondRunner)->leftatom == *AtomRunner) { const bond::ptr Binder = *BondRunner; if (DoLog(2)) { std::stringstream output; output << ((Binder->Type == GraphEdge::TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <"; output << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp."; Binder->leftatom->OutputComponentNumber(&output); output << " === "; output << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp."; Binder->rightatom->OutputComponentNumber(&output); output << ">."; LOG(2, output.str()); } if (Binder->Cyclic) // cyclic ?? LOG(3, "Lowpoint at each side are equal: CYCLIC!"); } } } unsigned int DepthFirstSearchAnalysis::CyclicBondAnalysis() const { unsigned int NoCyclicBonds = 0; World::ConstAtomComposite allatoms = const_cast(World::getInstance()). getAllAtoms(); for(World::ConstAtomComposite::const_iterator AtomRunner = allatoms.begin(); AtomRunner != allatoms.end(); ++AtomRunner) { const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds(); for(BondList::const_iterator BondRunner = ListOfBonds.begin(); BondRunner != ListOfBonds.end(); ++BondRunner) if ((*BondRunner)->leftatom == *AtomRunner) if ((*BondRunner)->rightatom->LowpointNr == (*BondRunner)->leftatom->LowpointNr) { // cyclic ?? (*BondRunner)->Cyclic = true; NoCyclicBonds++; } } return NoCyclicBonds; } void DepthFirstSearchAnalysis::SetWalkersGraphNr(atom *&Walker) { if (!BackStepping) { // if we don't just return from (8) Walker->GraphNr = CurrentGraphNr; Walker->LowpointNr = CurrentGraphNr; LOG(1, "Setting Walker[" << Walker->getName() << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << "."); AtomStack.push_front(Walker); CurrentGraphNr++; } } void DepthFirstSearchAnalysis::ProbeAlongUnusedBond(atom *&Walker, bond::ptr &Binder) { atom *OtherAtom = NULL; do { // (3) if Walker has no unused egdes, go to (5) BackStepping = false; // reset backstepping flag for (8) if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused Binder = FindNextUnused(Walker); if (Binder == NULL) break; LOG(2, "Current Unused Bond is " << *Binder << "."); // (4) Mark Binder used, ... Binder->MarkUsed(GraphEdge::black); OtherAtom = Binder->GetOtherAtom(Walker); LOG(2, "(4) OtherAtom is " << OtherAtom->getName() << "."); if (OtherAtom->GraphNr != -1) { // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3) Binder->Type = GraphEdge::BackEdge; BackEdgeStack.push_front(Binder); Walker->LowpointNr = (Walker->LowpointNr < OtherAtom->GraphNr) ? Walker->LowpointNr : OtherAtom->GraphNr; LOG(3, "(4a) Visited: Setting Lowpoint of Walker[" << Walker->getName() << "] to " << Walker->LowpointNr << "."); } else { // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2) Binder->Type = GraphEdge::TreeEdge; OtherAtom->Ancestor = Walker; Walker = OtherAtom; LOG(3, "(4b) Not Visited: OtherAtom[" << OtherAtom->getName() << "]'s Ancestor is now " << OtherAtom->Ancestor->getName() << ", Walker is OtherAtom " << OtherAtom->getName() << "."); break; } Binder.reset(); } while (1); // (3) } void DepthFirstSearchAnalysis::CheckForaNewComponent(atom *&Walker, ConnectedSubgraph &Subgraph) { atom *OtherAtom = NULL; // (5) if Ancestor of Walker is ... LOG(1, "(5) Number of Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "] is " << Walker->Ancestor->GraphNr << "."); if (Walker->Ancestor->GraphNr != Root->GraphNr) { // (6) (Ancestor of Walker is not Root) if (Walker->LowpointNr < Walker->Ancestor->GraphNr) { // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8) Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr; LOG(2, "(6) Setting Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << "."); } else { // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component Walker->Ancestor->SeparationVertex = true; LOG(2, "(7) Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s is a separating vertex, creating component."); SetNextComponentNumber(Walker->Ancestor, ComponentNumber); LOG(3, "(7) Walker[" << Walker->getName() << "]'s Ancestor's Compont is " << ComponentNumber << "."); SetNextComponentNumber(Walker, ComponentNumber); LOG(3, "(7) Walker[" << Walker->getName() << "]'s Compont is " << ComponentNumber << "."); do { ASSERT(!AtomStack.empty(), "DepthFirstSearchAnalysis_CheckForaNewComponent() - AtomStack is empty!"); OtherAtom = AtomStack.front(); AtomStack.pop_front(); Subgraph.push_back(OtherAtom); SetNextComponentNumber(OtherAtom, ComponentNumber); LOG(3, "(7) Other[" << OtherAtom->getName() << "]'s Compont is " << ComponentNumber << "."); } while (OtherAtom != Walker); ComponentNumber++; } // (8) Walker becomes its Ancestor, go to (3) LOG(2, "(8) Walker[" << Walker->getName() << "] is now its Ancestor " << Walker->Ancestor->getName() << ", backstepping. "); Walker = Walker->Ancestor; BackStepping = true; } } void DepthFirstSearchAnalysis::CleanRootStackDownTillWalker(atom *&Walker, bond::ptr &Binder, ConnectedSubgraph &Subgraph) { atom *OtherAtom = NULL; if (!BackStepping) { // coming from (8) want to go to (3) // (9) remove all from stack till Walker (including), these and Root form a component //AtomStack.Output(out); SetNextComponentNumber(Root, ComponentNumber); LOG(3, "(9) Root[" << Root->getName() << "]'s Component is " << ComponentNumber << "."); SetNextComponentNumber(Walker, ComponentNumber); LOG(3, "(9) Walker[" << Walker->getName() << "]'s Component is " << ComponentNumber << "."); do { ASSERT(!AtomStack.empty(), "DepthFirstSearchAnalysis::CleanRootStackDownTillWalker() - AtomStack is empty!"); OtherAtom = AtomStack.front(); AtomStack.pop_front(); Subgraph.push_back(OtherAtom); SetNextComponentNumber(OtherAtom, ComponentNumber); LOG(3, "(7) Other[" << OtherAtom->getName() << "]'s Component is " << ComponentNumber << "."); } while (OtherAtom != Walker); ComponentNumber++; // (11) Root is separation vertex, set Walker to Root and go to (4) Walker = Root; Binder = FindNextUnused(Walker); if (Binder != NULL) { // Root is separation vertex LOG(1, "(10) Walker is Root[" << Root->getName() << "], next Unused Bond is " << *Binder << "."); LOG(1, "(11) Root is a separation vertex."); Walker->SeparationVertex = true; } else { LOG(1, "(10) Walker is Root[" << Root->getName() << "], no next Unused Bond."); } } } const std::deque& DepthFirstSearchAnalysis::getBackEdgeStack() const { return BackEdgeStack; } void DepthFirstSearchAnalysis::operator()() { Info FunctionInfo("DepthFirstSearchAnalysis"); ListOfConnectedSubgraphs.clear(); int OldGraphNr = 0; atom *Walker = NULL; bond::ptr Binder; if (World::getInstance().numAtoms() == 0) return; Init(); LOG(0, "STATUS: Start walking the bond graph."); for(World::AtomIterator iter = World::getInstance().getAtomIter(); iter != World::getInstance().atomEnd();) { // don't advance, is done at the end Root = *iter; // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all AtomStack.clear(); // put into new subgraph molecule and add this to list of subgraphs ConnectedSubgraph CurrentSubgraph; CurrentSubgraph.push_back(Root); OldGraphNr = CurrentGraphNr; Walker = Root; do { // (10) do { // (2) set number and Lowpoint of Atom to i, increase i, push current atom SetWalkersGraphNr(Walker); ProbeAlongUnusedBond(Walker, Binder); if (Binder == NULL) { LOG(2, "No more Unused Bonds."); break; } else Binder.reset(); } while (1); // (2) // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished! if ((Walker == Root) && (Binder == NULL)) break; CheckForaNewComponent( Walker, CurrentSubgraph); CleanRootStackDownTillWalker(Walker, Binder, CurrentSubgraph); } while ((BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges ListOfConnectedSubgraphs.push_back(CurrentSubgraph); // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph std::stringstream output; output << CurrentSubgraph; LOG(1, "INFO: Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr-1 << ": " << output.str()); // step on to next root while (iter != World::getInstance().atomEnd()) { if ((*iter)->GraphNr != -1) { // if already discovered, step on iter++; } else { LOG(1,"Current next subgraph root candidate is " << (*iter)->getName() << " with GraphNr " << (*iter)->GraphNr << "."); break; } } } LOG(0, "STATUS: Done walking the bond graph."); // set cyclic bond criterium on "same LP" basis CyclicBondAnalysis(); OutputGraphInfoPerAtom(); OutputGraphInfoPerBond(); } void DepthFirstSearchAnalysis::UpdateMoleculeStructure() const { // remove all of World's molecules for (World::MoleculeIterator iter = World::getInstance().getMoleculeIter(); World::getInstance().getMoleculeIter() != World::getInstance().moleculeEnd(); iter = World::getInstance().getMoleculeIter()) { World::getInstance().destroyMolecule(*iter); } // instantiate new molecules molecule *newmol = NULL; for (ConnectedSubgraphList::const_iterator iter = ListOfConnectedSubgraphs.begin(); iter != ListOfConnectedSubgraphs.end(); ++iter) { newmol = (*iter).getMolecule(); if (DoLog(2)) { LOG(2, "STATUS: Creating new molecule:"); std::stringstream output; newmol->Output(&output); std::stringstream outstream(output.str()); std::string line; while (getline(outstream, line)) { LOG(2, "\t"+line); } } } } MoleculeLeafClass *DepthFirstSearchAnalysis::getMoleculeStructure() const { MoleculeLeafClass *Subgraphs = new MoleculeLeafClass(NULL); MoleculeLeafClass *MolecularWalker = Subgraphs; for (World::MoleculeIterator iter = World::getInstance().getMoleculeIter(); iter != World::getInstance().moleculeEnd(); ++iter) { // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include MolecularWalker = new MoleculeLeafClass(MolecularWalker); MolecularWalker->Leaf = (*iter); } return Subgraphs; }