/*******************************************************************

  Class for treating multi-component system (including binary)

  Written by k-fleak

*******************************************************************/


#include <vector>
#include <deque>
#include <iostream>
#include <iomanip>
#include <limits>
#include <set>
#include <map>
#include <iterator>
#include <math.h>

#include "atom.cpp"

#ifndef __MULTICOMPONENT_CPP
#define __MULTICOMPONENT_CPP


typedef std::map<std::vector<int>, std::vector<double> > SigmaSet;
typedef std::set<std::set<int> > clusterIndexSet;
typedef std::set<int> clusterIndex;
typedef std::vector<std::vector<double> > FuncSet;
typedef std::vector<double> Func;
typedef std::vector<std::vector<std::vector<int> > > Combi;
typedef std::vector<std::set<std::vector<int> > > SortCombi;



class MultiComponent{
private:
  FuncSet Psi; 
  FuncSet Sigma;
  Combi   Vindex;
  SortCombi SortVindex;
  double numericLimits;
  std::vector<double> spinSet;
  void setNumericLimits();
  template<typename Spin> void setSpin(std::vector<Spin> &spinSetIn);
  void setPsi0();
  void setSigma();
  void setVindex(int MaxBody);
  void setSortVindex();
  void setPsi(int n);
  SigmaSet setMultiCorr(int Nbody);
  template<class Cont1> void printData(Cont1 &data);
  void printMultiCorr(SigmaSet &MultiIn);

protected:
  template<class Dat, class Spin> Dat getNorm(std::vector<Dat> &V1, std::vector<Dat> &V2, std::vector<Spin> &spinSetIn);
  template<class Vec, class Coef> Vec getProduct(Coef coeff, Vec &V1);
  std::vector<int> convertDecimal(int base, int digit, int inputNum);
  
public:
  template<typename Spin> explicit MultiComponent(int MaxBody, std::vector<Spin> &spinSetIn);
  ~MultiComponent();
  SigmaSet getMultiCorrelation(clusterIndexSet &clusterPosition, std::vector<Atom> &atomArray);
  FuncSet getPsi();
  SortCombi getSortVindex();

  void printPsi();
  void printVindex();
  void printSortVindex();

};


template<typename Spin> MultiComponent::MultiComponent(int MaxBody,  std::vector<Spin> &spinSetIn){
  setNumericLimits();
  setSpin(spinSetIn);
  setPsi0();
  setSigma();
  setVindex(MaxBody);
  setSortVindex();
  for ( int i = 1; i < spinSet.size(); ++i){
    setPsi(i);
  }
}


MultiComponent::~MultiComponent(){
  Psi.clear();
}


SigmaSet MultiComponent::getMultiCorrelation(clusterIndexSet &clusterPosition, std::vector<Atom> &atomArray){
  SigmaSet MultiCorr;
  clusterIndexSet::iterator it;
  it = clusterPosition.begin();
  MultiCorr = setMultiCorr( (*it).size() );

  while ( it != clusterPosition.end()){
    std::vector<std::vector<int> >::iterator itV;
    itV = Vindex[(*it).size() - 1].begin();
    while ( itV != Vindex[(*it).size() - 1].end() ){
      clusterIndex::iterator it2;
      std::vector<int>::iterator itV2;
      double tcorr = 1.0;
      it2 = (*it).begin();
      itV2 = (*itV).begin();
      while ( it2 != (*it).end()) {
        int spin = atomArray[(*it2)].getSpin();
        double tcorr2 = 0.0;
        int i = 0;
        Func::iterator itPsi;
        itPsi = Psi[*itV2].begin();
        while ( itPsi != Psi[*itV2].end() ){
          tcorr2 += (*itPsi) * pow(spin, i);
          ++itPsi;
          ++i;
        }
        tcorr *= tcorr2;
        ++it2;
        ++itV2;
      }
      std::vector<int> vtmp(*itV);
      sort (vtmp.begin(), vtmp.end());
      MultiCorr[vtmp][0] += tcorr;
      MultiCorr[vtmp][1] += 1.0;
      ++itV;
    }
    ++it;
  }
  return MultiCorr;
}


SigmaSet MultiComponent::setMultiCorr( int Nbody ){
  SigmaSet Multi;
  std::set<std::vector<int> >::iterator it;
  it = SortVindex[ Nbody -1 ].begin();
  std::vector<double> vtmp(2, 0.0);
  while ( it != SortVindex[ Nbody -1 ].end()) {
    Multi.insert( std::pair<std::vector<int>, std::vector<double> >(*it, vtmp) );
    ++it;
  }
  return Multi;
}


inline FuncSet MultiComponent::getPsi(){
  return Psi;
}

inline SortCombi MultiComponent::getSortVindex(){
  return SortVindex;
}


void MultiComponent::setNumericLimits(){
  std::numeric_limits<double> limits;
  double edigit = pow(10.0, 5.0);
  numericLimits = limits.epsilon() * edigit;
}


template<typename Spin> void MultiComponent::setSpin(std::vector<Spin> &spinSetIn){
  for ( int i = 0; i < spinSetIn.size(); ++i){
    spinSet.push_back(spinSetIn[i]);
  }
}


void MultiComponent::setPsi0(){
  std::vector<double> tmpPsi;
  tmpPsi.resize(spinSet.size());
  for (int i = 0; i < tmpPsi.size(); ++i){
    tmpPsi[i] = 0.0;
  }
  tmpPsi[0] = 1.0;
  Psi.push_back(tmpPsi);
  tmpPsi.clear();
}


void MultiComponent::setSigma(){
  std::vector<double> tmpSigma;
  tmpSigma.resize(spinSet.size());
  for (int i = 0; i < tmpSigma.size(); ++i){
    tmpSigma[i] = 0.0;
  }
  for (int i = 0; i < spinSet.size(); ++i){
    Sigma.push_back(tmpSigma);
  }
  for (int j = 0; j < Sigma.size(); ++j){
    Sigma[j][j] = 1.0;
  }
}


void MultiComponent::setVindex(int MaxBody){
  for ( int i = 1; i <= MaxBody; ++i ){
    int numIndex = static_cast<int> (pow( (spinSet.size()-1), i ));
    std::vector<std::vector<int> > vtmp;
    for ( int j = 0; j < numIndex; ++j){
      std::vector<int> vtmp2;
      vtmp2 = convertDecimal( (spinSet.size() - 1), i, j);
      vtmp.push_back(vtmp2);
    }
    Vindex.push_back(vtmp);
  }
}


void MultiComponent::setSortVindex(){
  Combi::iterator it;
  it = Vindex.begin();
  while ( it != Vindex.end()){
    std::vector<std::vector<int> >::iterator it2;
    it2 = (*it).begin();
    std::set<std::vector<int> > stmp;
    while ( it2 != (*it).end()){
      std::vector<int> vtmp(*it2);
      sort(vtmp.begin(), vtmp.end());
      stmp.insert(vtmp);
      ++it2;
    }
    SortVindex.push_back(stmp);
    ++it;
  }
}
  


std::vector<int> MultiComponent::convertDecimal(int base, int digit, int inputNum){
  std::deque<int> decNum;
  std::vector<int> decNumV;
  int stmp = inputNum;
  for ( int i =0; i < digit; ++i){
    int atmp;
    atmp = stmp % base;
    stmp /= base;
    decNum.push_front(atmp);
  }
  
  std::deque<int>::iterator it;
  it = decNum.begin();
  while ( it != decNum.end()){
    decNumV.push_back( (*it) + 1 );
    ++it;
  }
  return decNumV;
}
  

template<class Dat, class Spin> Dat MultiComponent::getNorm(std::vector<Dat> &V1, std::vector<Dat> &V2, std::vector<Spin> &spinSets){
  Dat norm = 0.0;
  for ( int k = 0; k < spinSets.size(); ++k){
    for ( int i = 0; i < V1.size(); ++i){
      for ( int j = 0; j < V2.size(); ++j){
	double tmp;
	tmp = V1[i] * pow( spinSets[k], i)  *  V2[j] * pow( spinSets[k], j);
	norm += tmp * (1.0 / spinSets.size());
      }
    }
  }
  if ( norm <= numericLimits ){
    norm = 0.0;
  }
  return norm;
}


template<class Vec, class Coef> Vec MultiComponent::getProduct(Coef coeff, Vec &V1){
  Vec V2;
  V2.resize(V1.size());
  for ( int i = 0; i < V2.size(); ++i ){
    V2[i] = V1[i] * coeff;
  }
  return V2;
}


void MultiComponent::setPsi(int n){
  std::vector<double> gn;
  gn.resize(spinSet.size());
  for ( int i = 0; i < gn.size(); ++i){
    gn[i] = 0.0;
  }
  gn[n] += 1;
  for ( int i = 0; i < n; ++i ){
    double norm;
    std::vector<double> ntmp;
    norm = getNorm(Psi[i], Sigma[n], spinSet);
    ntmp = getProduct(norm, Psi[i]);
    for ( int j = 0; j < gn.size(); ++j){
      gn[j] -= ntmp[j];
    }
  }
  
  double gnNorm;
  std::vector<double> Psin;
  gnNorm = sqrt(getNorm(gn, gn, spinSet));
  Psin = getProduct( (1.0 / gnNorm), gn);
  Psi.push_back(Psin);

}


template<class Cont1> void MultiComponent::printData(Cont1 &data){
  typename Cont1::iterator it1;
  typedef typename Cont1::value_type Cont2;
  typedef typename Cont2::value_type Cont3;
  it1 = data.begin();
  std::cout << std::endl << std::endl;
  while ( it1 != data.end()){
    typename Cont2::iterator it2;
    it2 = (*it1).begin();
    while ( it2 != (*it1).end()){
      typename Cont3::const_iterator it3 = (*it2).begin();
      while ( it3 != (*it2).end()){
	std::cout << *it3 << " ";
	++it3;
      }
      std::cout << std::endl;
      ++it2;
    }
    std::cout << std::endl;
    ++it1;
  }
}

template<> void MultiComponent::printData(FuncSet &data){
  FuncSet::iterator it;
  it = data.begin();
  std::cout << std::endl << std::endl;
  while ( it != data.end()){
    Func::iterator it2;
    it2 = (*it).begin();
    while ( it2 != (*it).end()){
      std::cout << *it2 << " ";
      ++it2;
    }
    std::cout << std::endl;
    ++it;
  }
}


template<> void MultiComponent::printData(SigmaSet &data){
  SigmaSet::iterator it;
  it = data.begin();
  while ( it != data.end()){
    std::vector<int>::const_iterator it2;
    it2 = (*it).first.begin();
    while ( it2 != (*it).first.end()){
      std::cout << *it2 << " ";
      ++it2;
    }
    std::cout << " : ";
    std::vector<double>::iterator it3;
    it3 = (*it).second.begin();
    while ( it3 != (*it).second.end()){
      std::cout << std::setw(14) << std::right << *it3 << " ";
      ++it3;
    }
    std::cout << std::endl;
    ++it;
  }
}

void MultiComponent::printVindex(){
  printData(Vindex);
}

void MultiComponent::printSortVindex(){
  printData(SortVindex);
}

void MultiComponent::printPsi(){
  printData(Psi);
}

void MultiComponent::printMultiCorr(SigmaSet &MultiIn){
  printData(MultiIn);
}


#endif