source: src/Fragmentation/PowerSetGenerator.cpp@ 4464ef

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

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

#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"

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

/** Constructor of class PowerSetGenerator.
 *
 */
PowerSetGenerator::PowerSetGenerator(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, TouchedList_t &TouchedList, int SubOrder, int &TouchedIndex)
{
  LOG(1+verbosity, "Clearing touched list.");
  TouchedList.resize(SubOrder+1, -1);
  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::ptr > &BondsSet, TouchedList_t &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(1+verbosity, "Current Bond is " << BondsSet[j] << ", checking on " << *OtherWalker << ".");
      LOG(2+verbosity, "Adding " << *OtherWalker << " with nr " << OtherWalker->getNr() << ".");
      TestKeySetInsert = FragmentSet->insert(OtherWalker->getNr());
      if (TestKeySetInsert.second) {
        TouchedList[TouchedIndex++] = OtherWalker->getNr();  // note as added
        Added++;
      } else {
        LOG(2+verbosity, "This was item was already present in the keyset.");
      }
    } else {
      LOG(2+verbosity, "Not adding.");
    }
  }
  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::ptr >::const_iterator SetFirst, std::list<bond::ptr >::const_iterator SetLast, const TouchedList_t &TouchedList, int TouchedIndex)
{
  int SetDimension = 0;
  for( std::list<bond::ptr >::const_iterator Binder = SetFirst;
      Binder != SetLast;
      ++Binder) {
    for (TouchedList_t::const_iterator iter = TouchedList.begin();
        iter != TouchedList.end(); ++iter) {
      if ((*Binder)->ContainsNr(*iter))   // 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::ptr > &BondsList, std::list<bond::ptr >::const_iterator SetFirst, std::list<bond::ptr >::const_iterator SetLast, const TouchedList_t &TouchedList, int TouchedIndex)
{
  int SetDimension = 0;
  for( std::list<bond::ptr >::const_iterator Binder = SetFirst;
      Binder != SetLast;
      ++Binder) {
    for (TouchedList_t::const_iterator iter = TouchedList.begin();
        iter != TouchedList.end(); ++iter) {
      if ((*Binder)->leftatom->getNr() == *iter)   // 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, TouchedList_t &TouchedList, int &TouchedIndex)
{
  for (TouchedList_t::iterator iter = TouchedList.begin();
      iter != TouchedList.end();++iter) {
    const int Removal = *iter;
    LOG(2+verbosity, "Removing item nr. " << Removal << " from snake stack.");
    FragmentSet->erase(Removal);
    (*iter) = -1;
  }
  std::stringstream output;
  output << "Remaining local nr.s on snake stack are: ";
  for(KeySet::iterator runner = FragmentSet->begin(); runner != FragmentSet->end(); runner++)
    output << (*runner) << " ";
  LOG(2, output.str());
  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 treatment 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 HydrogenTreatment treatment)
{
  int Counter = FragmentSearch->FragmentCounter; // mark current value of counter

  LOG(0, std::endl);
  LOG(0, "Begin of PowerSetGenerator with order " << BondsPerSPList.getOrder() << " at Root " << *FragmentSearch->getRoot() << ".");

  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, treatment);

  // 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();
  LOG(0, "INFO: Total number of edges is " << SP << ".");
  {
    // start with root (push on fragment stack)
    LOG(0, "Starting fragment generation with " << *FragmentSearch->getRoot() << ", local nr is " << FragmentSearch->getRoot()->getNr() << ".");
    FragmentSearch->FragmentSet->clear();
    LOG(0, "Preparing subset for this root and calling generator.");

    // prepare the subset and call the generator
    std::vector<bond::ptr> 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());
  }

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

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

  // return list
  LOG(0, "End of PowerSetGenerator.");
  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::ptr > &BondsSet, int SetDimension, int SubOrder)
{
  Info info(__func__);
  int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
  int TouchedIndex, SubSetDimension, Added;
  int SpaceLeft;
  TouchedList_t TouchedList(SubOrder + 1, -1);
  KeySetTestPair TestKeySetInsert;

  const int 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(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).");

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

  // create every possible combination of the endpieces
  LOG(1+verbosity, "Going through all combinations of the power set.");
  for (int i=1;i<NumCombinations;i++) {  // sweep through all power set combinations (skip empty set!)
    // count the set bit of i
    const int bits = countBits(i, SetDimension);

    LOG(1+verbosity, "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << ".");
    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);

      // --2-- store the fragment
      const int order = BondsPerSPList.getOrder() - SubOrder + Added - 1;
      LOG(1+verbosity, "Storing fragment as order " << order << ".");
      // store fragment as a KeySet
      if (DoLog(2)) {
        std::stringstream output;
        output << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
        for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
          output << (*runner) << " ";
        LOG(2, output.str());
      }
      FragmentSearch->InsertFragmentIntoGraph(order);

      // --3-- if below SubOrder still, recurse
      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
      if ((SpaceLeft > 0) && (RootDistance < BondsPerSPList.getOrder())) {
        LOG(1+verbosity, "There's still some space left on stack: " << SpaceLeft << ".");
        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
          const int 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::ptr > BondsList;
          BondsList.resize(SubSetDimension);
          SubSetDimension = FillBondsList(BondsList, BondsPerSPList.BondsPerSPList[SP].begin(), BondsPerSPList.BondsPerSPList[SP].end(), TouchedList, TouchedIndex);

          // then iterate
          LOG(2+verbosity, "Calling subset generator " << SP << " away from root " << *FragmentSearch->getRoot() << " with sub set dimension " << SubSetDimension << ".");
          SPFragmentGenerator(SP, BondsList, SubSetDimension, SubOrder-bits);
        }
      } else {
        LOG(1+verbosity, "No more space or no shortest path levels, not recursing.");
      }

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