source: src/Fragmentation/PowerSetGenerator.cpp@ bbff92

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include "PowerSetGenerator.hpp"

#include "CodePatterns/Log.hpp"

#include "atom.hpp"
#include "Bond/bond.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/UniqueFragments.hpp"

/** Constructor of class PowerSetGenerator.
 *
 */
PowerSetGenerator::PowerSetGenerator(class UniqueFragments *_FragmentSearch, int Order) :
  BondsPerSPList(Order),
  FragmentSearch(_FragmentSearch)
{}

/** Destructor of class PowerSetGenerator.
 *
 */
PowerSetGenerator::~PowerSetGenerator()
{
  FragmentSearch = NULL;
}

/** Clears the touched list
 * \param verbosity verbosity level
 * \param *&TouchedList touched list
 * \param SubOrder current suborder
 * \param TouchedIndex currently touched
 */
void PowerSetGenerator::ClearingTouched(int verbosity, int *&TouchedList, int SubOrder, int &TouchedIndex)
{
  Log() << Verbose(1+verbosity) << "Clearing touched list." << endl;
  for (TouchedIndex=SubOrder+1;TouchedIndex--;)  // empty touched list
    TouchedList[TouchedIndex] = -1;
  TouchedIndex = 0;

}

/** Adds the current combination of the power set to the snake stack.
 * \param verbosity verbosity level
 * \param CurrentCombination
 * \param SetDimension maximum number of bits in power set
 * \param *FragmentSet snake stack to remove from
 * \param &BondsSet set of bonds
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 * \return number of set bits
 */
int PowerSetGenerator::AddPowersetToSnakeStack(int verbosity, int CurrentCombination, int SetDimension, KeySet *FragmentSet, std::vector<bond *> &BondsSet, int *&TouchedList, int &TouchedIndex)
{
  atom *OtherWalker = NULL;
  bool bit = false;
  KeySetTestPair TestKeySetInsert;

  int Added = 0;
  for (int j=0;j<SetDimension;j++) {  // pull out every bit by shifting
    bit = ((CurrentCombination & (1 << j)) != 0);  // mask the bit for the j-th bond
    if (bit) {  // if bit is set, we add this bond partner
      OtherWalker = BondsSet[j]->rightatom;  // rightatom is always the one more distant, i.e. the one to add
      //Log() << Verbose(1+verbosity) << "Current Bond is " << BondsSet[j] << ", checking on " << *OtherWalker << "." << endl;
      Log() << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->getNr() << "." << endl;
      TestKeySetInsert = FragmentSet->insert(OtherWalker->getNr());
      if (TestKeySetInsert.second) {
        TouchedList[TouchedIndex++] = OtherWalker->getNr();  // note as added
        Added++;
      } else {
        Log() << Verbose(2+verbosity) << "This was item was already present in the keyset." << endl;
      }
    } else {
      Log() << Verbose(2+verbosity) << "Not adding." << endl;
    }
  }
  return Added;
};

/** Counts the number of elements in a power set.
 * \param SetFirst begin iterator first bond
 * \param SetLast end iterator
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 * \return number of elements
 */
int PowerSetGenerator::CountSetMembers(std::list<bond *>::const_iterator SetFirst, std::list<bond *>::const_iterator SetLast, int *&TouchedList, int TouchedIndex)
{
  int SetDimension = 0;
  for( std::list<bond *>::const_iterator Binder = SetFirst;
      Binder != SetLast;
      ++Binder) {
    for (int k=TouchedIndex;k--;) {
      if ((*Binder)->Contains(TouchedList[k]))   // if we added this very endpiece
        SetDimension++;
    }
  }
  return SetDimension;
};

/** Fills a list of bonds from another
 * \param *BondsList bonds array/vector to fill
 * \param SetFirst begin iterator first bond
 * \param SetLast end iterator
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 * \return number of elements
 */
int PowerSetGenerator::FillBondsList(std::vector<bond *> &BondsList, std::list<bond *>::const_iterator SetFirst, std::list<bond *>::const_iterator SetLast, int *&TouchedList, int TouchedIndex)
{
  int SetDimension = 0;
  for( std::list<bond *>::const_iterator Binder = SetFirst;
      Binder != SetLast;
      ++Binder) {
    for (int k=0;k<TouchedIndex;k++) {
      if ((*Binder)->leftatom->getNr() == TouchedList[k])   // leftatom is always the closer one
        BondsList[SetDimension++] = (*Binder);
    }
  }
  return SetDimension;
};

/** Remove all items that were added on this SP level.
 * \param verbosity verbosity level
 * \param *FragmentSet snake stack to remove from
 * \param *&TouchedList touched list
 * \param TouchedIndex currently touched
 */
void PowerSetGenerator::RemoveAllTouchedFromSnakeStack(int verbosity, KeySet *FragmentSet, int *&TouchedList, int &TouchedIndex)
{
  int Removal = 0;
  for(int j=0;j<TouchedIndex;j++) {
    Removal = TouchedList[j];
    Log() << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
    FragmentSet->erase(Removal);
    TouchedList[j] = -1;
  }
  DoLog(2) && (Log() << Verbose(2) << "Remaining local nr.s on snake stack are: ");
  for(KeySet::iterator runner = FragmentSet->begin(); runner != FragmentSet->end(); runner++)
    DoLog(0) && (Log() << Verbose(0) << (*runner) << " ");
  DoLog(0) && (Log() << Verbose(0) << endl);
  TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
};


/** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment for a given root vertex in the context of \a this molecule.
 * -# initialises UniqueFragments structure
 * -# fills edge list via BFS
 * -# creates the fragment by calling recursive function SPFragmentGenerator with UniqueFragments structure, 0 as
 root distance, the edge set, its dimension and the current suborder
 * -# Free'ing structure
 * Note that we may use the fact that the atoms are SP-ordered on the atomstack. I.e. when popping always the last, we first get all
 * with SP of 2, then those with SP of 3, then those with SP of 4 and so on.
 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
 * \param saturation whether to treat hydrogen special or not
 * \return number of inserted fragments
 * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
 */
int PowerSetGenerator::operator()(KeySet &RestrictedKeySet, const enum HydrogenSaturation saturation)
{
  int Counter = FragmentSearch->FragmentCounter; // mark current value of counter

  DoLog(0) && (Log() << Verbose(0) << endl);
  DoLog(0) && (Log() << Verbose(0) << "Begin of PowerSetGenerator with order " << BondsPerSPList.getOrder() << " at Root " << *FragmentSearch->getRoot() << "." << endl);

  BondsPerSPList.SetSPList(FragmentSearch->getRoot());

  // do a BFS search to fill the SP lists and label the found vertices
  BondsPerSPList.FillSPListandLabelVertices(FragmentSearch->getRoot()->GetTrueFather()->getNr(), RestrictedKeySet, saturation);

  // outputting all list for debugging
  BondsPerSPList.OutputSPList();

  // creating fragments with the found edge sets  (may be done in reverse order, faster)
  int SP = BondsPerSPList.CountNumbersInBondsList();
  DoLog(0) && (Log() << Verbose(0) << "Total number of edges is " << SP << "." << endl);
  if (SP >= (BondsPerSPList.getOrder()-1)) {
    // start with root (push on fragment stack)
    DoLog(0) && (Log() << Verbose(0) << "Starting fragment generation with " << *FragmentSearch->getRoot() << ", local nr is " << FragmentSearch->getRoot()->getNr() << "." << endl);
    FragmentSearch->FragmentSet->clear();
    DoLog(0) && (Log() << Verbose(0) << "Preparing subset for this root and calling generator." << endl);

    // prepare the subset and call the generator
    std::vector<bond*> BondsList;
    BondsList.resize(BondsPerSPList.BondsPerSPCount[0]);
    ASSERT(BondsPerSPList.BondsPerSPList[0].size() != 0,
        "molecule::PowerSetGenerator() - BondsPerSPList.BondsPerSPList[0] contains no root bond.");
    BondsList[0] = (*BondsPerSPList.BondsPerSPList[0].begin());  // on SP level 0 there's only the root bond

    SPFragmentGenerator(0, BondsList, BondsPerSPList.BondsPerSPCount[0], BondsPerSPList.getOrder());
  } else {
    DoLog(0) && (Log() << Verbose(0) << "Not enough total number of edges to build " << BondsPerSPList.getOrder() << "-body fragments." << endl);
  }

  // as FragmentSearch structure is used only once, we don't have to clean it anymore
  // remove root from stack
  DoLog(0) && (Log() << Verbose(0) << "Removing root again from stack." << endl);
  FragmentSearch->FragmentSet->erase(FragmentSearch->getRoot()->getNr());

  // free'ing the bonds lists
  BondsPerSPList.ResetSPList();

  // return list
  DoLog(0) && (Log() << Verbose(0) << "End of PowerSetGenerator." << endl);
  return (FragmentSearch->FragmentCounter - Counter);
};

/** From a given set of Bond sorted by Shortest Path distance, create all possible fragments of size \a SetDimension.
 * -# loops over every possible combination (2^dimension of edge set)
 *  -# inserts current set, if there's still space left
 *    -# yes: calls SPFragmentGenerator with structure, created new edge list and size respective to root dist
ance+1
 *    -# no: stores fragment into keyset list by calling InsertFragmentIntoGraph
 *  -# removes all items added into the snake stack (in UniqueFragments structure) added during level (root
distance) and current set
 * \param RootDistance current shortest path level, whose set of edges is represented by **BondsSet
 * \param BondsSet array of bonds to check
 * \param SetDimension Number of possible bonds on this level (i.e. size of the array BondsSet[])
 * \param SubOrder remaining number of allowed vertices to add
 */
void PowerSetGenerator::SPFragmentGenerator(int RootDistance, std::vector<bond *> &BondsSet, int SetDimension, int SubOrder)
{
  int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
  int NumCombinations;
  int bits, TouchedIndex, SubSetDimension, SP, Added;
  int SpaceLeft;
  int *TouchedList = new int[SubOrder + 1];
  KeySetTestPair TestKeySetInsert;

  NumCombinations = 1 << SetDimension;

  // here for all bonds of Walker all combinations of end pieces (from the bonds)
  // have to be added and for the remaining ANOVA order GraphCrawler be called
  // recursively for the next level

  Log() << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
  Log() << Verbose(1+verbosity) << "We are " << RootDistance << " away from Root, which is " << *FragmentSearch->getRoot() << ", SubOrder is " << SubOrder << ", SetDimension is " << SetDimension << " and this means " <<  NumCombinations-1 << " combination(s)." << endl;

  // initialised touched list (stores added atoms on this level)
  ClearingTouched(verbosity, TouchedList, SubOrder, TouchedIndex);

  // create every possible combination of the endpieces
  Log() << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
  for (int i=1;i<NumCombinations;i++) {  // sweep through all power set combinations (skip empty set!)
    // count the set bit of i
    bits = 0;
    for (int j=SetDimension;j--;)
      bits += (i & (1 << j)) >> j;

    Log() << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
    if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
      // --1-- add this set of the power set of bond partners to the snake stack
      Added = AddPowersetToSnakeStack(verbosity, i, SetDimension, FragmentSearch->FragmentSet, BondsSet, TouchedList, TouchedIndex);

      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
      if (SpaceLeft > 0) {
        Log() << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
        if (SubOrder > 1) {    // Due to Added above we have to check extra whether we're not already reaching beyond the desired Order
          // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
          SP = RootDistance+1;  // this is the next level

          // first count the members in the subset
          SubSetDimension = CountSetMembers(BondsPerSPList.BondsPerSPList[SP].begin(), BondsPerSPList.BondsPerSPList[SP].end(), TouchedList, TouchedIndex);

          // then allocate and fill the list
          std::vector<bond *> BondsList;
          BondsList.resize(SubSetDimension);
          SubSetDimension = FillBondsList(BondsList, BondsPerSPList.BondsPerSPList[SP].begin(), BondsPerSPList.BondsPerSPList[SP].end(), TouchedList, TouchedIndex);

          // then iterate
          Log() << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->getRoot() << " with sub set dimension " << SubSetDimension << "." << endl;
          SPFragmentGenerator(SP, BondsList, SubSetDimension, SubOrder-bits);
        }
      } else {
        // --2-- otherwise store the complete fragment
        Log() << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
        // store fragment as a KeySet
        DoLog(2) && (Log() << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ");
        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
          DoLog(0) && (Log() << Verbose(0) << (*runner) << " ");
        DoLog(0) && (Log() << Verbose(0) << endl);
        FragmentSearch->InsertFragmentIntoGraph();
      }

      // --3-- remove all added items in this level from snake stack
      Log() << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
      RemoveAllTouchedFromSnakeStack(verbosity, FragmentSearch->FragmentSet, TouchedList, TouchedIndex);
    } else {
      Log() << Verbose(2+verbosity) << "More atoms to add for this set (" << bits << ") than space left on stack " << SubOrder << ", skipping this set." << endl;
    }
  }
  delete[](TouchedList);
  Log() << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->getRoot() << " and SubOrder is " << SubOrder << "." << endl;
};