source: src/Fragmentation/BondsPerShortestPath.cpp@ 20bb3b

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

/*
 * BondsPerShortestPath.cpp
 *
 *  Created on: Oct 18, 2011
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "BondsPerShortestPath.hpp"

#include <sstream>

#include "CodePatterns/Log.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Element/element.hpp"
#include "Fragmentation/KeySet.hpp"

BondsPerShortestPath::BondsPerShortestPath(int _Order) :
  Order(_Order)
{
  InitialiseSPList();
}

BondsPerShortestPath::~BondsPerShortestPath()
{
  // free Order-dependent entries of UniqueFragments structure for next loop cycle
  FreeSPList();
}

/** Allocates memory for BondsPerShortestPath::BondsPerSPList.
 * \sa BondsPerShortestPath::FreeSPList()
 */
void BondsPerShortestPath::InitialiseSPList()
{
  BondsPerSPList.resize(Order);
  BondsPerSPCount = new int[Order];
  for (int i=Order;i--;) {
    BondsPerSPCount[i] = 0;
  }
};

/** Free's memory for for BondsPerShortestPath::BondsPerSPList.
 * \sa BondsPerShortestPath::InitialiseSPList()
 */
void BondsPerShortestPath::FreeSPList()
{
  delete[](BondsPerSPCount);
};

/** Sets FragmenSearch to initial value.
 * Sets BondsPerShortestPath::ShortestPathList entries to zero, BondsPerShortestPath::BondsPerSPCount to zero (except zero level to 1) and
 * adds initial bond BondsPerShortestPath::Root to BondsPerShortestPath::Root to BondsPerShortestPath::BondsPerSPList
 * \param *_Root root node, self loop becomes first bond
 * \sa BondsPerShortestPath::FreeSPList()
 */
void BondsPerShortestPath::SetSPList(atom *_Root)
{
  // prepare root level (SP = 0) and a loop bond denoting Root
  for (int i=Order;i--;)
    BondsPerSPCount[i] = 0;
  BondsPerSPCount[0] = 1;
  bond::ptr Binder(new bond(_Root, _Root));
  BondsPerSPList[0].push_back(Binder);
};

/** Resets BondsPerShortestPath::ShortestPathList and cleans bonds from BondsPerShortestPath::BondsPerSPList.
 * \sa BondsPerShortestPath::InitialiseSPList()
 */
void BondsPerShortestPath::ResetSPList()
{
  LOG(1, "INFO: Free'ing all found lists and resetting index lists");
  for(int i=Order;i--;) {
    std::stringstream output;
    output << "Current SP level is " << i << ": ";
    // delete added bonds
    for (BondsPerSP::iterator iter = BondsPerSPList[i].begin();
        iter != BondsPerSPList[i].end();
        ++iter) {
      // TODO: Hack because we have not registered bond's in BondsPerSPList with atoms
      (*iter)->leftatom = NULL;
      (*iter)->rightatom = NULL;
    }
    BondsPerSPList[i].clear();
    // also start and end node
    output << "cleaned.";
    LOG(1, output.str());
  }
};


/** Fills the Bonds per Shortest Path List and set the vertex labels.
 * \param _RootKeyNr index of root node
 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
 * \param treatment this tells whether to treat hydrogen special or not.
 */
void BondsPerShortestPath::FillSPListandLabelVertices(int _RootKeyNr, KeySet &RestrictedKeySet, const enum HydrogenTreatment treatment)
{
  // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
  // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
  // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
  // (EdgeinSPLevel) of this tree ...
  // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
  // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
//  int AtomKeyNr = -1;
  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  atom *Predecessor = NULL;
  bond::ptr Binder;
  int RootKeyNr = _RootKeyNr;
  int RemainingWalkers = -1;
  int SP = -1;

  LOG(0, "INFO: Starting BFS analysis ...");
  for (SP = 0; SP < (Order-1); SP++) {
    {
      std::stringstream output;
      output << "New SP level reached: " << SP << ", creating new SP list with " << BondsPerSPCount[SP] << " item(s)";
      if (SP > 0) {
        output << ", old level closed with " << BondsPerSPCount[SP-1] << " item(s).";
        BondsPerSPCount[SP] = 0;
      } else
        output << ".";
      LOG(1, output.str());
    }

    RemainingWalkers = BondsPerSPCount[SP];
    for (BondsPerSP::const_iterator CurrentEdge = BondsPerSPList[SP].begin();
        CurrentEdge != BondsPerSPList[SP].end();
        ++CurrentEdge) { /// start till end of this SP level's list
      RemainingWalkers--;
      Walker = (*CurrentEdge)->rightatom;    // rightatom is always the one more distant
      Predecessor = (*CurrentEdge)->leftatom;    // ... and leftatom is predecessor
//      AtomKeyNr = Walker->getNr();
      LOG(1, "INFO: Current Walker is: " << *Walker << " with nr " << Walker->getNr() << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level.");
      // check for new sp level
      // go through all its bonds
//      LOG(2, "Going through all bonds of Walker.");
      const BondList& ListOfBonds = Walker->getListOfBonds();
      for (BondList::const_iterator Runner = ListOfBonds.begin();
          Runner != ListOfBonds.end();
          ++Runner) {
        OtherWalker = (*Runner)->GetOtherAtom(Walker);
        if ((RestrictedKeySet.find(OtherWalker->getNr()) != RestrictedKeySet.end())
            // skip hydrogens if desired and restrict to fragment
            && ((treatment == IncludeHydrogen) || (OtherWalker->getType()->getAtomicNumber() != 1))) {
          LOG(2, "Current partner is " << *OtherWalker << " with nr " << OtherWalker->getNr() << " in bond " << *(*Runner) << ".");
          // set the label if not set (and push on root stack as well)
          if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->getNr() > RootKeyNr)) { // only pass through those with label bigger than Root's
            // add the bond in between to the SP list
            Binder.reset(new bond(Walker, OtherWalker)); // create a new bond in such a manner, that bond::rightatom is always the one more distant
            BondsPerSPList[SP+1].push_back(Binder);
            BondsPerSPCount[SP+1]++;
            LOG(3, "DEBUG: Added its bond to SP list, having now " << BondsPerSPCount[SP+1] << " item(s).");
          } else {
            if (OtherWalker != Predecessor)
              LOG(3, "DEBUG: Not passing on, as index of " << *OtherWalker << " " << OtherWalker->GetTrueFather()->getNr() << " is smaller than that of Root " << RootKeyNr << ".");
            else
              LOG(3, "DEBUG: This is my predecessor " << *Predecessor << ".");
          }
        } else LOG(2, "DEBUG: Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << ".");
      }
    }
  }
};

/** prints the Bonds per Shortest Path list in BondsPerShortestPath.
 */
void BondsPerShortestPath::OutputSPList()
{
  LOG(2, "DEBUG: Printing all found lists.");
  for(int i=1;i<Order;i++) {    // skip the root edge in the printing
    LOG(2, "\tCurrent SP level is " << i << ".");
    std::stringstream output;
    for (BondsPerShortestPath::BondsPerSP::const_iterator Binder = BondsPerSPList[i].begin();
        Binder != BondsPerSPList[i].end();
        ++Binder) {
      output << "\t" << *Binder;
    }
    output << std::endl;
    LOG(2, output.str());
  }
};

/** Simply counts all bonds in all BondsPerShortestPath::BondsPerSPList lists.
 */
int BondsPerShortestPath::CountNumbersInBondsList()
{
  int SP = -1;  // the Root <-> Root edge must be subtracted!
  for(int i=Order;i--;) { // sum up all found edges
    for (BondsPerShortestPath::BondsPerSP::const_iterator Binder = BondsPerSPList[i].begin();
        Binder != BondsPerSPList[i].end();
        ++Binder) {
      SP++;
    }
  }
  return SP;
};