source: src/Fragmentation/Summation/SubsetMap.cpp@ e13990b

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. 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 <algorithm>
#include <bitset>
#include <set>

#include "CodePatterns/Log.hpp"

SubsetMap::SubsetMap(IndexSetContainer &_container) :
  tempset(new IndexSet())
{
  const IndexSetContainer::Container_t &container = _container.getContainer();

// we don't need this level-wise sorting as this is automatically achieved by specific sorting
//  /// place in level-wise multimap
//  typedef size_t level_t;
//  // note: We use an IndexSetContainer as it is not accessed while filling it
//  // also, as insert uses push_back, the IndexSetContainer is actually created
//  // as sorted. Hence, the call to sort() just needs linear time (and does
//  // nothing)
//  typedef std::map<level_t, IndexSetContainer > IndexSetPerLevel_t;
//  IndexSetPerLevel_t IndexSetPerLevel;
//  for(IndexSetContainer::Container_t::const_iterator iter = container.begin();
//      iter != container.end(); ++iter) {
//    // key present in map?
//    const size_t level = (*iter)->size();
//    LOG(1, "INFO: Current set is " << **iter << " with size " << level << ".");
//    IndexSetPerLevel_t::iterator leveliter = IndexSetPerLevel.find(level);
//    // we have to explicitly copy the shared_ptr for insertion due to const container
//    IndexSet_ptr ptr(*iter);
//    if (leveliter == IndexSetPerLevel.end()) {
//      LOG(2, "DEBUG: Level not present in Lookup, creating new container");
//      IndexSetPerLevel.insert( std::make_pair( level, IndexSetContainer(ptr)) );
//    } else {
//      LOG(2, "DEBUG: Level present in Lookup, appending to container");
//      leveliter->second.insert(ptr);
//    }
//  }

  /// go through all IndexSets
  std::for_each( container.begin(), container.end(),
      boost::bind(&SubsetMap::gatherSubsets, this, _1));
}

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(1, "INFO: Current set to gather all subsets for is " << *ptr << " with size " << SetSize << ".");
  const size_t NoPowerSubsets = getNoPowerSets(SetSize);
  LOG(1, "INFO: 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(2, "INFO: Current subset is " << *tempset << " at " << index << ".");
    // check whether IndexSet is contained and known
    if (ptr->contains(*tempset)) {
      // mark as contained
      *iter = Contained;
      LOG(3, "DEBUG: Marking subset as Contained.");
      if (Lookup.count(tempset)) {
        LOG(2, "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(2, "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(4, "INFO: Current subset of subset is " << **containeriter << " at " << subindex << ".");
            if (sets[subindex] != Contained) {
              LOG(5, "DEBUG: Setting power subset to Contained.");
              sets[subindex] = Contained;
              returnsets.push_back(**containeriter);
            }
          }
        }
      }
    }
  }
  LOG(1, "INFO: 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;
}