source: src/Graph/CyclicStructureAnalysis.cpp@ f68c68

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. 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 <http://www.gnu.org/licenses/>.
 */

/*
 * CyclicStructureAnalysis.cpp
 *
 *  Created on: Feb 16, 2011
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "CyclicStructureAnalysis.hpp"

#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 "Element/element.hpp"
#include "molecule.hpp"

CyclicStructureAnalysis::CyclicStructureAnalysis(const enum HydrogenSaturation _saturation) :
  saturation(_saturation)
{}

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::ptr &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);
        if ((saturation == DontSaturate) || (OtherAtom->getType()->getAtomicNumber() != 1)) {
          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;
        } else {
          LOG(2, "INFO: Skipping hydrogen atom " << *OtherAtom << "." << endl);
          ColorList[OtherAtom->getNr()] = GraphEdge::black;
        }
      } 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<atom *>::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::ptr &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<bond::ptr > * BackEdgeStack)
{
  Info FunctionInfo("CyclicStructureAnalysis");
  atom *Walker = NULL;
  atom *OtherAtom = NULL;
  bond::ptr BackEdge;
  int NumCycles = 0;
  int MinRingSize = -1;

  //std::stringstream output;
  //output << "Back edge list - ";
  //BackEdgeStack->Output(output);
  //LOG(0, output.str());

  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<atomId_t, int >& CyclicStructureAnalysis::getMinimumRingSize() const
{
  return MinimumRingSize;
}