source: src/Fragmentation/AdaptivityMap.cpp@ d4d7a1

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

/*
 * AdaptivityMap.cpp
 *
 *  Created on: Oct 20, 2011
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "AdaptivityMap.hpp"

#include <fstream>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

#include "Atom/atom.hpp"
#include "Fragmentation/AtomMask.hpp"
#include "Helpers/defs.hpp"
#include "Helpers/helpers.hpp"
#include "molecule.hpp"

/** Constructor of class AdaptivityMap.
 *
 */
AdaptivityMap::AdaptivityMap()
{}

/** Destructor of class AdaptivityMap.
 *
 */
AdaptivityMap::~AdaptivityMap()
{}

/** Inserts a (\a No, \a value) pair into the list, overwriting present one.
 * Note if values are equal, No will decided on which is first
 * \param *out output stream for debugging
 * \param &AdaptiveCriteriaList list to insert into
 * \param &IndexedKeySetList list to find key set for a given index \a No
 * \param FragOrder current bond order of fragment
 * \param No index of keyset
 * \param value energy value
 */
void AdaptivityMap::InsertIntoAdaptiveCriteriaList(int FragOrder, int No, double Value)
{
  ASSERT( AdaptiveCriteriaList != NULL,
      "AdaptivityMap::InsertIntoAdaptiveCriteriaList() - AdaptiveCriteriaList is not allocated yet.");
  const_iterator marker = find(No);    // find keyset to Frag No.
  if (marker != end()) {  // if found
    Value *= 1 + MYEPSILON*(*((*marker).second.begin()));     // in case of equal energies this makes them not equal without changing anything actually
    // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
    std::pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList->insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
    std::map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
    if (!InsertedElement.second) { // this root is already present
      if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
        //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
        {  // if value is smaller, update value and order
        (*PresentItem).second.first = fabs(Value);
        (*PresentItem).second.second = FragOrder;
        LOG(2, "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "]).");
      } else {
        LOG(2, "Did not update element " <<  (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << ".");
      }
    } else {
      LOG(2, "Inserted element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "]).");
    }
  } else {
    LOG(1, "No Fragment under No. " << No << "found.");
  }
};


/** Scans the adaptive order file and insert (index, value) into map.
 * \param &path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
 */
void AdaptivityMap::ScanAdaptiveFileIntoMap(std::string &path)
{
  int No = 0, FragOrder = 0;
  double Value = 0.;
  char buffer[MAXSTRINGSIZE];
  std::string filename = path + ENERGYPERFRAGMENT;
  std::ifstream InputFile(filename.c_str());

  if (InputFile.fail()) {
    ELOG(1, "Cannot find file " << filename << ".");
    return;
  }

  if (CountLinesinFile(InputFile) > 0) {
    // each line represents a fragment root (Atom::Nr) id and its energy contribution
    InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
    InputFile.getline(buffer, MAXSTRINGSIZE);
    while(!InputFile.eof()) {
      InputFile.getline(buffer, MAXSTRINGSIZE);
      if (strlen(buffer) > 2) {
        //LOG(2, "Scanning: " << buffer);
        stringstream line(buffer);
        line >> FragOrder;
        line >> ws >> No;
        line >> ws >> Value; // skip time entry
        line >> ws >> Value;
        No -= 1;  // indices start at 1 in file, not 0
        //LOG(2, " - yields (" << No << "," << Value << ", " << FragOrder << ")");

        // clean the list of those entries that have been superceded by higher order terms already
        InsertIntoAdaptiveCriteriaList(FragOrder, No, Value);
      }
    }
    // close and done
    InputFile.close();
    InputFile.clear();
  }
};

/** Maps adaptive criteria list back onto (Value, (Root Nr., Order))
 * (i.e. sorted by value to pick the highest ones)
 * \param *mol molecule with atoms
 * \return remapped list
 */
void AdaptivityMap::ReMapAdaptiveCriteriaListToValue(molecule *mol)
{
  atom *Walker = NULL;
  ASSERT( AdaptiveCriteriaList != NULL,
      "AdaptivityMap::ReMapAdaptiveCriteriaListToValue() - AdaptiveCriteriaList is not allocated yet.");
  FinalRootCandidates = new AdaptiveCriteriaValueMap;
  LOG(1, "Root candidate list is: ");
  for(AdaptiveCriteriaIndexMap::const_iterator runner = AdaptiveCriteriaList->begin(); runner != AdaptiveCriteriaList->end(); runner++) {
    Walker = mol->FindAtom((*runner).first);
    if (Walker != NULL) {
      //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
      if (!Walker->MaxOrder) {
        LOG(2, "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])");
        FinalRootCandidates->insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
      } else {
        LOG(2, "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order.");
      }
    } else {
      ELOG(0, "Atom No. " << (*runner).second.first << " was not found in this molecule.");
      performCriticalExit();
    }
  }
};

/** Counts lines in file.
 * Note we are scanning lines from current position, not from beginning.
 * \param InputFile file to be scanned.
 */
int AdaptivityMap::CountLinesinFile(std::ifstream &InputFile) const
{
  char *buffer = new char[MAXSTRINGSIZE];
  int lines=0;

  int PositionMarker = InputFile.tellg();  // not needed as Inputfile is copied, given by value, not by ref
  // count the number of lines, i.e. the number of fragments
  InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
  InputFile.getline(buffer, MAXSTRINGSIZE);
  while(!InputFile.eof()) {
    InputFile.getline(buffer, MAXSTRINGSIZE);
    lines++;
  }
  InputFile.seekg(PositionMarker, ios::beg);
  delete[](buffer);
  return lines;
};

/** Marks all candidate sites for update if below adaptive threshold.
 * Picks a given number of highest values and set *AtomMask to true.
 * \param AtomMask defines true/false per global Atom::Nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
 * \param Order desired order
 * \param *mol molecule with atoms
 * \return true - if update is necessary, false - not
 */
bool AdaptivityMap::MarkUpdateCandidates(AtomMask_t &AtomMask, int Order, molecule *mol) const
{
  atom *Walker = NULL;
  int No = -1;
  bool status = false;
  ASSERT( FinalRootCandidates != NULL,
      "AdaptivityMap::MarkUpdateCandidates() - FinalRootCandidates is not allocated yet.");
  for(AdaptiveCriteriaValueMap::const_iterator runner = FinalRootCandidates->upper_bound(pow(10.,Order)); runner != FinalRootCandidates->end(); runner++) {
    No = (*runner).second.first;
    Walker = mol->FindAtom(No);
    //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->getNr()]) {
      LOG(2, "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true.");
      AtomMask.setTrue(No);
      status = true;
    //} else
      //LOG(2, "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->getNr()] << " does not allow further adaptive increase.");
  }
  return status;
};

/** Checks whether there are any adaptive items currently.
 *
 * @return true - there are items for adaptive refinement, false - there are none
 */
bool AdaptivityMap::IsAdaptiveCriteriaListEmpty() const
{
  ASSERT( AdaptiveCriteriaList != NULL,
      "AdaptivityMap::ReMapAdaptiveCriteriaListToValue() - AdaptiveCriteriaList is not allocated yet.");
  return AdaptiveCriteriaList->empty();
}