#include "freq.h"
#include <cmath>

using namespace std;

inline double entropy(const double p)
{
   if(p==0) return 0;
   double v = abs(p*log(p)/log(2.0)); 
   if(v<1e-10) v= 0;
   return v;
}
template <class T>
double Frequencies<T>::getCountInRange(const uint start, const uint end)
{ 
   if(m_fr.empty() || start>=end || end>m_fr.size()) return 0;
   uint leftVal = start ? m_fr[start-1] : 0;
   uint rightVal = m_fr[end-1];

   return (rightVal-leftVal+.0);
}

template <class T>
double Frequencies<T>::getFreqInRange(const uint start, const uint end)
{ 

   return getCountInRange(start, end)/(end-start);
}

template <class T>
void JensenShannonDivergence<T>::cleanup()
{
   for(uint j = 0; j<m_freq.size(); ++j){
      delete m_freq[j];
   }
}

template <class T>
bool JensenShannonDivergence<T>::calculate(const vector<T>& seq)
{
   cleanup();
   for(uint i = 0; i<seq.size(); ++i){
      T val = seq[i];
      uint idx = find(val);
      if(idx != uint(-1)){
         m_freq[idx]->increment(i);
      }
      else{
         Frequencies<T> *pFreq = new Frequencies<T>(val);
         if(!pFreq){
            cleanup(); 
            return false;
         }
         pFreq->increment(i);
         m_freq.push_back(pFreq);
      }
   }
   for(uint i = 0;i<m_freq.size(); ++i){
      m_freq[i]->increment(seq.size()-1, true);      
   }
   //dumpData<T>("Seq", seq, cout);      
   for(uint i = 0;i<m_freq.size(); ++i){
      //cout<<m_freq[i]->get()<<endl;
      m_freq[i]->dump();      
   }
   return true;
}

template <class T>
uint JensenShannonDivergence<T>::find(const T& val)
{
   for(uint j = 0; j<m_freq.size();++j){
      if(m_freq[j]->isEqual(val)){
         return j;
      }
   }
   return uint(-1);
}
template <class T>
double JensenShannonDivergence<T>::getEntropie(const uint idx, const uint start, const uint end)
{
   double p = m_freq[idx]->getFreqInRange(start, end);
   return p ? p*log(p)/log(2.0) : 0;
}

template <class T>
double JensenShannonDivergence<T>::getJointEntropie(
      const uint idx, 
      const uint start1, const uint end1,
      const uint start2, const uint end2
   )
{
   if(m_size==0 || idx>m_freq.size()-1) return 0;
   double w1 = (end1-start1 + .0) / m_size;
   double w2 = (end2-start2 + .0) / m_size;

   double p1 = m_freq[idx]->getFreqInRange(start1, end1);
   double p2 = m_freq[idx]->getFreqInRange(start2, end2);
   double p = (w1*p1+w2*p2);
   return entropy(p);   
}

template <class T>
double JensenShannonDivergence<T>::getAt(const uint pos)
{
   double ret = 0;
   double sumJoint = 0, sumEnt1 = 0, sumEnt2 = 0;
   double w1 = (pos+1.0) / m_size;
   double w2 = (m_size-pos-1.0) / m_size;
   
   for(uint i = 0; i<m_freq.size(); ++i){
      double p1 = m_freq[i]->getFreqInRange(0, pos+1);
      double p2 = m_freq[i]->getFreqInRange(pos+1, m_size);
      double p = w1*p1+w2*p2;
      ret += entropy(p)-w1*entropy(p1) -w2*entropy(p2);
   }
   ret = abs(ret);
   ret = ret<1e-10 ?
      0 : ret;
   return ret;
}

template <class T>
double JensenShannonDivergence<T>::getSignificance(const uint k, const uint N)
{
   
}

template <class T>
double JensenShannonDivergence<T>::LnGamma(double xx)
{
   double x = .0,y = .0,tmp = .0,ser = .0;
   static double cof[6]={
      76.18009172947146,
      -86.50532032941677,
      24.01409824083091,
      -1.231739572450155,
      0.1208650973866179e-2,
      -0.5395239384953e-5
   };
   if(xx<0) return -1;
   y=x=xx;
   tmp=x+5.5;
   tmp -= (x+0.5)*log(tmp);
   ser=1.000000000190015;
   for (uint j=0;j<=5;j++) ser += cof[j]/++y;
   return -tmp+log(2.5066282746310005*ser/x);
}

template <class T>
double JensenShannonDivergence<T>::IncGamma(double a, double x)
{
   double gamser,gammcf,gln;
   if (x < 0.0 || a <= 0.0) return -1;
   if (x < (a+1.0)) { 
      gser(&gamser,a,x,&gln);
      return gamser;
   } 
   else { 
      gcf(&gammcf,a,x,&gln);
      return 1.0-gammcf;   
   }
}

#define ITMAX 100 
#define EPS 3.0e-7
#define FPMIN 1.0e-30 

template <class T>
void JensenShannonDivergence<T>::gser(double *gamser, double a, double x, double *gln)
{
   int n = 0;
   double sum = 0,del = 0,ap = 0;
   *gln = LnGamma(a);
   if (x <= 0.0) {
      *gamser = 0.0;
      return;
   } 
   else {
      ap=a;
      del=sum=1.0/a;
      for (n=1;n<=ITMAX;n++) {
         ++ap;
         del *= x/ap;
         sum += del;
         if (fabs(del) < fabs(sum)*EPS) {
            *gamser=sum*exp(-x+a*log(x)-(*gln));
            return;
         }
      }
      return;
   }
}
template <class T>
void JensenShannonDivergence<T>::gcf(double *gammcf, double a, double x, double *gln)
{
   double an,b,c,d,del,h;
   *gln=LnGamma(a);
   b = x+1.0-a;
   c = 1.0/FPMIN;
   d = 1.0/b;
   h = d;
   for ( int i = 1;i <= ITMAX;i++) { 
      an = -i*(i-a);
      b += 2.0;
      d = an*d+b;
      if (fabs(d) < FPMIN) d=FPMIN;
      c = b + an/c;
      if (fabs(c) < FPMIN) c=FPMIN;
      d = 1.0/d;
      del = d*c;
      h *= del;
      if (fabs(del-1.0) < EPS) break;
   }
   *gammcf=exp(-x+a*log(x)-(*gln))*h; 
}