/* * 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 . */ /* * AdjacencyList.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 "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 * const Walker = const_cast(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(); if (ListOfBonds.size() > 0) { // 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) ); } } else { // insert self-bond if atom is lonesome cowboy atombondmap.insert( std::make_pair(WalkerId, WalkerId) ); } } } 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 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; }