/***************************************************************************
                          Frequencies.h-  description
                             -------------------
    begin                : Mi Sept 13 2006
    copyright            : (C) 2006 by klaus frieler
    email                : kgf@omniversum.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#ifndef _FREQ_H
#define _FREQ_H
#include "Types.h"
#include <vector>

template <class T>
class Frequencies 
{
public:
   Frequencies(const T val):m_val(val){} 
   void  increment(const uint pos, const bool bFill = false)
   {
      if( m_fr.size()>pos) return;
      uint nCount = m_fr.size()? *m_fr.rbegin() : 0; 
      int nFill =  pos - m_fr.size();
      if(nFill>0) for(int i = 0; i<nFill; ++i)  m_fr.push_back(nCount);
      m_fr.push_back(bFill ? nCount : nCount+1);
   }
   uint getCountAt(const uint idx)
   { 
      return getCountInRange(idx, idx+1);
   }
   double getFreqAt(const uint idx)
   { 
      return getFreqInRange(idx, idx+1);
   }
   double getCountInRange(const uint start, const uint end);
   double getFreqInRange(const uint start, const uint end);

   double getEntropy(const uint idx){ 
      if(m_fr.empty()) return 0;
      double p = getFreqAt(idx);
      if(p== 0) return 0;
      return p*log(p)/log(2);
   }
   bool operator==(const T& v){return m_val == v;}
   bool isEqual(const T& v){return m_val == v;}
   const T& get(){ return m_val;}
   void dump()
   { 
      //dumpData<uint>("", m_fr, cout);
   }
protected:
   T m_val;
   std::vector<uint> m_fr;
};

template<class T>
class JensenShannonDivergence
{
public:
   JensenShannonDivergence(const std::vector<T>& seq){ calculate(seq); m_size = seq.size();}
   double getAt(const uint pos);
   uint size(){return m_size;}
   double Gamma(double x){ return exp(LnGamma(x));}
   double LnGamma(double xx);
   double IncGamma(double a, double x);
   double getSignificance(const uint k, const uint N);
protected:
   double getEntropie(const uint idx, const uint start, const uint end);
   double getJointEntropie(
      const uint idx, 
      const uint start1, const uint end1,
      const uint start2, const uint end2
   );
   bool calculate(const std::vector<T>& seq);
   vector<Frequencies<T>* > m_freq;
   void cleanup();
   uint find(const T& val);
   void gser(double *gamser, double a, double x, double *gln);
   void gcf(double *gammcf, double a, double x, double *gln);

   uint m_size;
};
#endif