/* * 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 . */ /* * CheckAgainstAdjacencyFile.cpp * * Created on: Mar 3, 2011 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include #include #include #include #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 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; } }