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

#include "Helpers/MemDebug.hpp"

#include <iostream>

#include "analysis_correlation.hpp"
#include "element.hpp"
#include "info.hpp"
#include "log.hpp"
#include "molecule.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "triangleintersectionlist.hpp"
#include "vector.hpp"
#include "verbose.hpp"
#include "World.hpp"


/** Calculates the pair correlation between given elements.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate
 * \return Map of doubles with values the pair of the two atoms.
 */
PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements)
{
  Info FunctionInfo(__func__);
  PairCorrelationMap *outmap = NULL;
  double distance = 0.;
  double *domain = World::getInstance().getDomain();

  if (molecules->ListOfMolecules.empty()) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();

  // create all possible pairs of elements
  set <pair<element *, element *> > PairsOfElements;
  if (elements.size() >= 2) {
    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
        if (type1 != type2) {
          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
        }
  } else if (elements.size() == 1) { // one to all are valid
    element *elemental = *elements.begin();
    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
  } else { // all elements valid
    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
  }

  outmap = new PairCorrelationMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    if ((*MolWalker)->ActiveFlag) {
      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
        for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
          if ((*MolOtherWalker)->ActiveFlag) {
            DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
            for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
              DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
              if ((*iter)->getId() < (*runner)->getId()){
                for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                  if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
                    distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
                    //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
                  }
              }
            }
          }
        }
      }
    }
  }
  return outmap;
};

/** Calculates the pair correlation between given elements.
 * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate
 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
 * \return Map of doubles with values the pair of the two atoms.
 */
PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const int ranges[NDIM] )
{
  Info FunctionInfo(__func__);
  PairCorrelationMap *outmap = NULL;
  double distance = 0.;
  int n[NDIM];
  Vector checkX;
  Vector periodicX;
  int Othern[NDIM];
  Vector checkOtherX;
  Vector periodicOtherX;

  if (molecules->ListOfMolecules.empty()) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();

  // create all possible pairs of elements
  set <pair<element *, element *> > PairsOfElements;
  if (elements.size() >= 2) {
    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
        if (type1 != type2) {
          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
        }
  } else if (elements.size() == 1) { // one to all are valid
    element *elemental = *elements.begin();
    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
  } else { // all elements valid
    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
  }

  outmap = new PairCorrelationMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    if ((*MolWalker)->ActiveFlag) {
      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
      double * FullInverseMatrix = InverseMatrix(FullMatrix);
      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
        periodicX = *(*iter)->node;
        periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
        // go through every range in xyz and get distance
        for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
          for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
            for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
              checkX = Vector(n[0], n[1], n[2]) + periodicX;
              checkX.MatrixMultiplication(FullMatrix);
              for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
                if ((*MolOtherWalker)->ActiveFlag) {
                  DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
                  for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
                    DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
                    if ((*iter)->getId() < (*runner)->getId()){
                      for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                        if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
                          periodicOtherX = *(*runner)->node;
                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
                          // go through every range in xyz and get distance
                          for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
                            for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
                              for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
                                checkOtherX.MatrixMultiplication(FullMatrix);
                                distance = checkX.distance(checkOtherX);
                                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
                              }
                        }
                    }
                  }
                }
              }
            }
      }
      delete[](FullMatrix);
      delete[](FullInverseMatrix);
    }
  }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a point.
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate with point
 * \param *point vector to the correlation point
 * \return Map of dobules with values as pairs of atom and the vector
 */
CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point )
{
  Info FunctionInfo(__func__);
  CorrelationToPointMap *outmap = NULL;
  double distance = 0.;
  double *cell_size = World::getInstance().getDomain();

  if (molecules->ListOfMolecules.empty()) {
    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToPointMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
          if ((*type == NULL) || ((*iter)->type == *type)) {
            distance = (*iter)->node->PeriodicDistance(*point, cell_size);
            DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
            outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
          }
      }
    }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate to point
 * \param *point vector to the correlation point
 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
 * \return Map of dobules with values as pairs of atom and the vector
 */
CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
{
  Info FunctionInfo(__func__);
  CorrelationToPointMap *outmap = NULL;
  double distance = 0.;
  int n[NDIM];
  Vector periodicX;
  Vector checkX;

  if (molecules->ListOfMolecules.empty()) {
    DoLog(1) && (Log() << Verbose(1) <<"No molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToPointMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
      double * FullInverseMatrix = InverseMatrix(FullMatrix);
      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
          if ((*type == NULL) || ((*iter)->type == *type)) {
            periodicX = *(*iter)->node;
            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
            // go through every range in xyz and get distance
            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
                  checkX.MatrixMultiplication(FullMatrix);
                  distance = checkX.distance(*point);
                  DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
                  outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
                }
          }
      }
      delete[](FullMatrix);
      delete[](FullInverseMatrix);
    }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element and a surface.
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate to surface
 * \param *Surface pointer to Tesselation class surface
 * \param *LC LinkedCell structure to quickly find neighbouring atoms
 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
 */
CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
{
  Info FunctionInfo(__func__);
  CorrelationToSurfaceMap *outmap = NULL;
  double distance = 0;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;

  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    DoeLog(1) && (eLog()<< Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToSurfaceMap;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
      if ((*MolWalker)->empty())
        DoLog(1) && (1) && (Log() << Verbose(1) << "\t is empty." << endl);
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
          if ((*type == NULL) || ((*iter)->type == *type)) {
            TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
            distance = Intersections.GetSmallestDistance();
            triangle = Intersections.GetClosestTriangle();
            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
          }
      }
    } else {
      DoLog(1) && (Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl);
    }

  return outmap;
};

/** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
 * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
 * I.e. We multiply the atom::node with the inverse of the domain matrix, i.e. transform it to \f$[0,0^3\f$, then add per
 * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
 * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
 * \param *molecules list of molecules structure
 * \param &elements vector of elements to correlate to surface
 * \param *Surface pointer to Tesselation class surface
 * \param *LC LinkedCell structure to quickly find neighbouring atoms
 * \param ranges[NDIM] interval boundaries for the periodic images to scan also
 * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
 */
CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
{
  Info FunctionInfo(__func__);
  CorrelationToSurfaceMap *outmap = NULL;
  double distance = 0;
  class BoundaryTriangleSet *triangle = NULL;
  Vector centroid;
  int n[NDIM];
  Vector periodicX;
  Vector checkX;

  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    DoLog(1) && (Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl);
    return outmap;
  }
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    (*MolWalker)->doCountAtoms();
  outmap = new CorrelationToSurfaceMap;
  double ShortestDistance = 0.;
  BoundaryTriangleSet *ShortestTriangle = NULL;
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    if ((*MolWalker)->ActiveFlag) {
      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
      double * FullInverseMatrix = InverseMatrix(FullMatrix);
      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
          if ((*type == NULL) || ((*iter)->type == *type)) {
            periodicX = *(*iter)->node;
            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
            // go through every range in xyz and get distance
            ShortestDistance = -1.;
            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
                  checkX.MatrixMultiplication(FullMatrix);
                  TriangleIntersectionList Intersections(&checkX,Surface,LC);
                  distance = Intersections.GetSmallestDistance();
                  triangle = Intersections.GetClosestTriangle();
                  if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
                    ShortestDistance = distance;
                    ShortestTriangle = triangle;
                  }
                }
            // insert
            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
            //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
          }
      }
      delete[](FullMatrix);
      delete[](FullInverseMatrix);
    }

  return outmap;
};

/** Returns the index of the bin for a given value.
 * \param value value whose bin to look for
 * \param BinWidth width of bin
 * \param BinStart first bin
 */
int GetBin ( const double value, const double BinWidth, const double BinStart )
{
  Info FunctionInfo(__func__);
  int bin =(int) (floor((value - BinStart)/BinWidth));
  return (bin);
};


/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tCount" << endl;
  for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
  }
};

/** Prints correlation (double, (atom*,atom*) ) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tAtom1\tAtom2" << endl;
  for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
  }
};

/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tAtom::x[i]-point.x[i]" << endl;
  for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    *file << runner->first;
    for (int i=0;i<NDIM;i++)
      *file << "\t" << setprecision(8) << (runner->second.first->node->at(i) - runner->second.second->at(i));
    *file << endl;
  }
};

/** Prints correlation (double, int) pairs to file.
 * \param *file file to write to
 * \param *map map to write
 */
void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
{
  Info FunctionInfo(__func__);
  *file << "BinStart\tTriangle" << endl;
  if (!map->empty())
    for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
      *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
    }
};