source: src/Fragmentation/Summation/SubsetMap.cpp@ 13e5be

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  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/>.
 */

/*
 * SubsetMap.cpp
 *
 *  Created on: Jul 3, 2012
 *      Author: heber
 */


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

//#include "CodePatterns/MemDebug.hpp"

#include "SubsetMap.hpp"

#include <boost/bind.hpp>
#include <boost/foreach.hpp>
#include <boost/lambda/lambda.hpp>
#include <algorithm>
#include <bitset>
#include <set>

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

SubsetMap::SubsetMap(IndexSetContainer &_container) :
  tempset(new IndexSet())
{
  /// go through all IndexSets
  const IndexSetContainer::Container_t &container = _container.getContainer();
  for (IndexSetContainer::Container_t::const_iterator iter = container.begin();
      iter != container.end(); ++iter) {
    // check whether its the super set
    if (*iter == _container.getSuperSet()) {
      // place all other sets as its subsets
      std::vector<IndexSet> returnsets;
      returnsets.reserve(container.size());
      for (IndexSetContainer::Container_t::const_iterator insertiter = container.begin();
            insertiter != container.end(); ++insertiter)
        if (iter != insertiter) {
          LOG(2, "INFO: Current subset is " << **insertiter << " for super set.");
          ASSERT( _container.getSuperSet()->contains(**insertiter),
              "SubsetMap::SubsetMap() - "+toString(**insertiter)+" is not contained in super set "
              +toString(*_container.getSuperSet())+".");
          returnsets.push_back(**insertiter);
        }
      Lookup.insert( std::make_pair(*iter, IndexSetContainer::ptr(new IndexSetContainer(returnsets))) );
    } else {
      gatherSubsets(*iter);
    }
  }
}

void SubsetMap::gatherSubsets(const IndexSet_ptr &ptr) {
  // get the sets size and the size of the power sets
  const size_t SetSize = ptr->size();
  LOG(2, "DEBUG: Current set to gather all subsets for is " << *ptr << " with size " << SetSize << ".");
  const size_t NoPowerSubsets = getNoPowerSets(SetSize);
  LOG(2, "DEBUG: There are " << NoPowerSubsets << " sets to check for containment.");
  std::vector<IndexSet> returnsets;
  // sets is the number of possible subsets to check: NotContained means still
  // needs check or is not contained; Contained means is contained, or has
  // already been treated via a super(sub)set.
  typedef std::vector<enum ContainedState> StateOfSubsets;
  StateOfSubsets sets( NoPowerSubsets, NotContained);
  // we have to traverse backwards, i.e. from greatest subsets to smallest to
  // make the best possible use of the subset knowledge we already have for these
  // subsets of the subset.
  // WARNING: last set is the set itself which must not be in IndexSetContainer due to
  // cyclic structure which shared_ptr cannot resolve.
  StateOfSubsets::reverse_iterator iter = sets.rbegin()+1;
  for (;iter != sets.rend(); iter = std::find(iter+1, sets.rend(), NotContained)) {
    // get IndexSet (get index via distance)
    // NOTE: reverse iterator naturally give reverse distance, hence flip both arguments
    // also, reverse iterator is off by one, see ASSERT
    const size_t index = std::distance(iter, sets.rend())-1;
    ASSERT( sets.at(index) == *iter,
        "SubsetMap::gatherSubsets() - distance from begin to iterator is not correct.");
    *tempset = getSubset(*ptr, index );
    LOG(3, "DEBUG: Current subset is " << *tempset << " at " << index << ".");
    // check whether IndexSet is contained and known
    if (ptr->contains(*tempset)) {
      // mark as contained
      *iter = Contained;
      LOG(4, "DEBUG: Marking subset as Contained.");
      if (Lookup.count(tempset)) {
        LOG(3, "DEBUG: Subset is present in Lookup, adding.");
        returnsets.push_back(*tempset);
        // if so, also add its contained sets and mark them out
        const IndexSetContainer::Container_t &container = getSubsets(tempset)->getContainer();
        if (container.size()) {
          LOG(3, "DEBUG: Subset is also present in Lookup, adding all its subsets");
          for(IndexSetContainer::Container_t::const_iterator containeriter = container.begin();
              containeriter != container.end();
              ++containeriter) {
            const size_t subindex = getPowerSetIndex(ptr, **containeriter);
            LOG(5, "INFO: Current subset of subset is " << **containeriter << " at " << subindex << ".");
            if (sets[subindex] != Contained) {
              LOG(6, "DEBUG: Setting power subset to Contained.");
              sets[subindex] = Contained;
              returnsets.push_back(**containeriter);
            }
          }
        }
      }
    }
  }
  LOG(2, "DEBUG: Adding " << returnsets.size() << " sets to Lookup for set " << *ptr << ".");
  Lookup.insert( std::make_pair(ptr, IndexSetContainer::ptr(new IndexSetContainer(returnsets))) );
}

const size_t SubsetMap::getPowerSetIndex(const IndexSet_ptr &ptr, const IndexSet &subset) {
  ASSERT( (ptr->size() < MAX_SETSIZE) && (subset.size() < MAX_SETSIZE),
      "SubsetMap::getPowerSetIndex() - power sets of subsets must not be bigger than "
      +toString((int)MAX_SETSIZE)+".");
  ASSERT( ptr->contains(subset),
      "SubsetMap::getPowerSetIndex() - "+(toString(*ptr))+" does not contain "
          +toString(subset)+".");
  std::bitset<MAX_SETSIZE> bits(0);
  // go through each and search for indices present in both
  IndexSet::iterator indexiter = ptr->begin();
  IndexSet::iterator subindexiter = subset.begin();
  for (; subindexiter != subset.end();) {
    if (*indexiter == *subindexiter) {
      // if matching, set bit and advance both iterators
      bits.set( (size_t)std::distance(ptr->begin(), indexiter) );
      ++indexiter;
      ++subindexiter;
      // if not, advance the iterator lacking behind
    } else if ( *indexiter < *subindexiter) {
      ++indexiter;
    } else {
      ++subindexiter;
    }
  }
  return bits.to_ulong();
}

IndexSet SubsetMap::getSubset(const IndexSet &_set, size_t PowerSetNo) {
  ASSERT( _set.size() < MAX_SETSIZE,
      "SubsetMap::getSubset() - power sets of subsets must not be bigger than "
      +toString((int)MAX_SETSIZE)+".");
  ASSERT((PowerSetNo >= 0) && (PowerSetNo < getNoPowerSets(_set.size())),
      "SubsetMap::getSubset() - Power set index "+toString(PowerSetNo)+" is out of range.");
  std::bitset<MAX_SETSIZE> bits(PowerSetNo);
  // place index into container with random access
  std::vector<IndexSet::value_type> indices(_set.begin(), _set.end());
  IndexSet indexset;
  // go through each bit and insert if 1, not if 0
  for (size_t index=0; index < indices.size(); ++index)
    if (bits[index])
      indexset.insert(indices[index]);
  return indexset;
}

size_t SubsetMap::getMaximumSetLevel() const
{
  if (Lookup.empty())
    return 0;
  // last one is super set
  Lookup_t::const_iterator iter = --Lookup.end();
  return iter->first->size();
}