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

// STL headers
#include <map>

#include "defs.hpp"

#include "atom.hpp"

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

class BoundaryTriangleSet;
class element;
class LinkedCell;
class MoleculeListClass;
class Tesselation;
class Vector;

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

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

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

PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const  type2 );
CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point );
CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC );
PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const  type2, const int ranges[NDIM] );
CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] );
CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] );
int GetBin ( const double value, const double BinWidth, const double BinStart );
void OutputCorrelation( ofstream * const file, const BinPairMap * const map );
void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map );
void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map );
void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map );


/** 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) {
    DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to min/max, map is NULL!" << endl);
    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) {
    DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to bin, is NULL!" << endl);
    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;
  }
  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;
    }
  }

  return outmap;
};

#endif /* ANALYSIS_HPP_ */