source: src/Graph/AdjacencyList.cpp@ 97a858

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

/*
 * AdjacencyList.cpp
 *
 *  Created on: Mar 3, 2011
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <iostream>
#include <map>
#include <set>
#include <utility>

#include "AdjacencyList.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Helpers/defs.hpp"
#include "World.hpp"

/** Constructor of class AdjacencyList.
 *
 * \param File file to parser
 */
AdjacencyList::AdjacencyList(std::istream &File)
{
  bool status = ParseFromFile(File);
  ASSERT( status,
      "AdjacencyList::AdjacencyList() - File's contents was nor parsable");
  if (!status) // remove map if failed to parse
    atombondmap.clear();
}

/** Constructor of class AdjacencyList.
 *
 * \param File file to parser
 */
AdjacencyList::AdjacencyList(const atomids_t &atomids)
{
  CreateMap(atomids);
}

AdjacencyList::~AdjacencyList()
{}

/** Parses the bond partners of each atom from an external file into AdjacencyList::atombondmap.
 *
 * @param File file to parse
 * @return true - everything ok, false - error while parsing
 */
bool AdjacencyList::ParseFromFile(std::istream &File)
{
  if (File.fail()) {
    LOG(1, "STATUS: Adjacency file not found." << endl);
    return false;
  }

  atombondmap.clear();
  char buffer[MAXSTRINGSIZE];
  int tmp;
  // Parse the file line by line and count the bonds
  while (File.good()) {
    File.getline(buffer, MAXSTRINGSIZE);
    stringstream line;
    line.str(buffer);
    int AtomNr = -1;
    line >> AtomNr;
    // parse into structure
    if (AtomNr > 0) {
      const atomId_t WalkerId = AtomNr-1;
      // parse bond partner ids associated to AtomNr
      while (line >> ws >> tmp) {
        LOG(3, "INFO: Recognized bond partner " << tmp-1 << " for " << WalkerId << ".");
        atombondmap.insert( std::make_pair(WalkerId, tmp-1) );
      }
    } else {
      if (AtomNr != -1) {
        ELOG(2, AtomNr << " is negative.");
        return false;
      }
    }
  }
  return true;
}

/** Fills the AdjacencyList::atombondmap from the atoms given by the two iterators.
 *
 * @param atomids set of atomic ids to check (must be global ids, i.e. from atom::getId())
 */
void AdjacencyList::CreateMap(atomids_t atomids)
{
  atombondmap.clear();
  std::sort(atomids.begin(), atomids.end());
  // go through each atom in the list
  for (atomids_t::const_iterator iter = atomids.begin(); iter != atomids.end(); ++iter) {
    const atomId_t WalkerId = *iter;
    ASSERT(WalkerId != (size_t)-1,
        "AdjacencyList::CreateMap() - Walker has no id.");
    const atom *Walker = World::getInstance().getAtom(AtomById(WalkerId));
    ASSERT( Walker != NULL,
        "AdjacencyList::CreateMap() - Walker id "+toString(*iter)
        +" is not associated to any of World's atoms.");
    const BondList& ListOfBonds = Walker->getListOfBonds();
    // go through each of its bonds
    for (BondList::const_iterator Runner = ListOfBonds.begin();
        Runner != ListOfBonds.end();
        ++Runner) {
      const atomId_t id = (*Runner)->GetOtherAtom(Walker)->getId();
      ASSERT(id != (size_t)-1,
          "AdjacencyList::CreateMap() - OtherAtom has not id.");
      if (std::binary_search(atomids.begin(), atomids.end(), id))
        atombondmap.insert( std::make_pair(WalkerId, id) );
    }
  }
}

AdjacencyList::KeysSet AdjacencyList::getKeys(const AdjacencyList::AtomBondRange &_range) const
{
  KeysSet Keys;
  for (AtomBondMap::const_iterator iter = _range.first;
      iter != _range.second;
      ++iter) {
    Keys.insert( iter->first );
  }
  return Keys;
}

AdjacencyList::ValuesSet AdjacencyList::getValues(const AdjacencyList::AtomBondRange&_range) const
{
  ValuesSet Values;
  for (AtomBondMap::const_iterator iter = _range.first;
      iter != _range.second;
      ++iter) {
    Values.insert( iter->second );
  }
  return Values;
}

/** Counts the number of items in the second set not present in the first set.
 *
 * \note We assume that the sets are sorted.
 *
 * @param firstset check set
 * @param secondset reference set
 * @return number of items in the first set that are missing in the second
 */
template <class T>
size_t getMissingItems(const T &firstset, const T &secondset)
{
  size_t Mismatch = 0;
  typename T::const_iterator firstiter = firstset.begin();
  typename T::const_iterator seconditer = secondset.begin();
  for (; (firstiter != firstset.end()) && (seconditer != secondset.end());) {
    if (*firstiter > *seconditer)
      ++seconditer;
    else {
      if (*firstiter < *seconditer)
        ++Mismatch;
      ++firstiter;
    }
  }
  return Mismatch;
}

/** Compares whether AdjacencyList::atombondmap in this instance is a subset of
 * the one in \a other.
 *
 * @param other other instance of an adjacency list to compare against
 * @return true - this atombondmap is subset, false - else
 */
bool AdjacencyList::operator<(const AdjacencyList &other) const
{
  int NonMatchNumber = 0;
  // extract keys and check whether they match
  const AtomBondRange Intrange(atombondmap.begin(), atombondmap.end());
  const AtomBondRange Extrange(other.atombondmap.begin(), other.atombondmap.end());
  KeysSet InternalKeys( getKeys(Intrange) );
  KeysSet ExternalKeys( getKeys(Extrange) );

//  std::cout << "InternalKeys: " << InternalKeys << std::endl;
//  std::cout << "ExternalKeys: " << ExternalKeys << std::endl;

  // check amount of keys
  if (ExternalKeys.size() < InternalKeys.size()) {
    NonMatchNumber = (int)InternalKeys.size() - (int)ExternalKeys.size();
    LOG(2, "INFO: Number of internal keys exceeds external one by " << NonMatchNumber << ".");
    return false;
  }

  // check which keys are missing in the internal set
  NonMatchNumber = getMissingItems(InternalKeys, ExternalKeys);

  if (NonMatchNumber != 0) {
    LOG(2, "INFO: " << NonMatchNumber << " keys are not the same.");
    return false;
  }

  // now check each map per key
  for (KeysSet::const_iterator keyIter = InternalKeys.begin();
      keyIter != InternalKeys.end();
      ++keyIter) {
//    std::cout << "Current key is " << *keyIter << std::endl;
    const AtomBondRange IntRange( atombondmap.equal_range(*keyIter) );
    const AtomBondRange ExtRange( other.atombondmap.equal_range(*keyIter) );
    ValuesSet InternalValues( getValues(IntRange) );
    ValuesSet ExternalValues( getValues(ExtRange) );
    // throw out all values not present in ExternalKeys
    ValuesSet ExternalValues_temp( ExternalValues );
    for(KeysSet::const_iterator iter = InternalKeys.begin();
        iter != InternalKeys.end(); ++iter)
      ExternalValues_temp.erase(*iter);
    // all remaining values must be masked out
    for (ValuesSet::const_iterator iter = ExternalValues_temp.begin();
        iter != ExternalValues_temp.end(); ++iter)
      ExternalValues.erase(*iter);
//    std::cout << "InternalValues: " << InternalValues << std::endl;
//    std::cout << "ExternalValues: " << ExternalValues << std::endl;
    NonMatchNumber += getMissingItems(InternalValues, ExternalValues);
  }
  if (NonMatchNumber != 0) {
    LOG(2, "INFO: " << NonMatchNumber << " keys are not the same.");
    return false;
  } else {
    LOG(2, "INFO: All keys are the same.");
    return true;
  }
}

/** Storing the bond structure of a molecule to file.
 * Simply stores Atom::Nr and then the Atom::Nr of all bond partners per line.
 * @param File output stream to write to
 * \return true - file written successfully, false - writing failed
 */
bool AdjacencyList::StoreToFile(std::ostream &File) const
{
  bool status = true;
  std::stringstream output;
  output << "Saving adjacency list ... ";
  if (!File.bad()) {
    //File << "m\tn" << endl;
    AtomBondMap::const_iterator advanceiter = atombondmap.begin();
    for (AtomBondMap::const_iterator iter = atombondmap.begin();
        iter != atombondmap.end();
        iter = advanceiter) {
      advanceiter = atombondmap.upper_bound(iter->first);
      File << iter->first+1;
      for (AtomBondMap::const_iterator adjacencyiter = iter;
          adjacencyiter != advanceiter;
          ++adjacencyiter)
        File << " " << adjacencyiter->second+1;
      File << std::endl;
    }
    output << "done.";
  } else {
    output << "given stream is not good.";
    status = false;
  }
  LOG(1, output.str());

  return status;
}