source: src/Graph/CheckAgainstAdjacencyFile.cpp@ b78dfd

/*
 * 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) :
  status(true),
  NonMatchNumber(0)
{
  ParseInExternalMap(File);
}

CheckAgainstAdjacencyFile::~CheckAgainstAdjacencyFile()
{
  ExternalAtomBondMap.clear();
  InternalAtomBondMap.clear();
}

/** 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::ParseInExternalMap(std::istream &File)
{
  if (File.fail()) {
    LOG(1, "STATUS: Adjacency file not found." << endl);
    return false;
  }

  ExternalAtomBondMap.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));
      ASSERT(Walker != NULL,
          "CheckAgainstAdjacencyFile::ParseInExternalMap() - there is no atom with id "+toString(AtomNr-1)+".");
      if (Walker == NULL)
        return false;
      // parse bond partner ids associated to AtomNr
      while (line >> ws >> tmp) {
        LOG(3, "INFO: Recognized bond partner " << tmp-1);
        ExternalAtomBondMap.insert( std::make_pair(Walker->getId(), 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 AtomMapBegin iterator pointing to begin of map (think of World's SelectionIterator)
 * @param AtomMapEnd iterator pointing past end of map (think of World's SelectionIterator)
 */
void CheckAgainstAdjacencyFile::CreateInternalMap(World::AtomSet::const_iterator AtomMapBegin, World::AtomSet::const_iterator AtomMapEnd)
{
  InternalAtomBondMap.clear();
  // go through each atom in the list
  for (World::AtomSet::const_iterator iter = AtomMapBegin; iter != AtomMapEnd; ++iter) {
    const atom *Walker = iter->second;
    const atomId_t WalkerId = Walker->getId();
    ASSERT(WalkerId != (size_t)-1,
        "CheckAgainstAdjacencyFile::CreateInternalMap() - Walker has no id.");
    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::CreateInternalMap() - OtherAtom has not id.");
      InternalAtomBondMap.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()(World::AtomSet::const_iterator AtomMapBegin, World::AtomSet::const_iterator AtomMapEnd)
{
  LOG(0, "STATUS: Looking at bond structure stored in adjacency file and comparing to present one ... ");

  bool status = true;

  if (InternalAtomBondMap.empty())
    CreateInternalMap(AtomMapBegin, AtomMapEnd);

  status = 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 mismatching items in each set.
 *
 * @param firstset first set
 * @param secondset second set
 * @return number of items that don't match between first and second set
 */
template <class T>
size_t getMismatchingItems(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());
      ++firstiter, ++seconditer) {
    if (*firstiter != *seconditer)
      ++Mismatch;
  }
  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;
  // check whether sizes match
  if (ExternalAtomBondMap.size() != InternalAtomBondMap.size()) {
    NonMatchNumber = abs((int)ExternalAtomBondMap.size() - (int)InternalAtomBondMap.size());
    LOG(2, "INFO: " << NonMatchNumber << " entries don't match.");
    return false;
  }
  // 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 (InternalKeys.size() != ExternalKeys.size()) {
    NonMatchNumber = abs((int)ExternalKeys.size() - (int)InternalKeys.size());
    LOG(2, "INFO: Number of keys don't match: "
        << InternalKeys.size() << " != " << ExternalKeys.size());
    return false;
  }

  // check items against one another
  NonMatchNumber = getMismatchingItems(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( InternalAtomBondMap.equal_range(*keyIter) );
    const AtomBondRange ExtRange( ExternalAtomBondMap.equal_range(*keyIter) );
    ValuesSet InternalValues( getValues(IntRange) );
    ValuesSet ExternalValues( getValues(ExtRange) );
//    std::cout << "InternalValues: " << InternalValues << std::endl;
//    std::cout << "ExternalValues: " << ExternalValues << std::endl;
    NonMatchNumber += getMismatchingItems(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;
  }
}