source: src/Graph/DepthFirstSearchAnalysis.cpp@ df5b8c

/*
 * 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 <http://www.gnu.org/licenses/>.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include "DepthFirstSearchAnalysis.hpp"

#include <algorithm>
#include <functional>

#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 "MoleculeListClass.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<bond::ptr > *&LocalStack) const
{
  bool status = true;
  if (BackEdgeStack.empty()) {
    ELOG(1, "Reference BackEdgeStack is empty!");
    return false;
  }
  std::deque<bond::ptr > 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<const World &>(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<const World &>(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<const World &>(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<bond::ptr >& 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()) {
    // TODO: remove when deprecated MoleculeListClass is gone
    World::getInstance().getMolecules()->erase(*iter);
    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;
}