source: src/Analysis/analysis_correlation.hpp@ beadd0

/*
 * 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<atom *> &atoms);
std::map<atomId_t, Vector> CalculateZeroAngularDipole(const std::vector<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(std::vector<molecule *> &molecules);
PairCorrelationMap *PairCorrelation(
    const World::AtomComposite &atoms_first,
    const World::AtomComposite &atoms_second,
    const double max_distance);
CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point );
CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC );
CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point, const int ranges[NDIM] );
CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<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_ */