source: src/Analysis/analysis_correlation.hpp@ 4464ef

/*
 * analysis.hpp
 *
 *  Created on: Oct 13, 2009
 *      Author: heber
 */

#ifndef ANALYSIS_HPP_
#define ANALYSIS_HPP_

using namespace std;

/*********************************************** includes ***********************************/

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

#include <cmath>
#include <iomanip>
#include <fstream>

// STL headers
#include <map>
#include <vector>

#include "Helpers/defs.hpp"

#include "Atom/atom.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Helpers/helpers.hpp"
#include "World.hpp"

/****************************************** forward declarations *****************************/

class BoundaryTriangleSet;
class Formula;
class element;
class LinkedCell_deprecated;
class Tesselation;
class Vector;

/********************************************** definitions *********************************/

typedef multimap<double, pair<const TesselPoint *, const TesselPoint *> > PairCorrelationMap;
typedef multimap<double, const atom * > DipoleAngularCorrelationMap;
typedef multimap<double, pair<const molecule *, const molecule *> > DipoleCorrelationMap;
typedef multimap<double, pair<const atom *, const Vector *> > CorrelationToPointMap;
typedef multimap<double, pair<const atom *, BoundaryTriangleSet *> > CorrelationToSurfaceMap;
typedef map<double, int> BinPairMap;

enum ResetWorldTime {
  DontResetTime,
  DoResetTime
};

/********************************************** declarations *******************************/

range<size_t> getMaximumTrajectoryBounds(
    const std::vector<const atom *> &atoms);
std::map<atomId_t, Vector> CalculateZeroAngularDipole(
    const std::vector<const molecule *> &molecules);

DipoleAngularCorrelationMap *DipoleAngularCorrelation(
    const Formula &DipoleFormula,
    const size_t timestep,
    const std::map<atomId_t, Vector> &ZeroVector,
    const enum ResetWorldTime DoTimeReset = DontResetTime);
DipoleCorrelationMap *DipoleCorrelation(
    const std::vector<const molecule *> &molecules);
PairCorrelationMap *PairCorrelation(
    const World::ConstAtomComposite &atoms_first,
    const World::ConstAtomComposite &atoms_second,
    const double max_distance);
CorrelationToPointMap *CorrelationToPoint(
    const std::vector<const molecule *> &molecules,
    const std::vector<const element *> &elements,
    const Vector *point );
CorrelationToSurfaceMap *CorrelationToSurface(
    const std::vector<const molecule *> &molecules,
    const std::vector<const element *> &elements,
    const Tesselation * const Surface,
    const LinkedCell_deprecated *LC );
CorrelationToPointMap *PeriodicCorrelationToPoint(
    const std::vector<const molecule *> &molecules,
    const std::vector<const element *> &elements,
    const Vector *point, const int ranges[NDIM] );
CorrelationToSurfaceMap *PeriodicCorrelationToSurface(
    const std::vector<const molecule *> &molecules,
    const std::vector<const element *> &elements,
    const Tesselation * const Surface,
    const LinkedCell_deprecated *LC,
    const int ranges[NDIM] );
int GetBin (
    const double value,
    const double BinWidth,
    const double BinStart );
void OutputCorrelation_Header(
    ofstream * const file );
void OutputCorrelation_Value(
    ofstream * const file,
    BinPairMap::const_iterator &runner );
void OutputDipoleAngularCorrelation_Header(
    ofstream * const file );
void OutputDipoleAngularCorrelation_Value(
    ofstream * const file,
    DipoleAngularCorrelationMap::const_iterator &runner );
void OutputDipoleCorrelation_Header(
    ofstream * const file );
void OutputDipoleCorrelation_Value(
    ofstream * const file,
    DipoleCorrelationMap::const_iterator &runner );
void OutputPairCorrelation_Header(
    ofstream * const file );
void OutputPairCorrelation_Value(
    ofstream * const file,
    PairCorrelationMap::const_iterator &runner );
void OutputCorrelationToPoint_Header(
    ofstream * const file );
void OutputCorrelationToPoint_Value(
    ofstream * const file,
    CorrelationToPointMap::const_iterator &runner );
void OutputCorrelationToSurface_Header(
    ofstream * const file );
void OutputCorrelationToSurface_Value(
    ofstream * const file,
    CorrelationToSurfaceMap::const_iterator &runner );

/** Searches for lowest and highest value in a given map of doubles.
 * \param *map map of doubles to scan
 * \param &min minimum on return
 * \param &max maximum on return
 */
template <typename T> void GetMinMax ( T *map, double &min, double &max)
{
  min = 0.;
  max = 0.;
  bool FirstMinFound = false;
  bool FirstMaxFound = false;

  if (map == NULL) {
    ELOG(0, "Nothing to min/max, map is NULL!");
    performCriticalExit();
    return;
  }

  for (typename T::iterator runner = map->begin(); runner != map->end(); ++runner) {
    if ((min > runner->first) || (!FirstMinFound)) {
      min = runner->first;
      FirstMinFound = true;
    }
    if ((max < runner->first) || (!FirstMaxFound)) {
      max = runner->first;
      FirstMaxFound = true;
    }
  }
};

/** Puts given correlation data into bins of given size (histogramming).
 * Note that BinStart and BinEnd are desired to fill the complete range, even where Bins are zero. If this is
 * not desired, give equal BinStart and BinEnd and the map will contain only Bins where the count is one or greater. If a
 * certain start value is desired, give BinStart and a BinEnd that is smaller than BinStart, then only BinEnd will be
 * calculated automatically, i.e. BinStart = 0. and BinEnd = -1. .
 * Also note that the range is given inclusive, i.e. [ BinStart, BinEnd ].
 * \param *map map of doubles to count
 * \param BinWidth desired width of the binds
 * \param BinStart first bin
 * \param BinEnd last bin
 * \return Map of doubles (the bins) with counts as values
 */
template <typename T> BinPairMap *BinData ( T *map, const double BinWidth, const double BinStart, const double BinEnd )
{
  BinPairMap *outmap = new BinPairMap;
  int bin = 0;
  double start = 0.;
  double end = 0.;
  pair <BinPairMap::iterator, bool > BinPairMapInserter;

  if (map == NULL) {
    ELOG(0, "Nothing to bin, is NULL!");
    performCriticalExit();
    return outmap;
  }

  if (BinStart == BinEnd) { // if same, find range ourselves
    GetMinMax( map, start, end);
  } else if (BinEnd < BinStart) { // if BinEnd smaller, just look for new max
    GetMinMax( map, start, end);
    start = BinStart;
  } else { // else take both values
    start = BinStart;
    end = BinEnd;
  }
  // limit precision of start and end bin to 6 digits
  const double precision = 1e-6/BinWidth;
  end = precision*floor(end/precision);
  start = precision*floor(start/precision);
  for (int runner = 0; runner <= ceil((end-start)/BinWidth); runner++)
    outmap->insert( pair<double, int> ((double)runner*BinWidth+start, 0) );

  for (typename T::iterator runner = map->begin(); runner != map->end(); ++runner) {
    bin = GetBin (runner->first, BinWidth, start);
    BinPairMapInserter = outmap->insert ( pair<double, int> ((double)bin*BinWidth+start, 1) );
    if (!BinPairMapInserter.second) {  // bin already present, increase
      BinPairMapInserter.first->second += 1;
    } else { // bin not present, then remove again
      outmap->erase(BinPairMapInserter.first);
    }
  }

  return outmap;
};

/** Prints correlation map of template type to file.
 * \param *file file to write to
 * \param *map map to write
 * \param (*header) pointer to function that adds specific part to header
 * \param (*value) pointer to function that prints value from given iterator
 */
template <typename T>
void OutputCorrelationMap(
    ofstream * const file,
    const T * const map,
    void (*header)( ofstream * const file),
    void (*value)( ofstream * const file, typename T::const_iterator &runner)
    )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tBinCenter\tBinEnd";
  header(file);
  *file << endl;
  typename T::const_iterator advancer = map->begin();
  typename T::const_iterator runner = map->begin();
  if (advancer != map->end())
    advancer++;
  for ( ;(advancer != map->end()) && (runner != map->end()); ++advancer, ++runner) {
    *file << setprecision(8)
        << runner->first << "\t"
        << (advancer->first + runner->first)/2. << "\t"
        << advancer->first << "\t";
    value(file, runner);
    *file << endl;
  }
  if(runner != map->end()) {
    *file << setprecision(8) << runner->first << "\t0\t0\t";
    value(file, runner);
    *file << endl;
  }
};


#endif /* ANALYSIS_HPP_ */