source: src/Graph/CyclicStructureAnalysis.cpp@ cceb8c

/*
 * 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 <algorithm>

#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 HydrogenTreatment _treatment) :
  treatment(_treatment)
{}

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)
{
  ColorList.clear();
  ShortestPathList.clear();
  ShortestPathList[Root->getNr()] = 0;
  PredecessorList.clear();
  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 all cycles.
 *
 * We exclude the back edge, hence the direct way is prohibited. But as it is a back edge,
 * there must be at least one more way to \a *Root. This leads us to all cycles for this
 * back edge.
 *
 * \param OtherAtom pointing to Root on return indicating found cycle
 * \param *&BackEdge the edge from root that we don't want to move along
 * \param MinRingSize set to minimum over all cycles encountered
 */
void CyclicStructureAnalysis::findAllCyclesforBackEdge(
    atom *&OtherAtom,
    bond::ptr &BackEdge,
    int &MinRingSize)
{
  size_t MaximumHorizon = (size_t)-1; // will overflow to largest number
  atom *Walker = NULL;
  do { // look for Root
    ASSERT(!BFSStack.empty(), "CyclicStructureAnalysis_CyclicBFSFromRootToRoot() - BFSStack is empty!");
    Walker = BFSStack.back();
    BFSStack.pop_back();
    LOG(2, "INFO: Current Walker is " << *Walker << ", we look for SP to Root " << *Root << ".");

    // check all edges/bonds from current Walker
    if (MaximumHorizon >=  (size_t)ShortestPathList[Walker->getNr()]) {
      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 ((treatment == IncludeHydrogen) || (OtherAtom->getType()->getAtomicNumber() != 1)) {
            LOG(2, "INFO: Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << ".");
            if (ColorList[OtherAtom->getNr()] == GraphEdge::white) {
              // we discovered a new node/atom
              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.");
              //if (ShortestPathList[OtherAtom->getNr()] < MinimumRingSize[Walker->GetTrueFather()->getNr()]) { // Check for maximum distance
              LOG(3, "ACCEPT: Putting OtherAtom " << OtherAtom->getName() << " into queue.");
              BFSStack.push_front(OtherAtom);
              //}
            } else {
              LOG(3, "REJECT: Not Adding, has already been visited.");
            }
            if (OtherAtom == Root)
              break;
          } else {
            LOG(2, "INFO: Skipping hydrogen atom " << *OtherAtom << ".");
            ColorList[OtherAtom->getNr()] = GraphEdge::black;
          }
        } else {
          LOG(2, "REJECT: Bond " << *(*Runner) << " not Visiting, is the back edge.");
        }
      }
      ColorList[Walker->getNr()] = GraphEdge::black;
      LOG(1, "INFO: Coloring Walker " << Walker->getName() << " " << GraphEdge::getColorName(ColorList[Walker->getNr()]) << ".");
    } else {
      LOG(1, "INFO: Skipping bonds for " << Walker->getName() << " as it resides on the horizon.");
    }

    // have we closed the cycle, i.e. stumbled upon Root by another mean than
    // the back edge?
    if (OtherAtom == Root) {
      // check whether this cycle wasn't already found beforehand by stepping
      // through predecessor list
      int RingSize = RetrieveCycleMembers(OtherAtom, BackEdge, MinRingSize);
      MaximumHorizon = std::min( MaximumHorizon, (size_t)RingSize );
      LOG(2, "INFO: Maximum horizon is set to " << MaximumHorizon);

      // remove last step and prepare for a possible yet another path to Root
      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.");
          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
    }
  } while ((!BFSStack.empty()) && (OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->getNr()] < MinimumRingSize[Walker->GetTrueFather()->getNr()])));
}

/** Extracts cycle from BFSAccounting::PredecessorList with given \a BackEdge and current \a Root.
 *
 * \param BackEdge back edge of this cycle
 */
CyclicStructureAnalysis::cycle_t CyclicStructureAnalysis::extractCurrentCycle(
    bond::ptr &BackEdge)
{
  CyclicStructureAnalysis::cycle_t currentcycle;
  atom *Walker = Root;
  currentcycle.insert(Walker->GetTrueFather()->getId());
  std::stringstream output;
  while (Walker != BackEdge->rightatom) { // leftatom is root
    Walker = PredecessorList[Walker->getNr()];
    Walker->GetTrueFather()->IsCyclic = true;
    output << Walker->getName() << " <-> ";
#ifndef NDEBUG
    std::pair< cycle_t::iterator, bool > inserter =
#endif
        currentcycle.insert(Walker->GetTrueFather()->getId());
    ASSERT( inserter.second,
        "CyclicStructureAnalysis::RetrieveCycleMembers() - we already inserted "
        +toString(Walker->GetTrueFather()->getId())+" into currentcycle.");
  }
  LOG(3, "INFO: " << output.str() << Root->getName());
  return currentcycle;
}

/** 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
 * \return size of found cycle, -1 - nothing found, something went wrong
 */
int CyclicStructureAnalysis::RetrieveCycleMembers(
    atom *&OtherAtom,
    bond::ptr &BackEdge,
    int &MinRingSize)
{
  int RingSize = -1;
  bool newcycle = false;

  if (OtherAtom == Root) {
    // now climb back the predecessor list and thus find the cycle members
    Root->GetTrueFather()->IsCyclic = true;

    {
      // check whether cycle is present already
      atom *Walker = Root;
      LOG(4, "DEBUG: Checking whether all predecessors are already marked cyclic ...");
      while (Walker != BackEdge->rightatom) {  // Note that leftatom is Root itself
        if (!Walker->GetTrueFather()->IsCyclic) { // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
          LOG(4, "\tDEBUG: Walker " << *Walker << " is not cyclic, breaking.");
          break;
        } else
          Walker = PredecessorList[Walker->getNr()];
      }
      LOG(4, "DEBUG: Checking done.");

      // if each atom in found cycle is cyclic, loop's been found before already

      // exctract cycle
      {
        const cycle_t currentcycle = extractCurrentCycle(BackEdge);
        if (Walker != BackEdge->rightatom) { // loop got round completely
          allcycles.push_back(currentcycle);
          newcycle = true;
        } else {
          LOG(3, "INFO: All bonds are marked cyclic already, checking allcycles whether cycle is already present.");
          const cycles_t::const_iterator iter =
              std::find(allcycles.begin(), allcycles.end(), currentcycle);
          if (iter == allcycles.end()) { // not found
            allcycles.push_back(currentcycle);
            newcycle = true;
          }
        }
        RingSize = currentcycle.size();
        if (newcycle) {
          LOG(0, "INFO: " << "Found ring contains: " << currentcycle << "  with a length of " << RingSize);
        } else {
          LOG(1, "INFO: cycle " << currentcycle << " is already present.");
        }
      }
    }

    if (newcycle) {
      // walk through all and set MinimumRingSize
      atom *Walker = Root;
      if ((MinimumRingSize.count(Walker->GetTrueFather()->getNr()) == 0)
          || (RingSize < MinimumRingSize[Walker->GetTrueFather()->getNr()])) {
        MinimumRingSize[Walker->GetTrueFather()->getNr()] = RingSize;
      } else {
        LOG(3, "INFO: Not setting MinimumRingSize of "<< *(Walker->GetTrueFather())
            << " to " << RingSize << " which is already set to "
            << MinimumRingSize[Walker->GetTrueFather()->getNr()] << ".");
      }
      while (Walker != BackEdge->rightatom) { // note that Root is leftatom
        Walker = PredecessorList[Walker->getNr()];
        if ((MinimumRingSize.count(Walker->GetTrueFather()->getNr()) == 0)
            || (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.");
  }
  return RingSize;
}

/** From a given node performs a BFS to touch the next cycle, for whose nodes \a MinimumRingSize is set and set it accordingly.
 * \param *&Walker 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 *Walker)
{
  atom *Root = 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 << ", BFS-stepping away from 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

      // 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
//      if ((ListOfBonds.size() != 1)) {
        OtherAtom = (*Runner)->GetOtherAtom(Walker);
        if ((treatment == IncludeHydrogen) || (OtherAtom->getType()->getAtomicNumber() != 1)) {
          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.");
            // distance is a locally optimal criterion (we have eliminated all
            // cycles already). Hence, we may assume that all set MinimumRingSize
            // correspond to shortest distances to cycles. I.e., as soon as we reach
            // as set MinimumRingSize we may use it and the current shortest path
            // distance to it
            if (MinimumRingSize.count(OtherAtom->GetTrueFather()->getNr())) {
              LOG(2, "SUCCESS: Found set MinimumRingSize at " << *OtherAtom
                  << ", walking back to Root " << *Root << ".");
              // set all predecessors
              const unsigned int shorttestpath = ShortestPathList[OtherAtom->getNr()];
              atom *Backwalker = OtherAtom;
              while (Backwalker != Root) {
                Backwalker = PredecessorList[Backwalker->getNr()];
                MinimumRingSize[Backwalker->GetTrueFather()->getNr()] =
                    (shorttestpath - ShortestPathList[Backwalker->getNr()])
                    + MinimumRingSize[OtherAtom->GetTrueFather()->getNr()];
                LOG(2, "Setting MinimumRingSize of " << *Backwalker << " to "
                    << MinimumRingSize[Backwalker->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 to hydrogen.");
        }
//      }
    }
    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 *&MinimumRingSize by BFS to next cycle.
 * \param *&MinimumRingSize array with minimum distance without encountering oneself for each atom
 * \param MinRingSize global minium distance
 */
void CyclicStructureAnalysis::AssignRingSizetoNonCycleMembers(const int MinRingSize)
{
  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) {
      Walker = *iter;

      if (MinimumRingSize.find(Walker->GetTrueFather()->getNr()) == MinimumRingSize.end()) { // check whether MinimumRingSize is set, if not BFS to next where it is
        LOG(1, "---------------------------------------------------------------------------------------------------------");
        BFSToNextCycle(Walker);
      }
      ASSERT(MinimumRingSize.find(Walker->GetTrueFather()->getNr()) != MinimumRingSize.end(),
          "CyclicStructureAnalysis::AssignRingSizetoNonCycleMembers() - BFSToNextCycle did not set MinimumRingSize of "
          +toString(*(Walker->GetTrueFather()))+".");
      LOG(1, "INFO: Minimum ring size of " << *Walker << " is " << MinimumRingSize[Walker->GetTrueFather()->getNr()] << ".");
    }
    LOG(1, "INFO: Minimum ring size is " << MinRingSize << ", over " << allcycles.size() << " cycle(s) total.");
  } else
    LOG(1, "INFO: No rings were detected in the molecular structure.");
}

/** 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 MinRingSize = -1;

  // clear cycle container
  allcycles.clear();

  {
    std::stringstream output;
    output << "Back edge list - ";
    for (std::deque<bond::ptr >::const_iterator iter = BackEdgeStack->begin();
        iter != BackEdgeStack->end(); ++iter)
      output << **iter << " ";
    LOG(0, output.str());
  }

  LOG(1, "STATUS: Analysing cycles ... ");
  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, "---------------------------------------------------------------------------------------------------------");
    OtherAtom = NULL;
    // find all cycles for current BackEdge
    findAllCyclesforBackEdge(OtherAtom, BackEdge, MinRingSize);

    CleanAllTouched();
  }
  AssignRingSizetoNonCycleMembers(MinRingSize);
}

/** 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;
}