source: src/Graph/CheckAgainstAdjacencyFile.cpp@ b4f72c

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

/*
 * CheckAgainstAdjacencyFile.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 "CheckAgainstAdjacencyFile.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 CheckAgainstAdjacencyFile.
 *
 * \param File file to parser
 */
CheckAgainstAdjacencyFile::CheckAgainstAdjacencyFile(std::istream &File) :
  NonMatchNumber(0)
{
  bool status = ParseInInternalMap(File);
  if (!status) // remove map if failed to parse
    InternalAtomBondMap.clear();
}

CheckAgainstAdjacencyFile::~CheckAgainstAdjacencyFile()
{}

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

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

/** Fills the InternalAtomBondMap 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 CheckAgainstAdjacencyFile::CreateExternalMap(const atomids_t &atomids)
{
  ExternalAtomBondMap.clear();
  // 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,
        "CheckAgainstAdjacencyFile::CreateExternalMap() - Walker has no id.");
    const atom *Walker = World::getInstance().getAtom(AtomById(WalkerId));
    ASSERT( Walker != NULL,
        "CheckAgainstAdjacencyFile::CreateExternalMap() - 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,
          "CheckAgainstAdjacencyFile::CreateExternalMap() - OtherAtom has not id.");
      ExternalAtomBondMap.insert( std::make_pair(WalkerId, id) );
    }
  }
}

/** Checks contents of adjacency file against bond structure in structure molecule.
 * \return true - structure is equal, false - not equivalence
 */
bool CheckAgainstAdjacencyFile::operator()(const atomids_t &atomids)
{
  LOG(0, "STATUS: Looking at bond structure of given ids and comparing against stored in adjacency file... ");

  // parse in external map
  CreateExternalMap(atomids);

  bool status = CompareInternalExternalMap();
  if (status) { // if equal we parse the KeySetFile
    LOG(0, "STATUS: Equal.");
  } else
    LOG(0, "STATUS: Not equal by " << NonMatchNumber << " atoms.");
  return status;
}

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

CheckAgainstAdjacencyFile::ValuesSet CheckAgainstAdjacencyFile::getValues(const CheckAgainstAdjacencyFile::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 InternalAtomBondMap and ExternalAtomBondMap and sets NonMatchNumber.
 *
 * @return true - both maps are the same, false - both maps diverge by NonMatchNumber counts.
 */
bool CheckAgainstAdjacencyFile::CompareInternalExternalMap()
{
  NonMatchNumber = 0;
  // extract keys and check whether they match
  const AtomBondRange Intrange(InternalAtomBondMap.begin(), InternalAtomBondMap.end());
  const AtomBondRange Extrange(ExternalAtomBondMap.begin(), ExternalAtomBondMap.end());
  KeysSet InternalKeys( getKeys(Intrange) );
  KeysSet ExternalKeys( getKeys(Extrange) );

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

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

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

  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 = ExternalKeys.begin();
      keyIter != ExternalKeys.end();
      ++keyIter) {
//    std::cout << "Current key is " << *keyIter << std::endl;
    const AtomBondRange IntRange( InternalAtomBondMap.equal_range(*keyIter) );
    const AtomBondRange ExtRange( ExternalAtomBondMap.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 = ExternalKeys.begin();
        iter != ExternalKeys.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(ExternalValues, InternalValues);
  }
  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;
  }
}