source: src/Graph/CyclicStructureAnalysis.cpp@ 0c72e2

/*
 * 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 a cycle.
 * \param OtherAtom pointing to Root on return indicating found cycle
 * \param *&BackEdge the edge from root that we don't want to move along
 */
void CyclicStructureAnalysis::CyclicBFSFromRootToRoot(atom *&OtherAtom, bond::ptr &BackEdge)
{
  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 << ".");
    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) {
            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()]) << ".");
    if (OtherAtom == Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
      // step through predecessor list
      LOG(4, "DEBUG: Checking whether all predecessors are already marked cyclic ...");
      while (OtherAtom != BackEdge->rightatom) {  // Note that leftatom is Root itself
        if (!OtherAtom->GetTrueFather()->IsCyclic) { // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
          LOG(4, "\tDEBUG: OtherAtom " << *OtherAtom << " is not cyclic, breaking.");
          break;
        } else
          OtherAtom = PredecessorList[OtherAtom->getNr()];
      }
      LOG(4, "DEBUG: Checking done.");
      // if each atom in found cycle is cyclic, loop's been found before already
      if (OtherAtom == BackEdge->rightatom) { // loop got round completely
        LOG(3, "INFO: All bonds are marked cyclic already, checking allcycles whether cycle is already present.");
        const cycle_t currentcycle = extractCurrentCycle(BackEdge);
        const cycles_t::const_iterator iter =
            std::find(allcycles.begin(), allcycles.end(), currentcycle);
        if (iter == allcycles.end()) { // not found
          OtherAtom = Root;
          LOG(2, "INFO: Cycle is not present.");
          LOG(2, "INFO: We have reached Root " << *OtherAtom << " and may extract the cycle.");
        } else {
          LOG(2, "INFO: Cycle is already present.");
          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
        }
      } else {
        OtherAtom = Root;
        LOG(2, "INFO: We have reached Root " << *OtherAtom << " and may extract the cycle.");
      }
    }
  } 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
 * \param &NumCyles number of cycles in graph
 */
void CyclicStructureAnalysis::RetrieveCycleMembers(
    atom *&OtherAtom,
    bond::ptr &BackEdge,
    int &MinRingSize,
    int &NumCycles)
{
  atom *Walker = NULL;
  int RingSize = -1;

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

    // exctract cycle
    {
      allcycles.push_back(extractCurrentCycle(BackEdge));
      CyclicStructureAnalysis::cycle_t &lastcycle = allcycles.back();
      RingSize = lastcycle.size();
      LOG(0, "INFO: " << "Found ring contains: " << lastcycle << "  with a length of " << RingSize);
    }

    // walk through all and set MinimumRingSize
    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.");
  }
}

/** 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
 * \param NumCyles number of cycles in graph
 */
void CyclicStructureAnalysis::AssignRingSizetoNonCycleMembers(const int MinRingSize, const int NumCycles)
{
  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 " << NumCycles << " 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 NumCycles = 0;
  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 ... ");
  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, "---------------------------------------------------------------------------------------------------------");
    OtherAtom = NULL;
    // go to next cycle via BFS
    CyclicBFSFromRootToRoot(OtherAtom, BackEdge);
    // get all member nodes of this cycle
    RetrieveCycleMembers(OtherAtom, BackEdge, MinRingSize, NumCycles);

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