/**
 * @file  Histogram.h
 * @brief Management of multidimensional histogram.
 *
 * Copyright (c) 2008, 2009, 2012, 2013 University of Pennsylvania.
 *
 * This file is part of DTI-DROID.
 *
 * DTI-DROID 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 3 of the License, or
 * (at your option) any later version.
 *
 * DTI-DROID is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with DTI-DROID.  If not, see <http://www.gnu.org/licenses/>.
 *
 * Contact: SBIA Group <sbia-software at uphs.upenn.edu>
 */

#ifndef _DTIDROID_HISTOGRM_H
#define _DTIDROID_HISTOGRM_H  


#include <iostream>
#include <map>
#include <vector>
#include <algorithm>
#include <utility>


/***
  * class: Histogram template, with specified data type and dimension.
***/
template < class T, int n>
class Histogram: public std::map< int, double>
{
public:
  ///@{
  /// default constructor
  Histogram() : _bRemoveMomentDC(false) 
  {_range[0]=0; _range[1]=255; _nBins.clear(); _nBins.resize(n,1); init();}
  /// constructor with initialization of one bin number for all dimensions
  Histogram(int nBin);
  /// constructor with initialization of a set of bin numbers for each dimension
  Histogram(const std::vector<int> &nBins);
  /// copy constructor
  Histogram(const Histogram<T,n> & rhs);
  
  /// destructor
  ~Histogram() {}
  ///@}
  
  /***
    * function: build the histogram with the given data
    * input: data is in vector format. outer vector denotes the dimension with size n.
    * return: true if build successfully, false otherwise.
    ***/
  bool build(const std::vector< std::vector<T> > & data);
  
  /// compute the entropy of this histogram
  double entropy();
  
  /// compute geometric moments
  /// input is the desired number of order for each dimension
  double moments(const std::vector<int> & orders);
  
  /// set the flag to remove moment DC or not
  void setRemoveMomentDC(bool bRemoveMomentDC) 
  { 
    _bRemoveMomentDC = bRemoveMomentDC; 
  }
  
  /// find the histogram value for a specified Id
  double find(const std::vector<int> & histID);
  
  /// clear all data for each entry of the histogram
  void clear();
  
  /// print out histogram
  void print();
  
  /// clear all histogram data and resize the histogram
  void resize(const std::vector<int> & nBins);
  
  /// set histogram value range
  void setValueRange(T valLow, T valHigh)
  {
    _range[0] = valLow; _range[1] = valHigh;
  }
  
  /// get histogram value range
  void getValueRange(T *valLow, T *valHigh)
  {
    *valLow = _range[0]; *valHigh = _range[1];
  }
  
  /// get number of bins
  const std::vector<int>& getBins() const {return _nBins;}
  const std::vector<int>& getBins() {return _nBins;}
  
private:
  /// number of bins
  std::vector<int> _nBins;
  
  /// value range
  T _range[2];
  
  /// initialize the histogram
  void init();
  
  /// transform a vector index to an integer index
  int vect2idx(const std::vector<int>& histID);
  
  /// transform an integer index to a vector index
  void idx2vect(int idx, std::vector<int>& histID);
  
  /// remove the DC value when computing geometric moment
  bool _bRemoveMomentDC;
};

/***
  * implementation
***/
template< class T, int n>
Histogram<T,n>::Histogram(int nBin)
{
  _nBins.resize(n);
  for (int i=0; i<n; i++) _nBins[i] = nBin;
  
  init();
}

template< class T, int n>
Histogram<T,n>::Histogram( const std::vector<int> &nBins)
{
  _nBins.resize(n);
  
  int n1 = nBins.size();
  n1 = std::min(n,n1);
  for (int i = 0; i < n1; i++) _nBins[i] = nBins[i];
  if (n1 < n)
  {
    for (int i=n1; i<n; i++) _nBins[i] = nBins[n1-1];
  }
  
  init();
}

template< class T, int n>
Histogram<T,n>::Histogram( const Histogram<T,n> & rhs)
{
  _nBins = rhs.getBins();
  insert(rhs.begin(), rhs.end());
}

template< class T, int n>
bool Histogram<T,n>::build(const std::vector< std::vector<T> > & data)
{
  if (data.size() < n)
  {
    std::cout << "ERROR: there are not enough data to build the histogram!" << std::endl;
    return false;
  }

  int dLen = data[0].size();
  
  for (int i=1; i<n; i++)
  {
    if (dLen != data[i].size())
    {
      std::cout <<"ERROR: the data length should be the same for each dimension!" << std::endl;
      return false;
    }
  }
  
  //vector<T> dMax(n, (T)0.0);
  //vector<T> dMin(n, (T)0.0);
  
  //for (int i=0; i<n; i++)
  //{
  //  dMax[i] = *max_element(data[i].begin(), data[i].end());
  //  dMin[i] = *min_element(data[i].begin(), data[i].end());
    
    //cout << "Max: " << (int)dMax[i] << "  Min: " << (int)dMin[i] << endl;
  //}
  
  // initialize the histogram  
  init();
  
  for (int m=0; m<dLen; m++)
  {
    // build the histogram
    // get the index for each data point
    std::vector<int> histId(n, -1);
    for (int i=0; i<n; i++)
    {
      if (_range[1] <= _range[0]) histId[i] = 0;
      else
      {
        if (data[i][m] > _range[1]) histId[i] = _nBins[i] - 1;
        else if (data[i][m] < _range[0]) histId[i] = 0;
        else 
        {
          histId[i] = (int)((double)_nBins[i]*(data[i][m] - _range[0])/(_range[1] - _range[0]));
          if (histId[i] == _nBins[i]) histId[i] = _nBins[i] -1;
        }
        //histId[i] = (int)((double)_nBins[i]*(data[i][m] - dMin[i])/(dMax[i] - dMin[i]));
        //if (histId[i] == _nBins[i]) histId[i] = _nBins[i] - 1;
      }
    }
    
    int idx = vect2idx(histId);   
    
    std::map<int,double>::iterator itr = std::map<int,double>::find(idx);
    if (itr == end())
    {
      std::cout << "ERROR: error during buiding the histogram. program stopped." << std::endl;
      return false;
    }
    itr->second++;
  }
  
  // normalize the histogram
  std::map<int,double>::iterator itr;
  for (itr=begin(); itr!=end(); itr++)
  {
    itr->second = itr->second / (double)dLen;
  }
    
  return true;
}

template< class T, int n>
double Histogram<T,n>::entropy()
{
  double en = 0.0;
  std::map<int,double>::iterator itr;
  for (itr=begin(); itr!=end(); itr++)
  {
    en += -itr->second*log(itr->second+1e-20)/log(2.0);
  }
  
  return en;
}

template< class T, int n>
double Histogram<T,n>::moments(const std::vector<int> & orders)
{
  if (orders.size() != n )
  {
    std::cout << "ERROR: the number of specified order for the moments is incorrect!" << std::endl;
    return 0.0;
  }
  double val = 0.0;
  
  int idx = 0;
  //if (_bRemoveMomentDC) idx = 1;
  
  std::map<int,double>::iterator itr=begin();
  //if (_bRemoveMomentDC) itr++;
  
  for (; itr!=end(); itr++, idx++)
  {
    std::vector<int> histID;
    idx2vect(idx, histID);
    
    double tval = 1.0;
    for (int i=0; i<n; i++)
      tval *= pow((double)(histID[i]+1), (double)orders[i]);
      
    val += tval*itr->second;
  }
  
  //if (_bRemoveMomentDC) 
  //{
  //  itr = begin();
  //  val /= (1.0 - itr->second);
  //}
  
  return val;
}

template< class T, int n>
double Histogram<T,n>::find(const std::vector<int> & histID)
{
  if (histID.size() != n)
  {
    std::cout << "ERROR: the index of the histogram is incorrect!" << std::endl;
    return 0.0;
  }
  
  for (int i=0; i<n; i++)
  {
    if (histID[i] >= _nBins[i])
    {
      std::cout << "ERROR: the histogram index is out of range!" << std::endl;
      return 0.0;
    }
  }
  
  int idx = vect2idx(histID);
  
  std::map<int, double>::iterator itr = find(idx);
  if (itr !=end())
    return(itr->second);
  else
    return(0.0);
}

template< class T, int n>
void Histogram<T,n>::clear()
{
 std::map<int, double>::iterator itr;
 for (itr=begin(); itr!=end(); itr++)
 {
   itr->second = 0.0;
 }
}

template< class T, int n>
void Histogram<T,n>::print()
{
  std::map<int, double>::iterator itr;
  int i = 0;
  
  for (itr=begin(); itr!=end(); i++,itr++)
  {
    std::cout << "[" << i << "]: " << itr->second << std::endl;
  }
  std::cout << std::endl;
}

template< class T, int n>
void Histogram<T,n>::resize(const std::vector<int> & nBins)
{
  std::map<int,double>::clear();
  
  _nBins.resize(n);
  
  int n1 = nBins.size();
  n1 = std::min(n,n1);
  for (int i = 0; i < n1; i++) _nBins[i] = nBins[i];
  if (n1 < n)
  {
    for (int i=n1; i<n; i++) _nBins[i] = nBins[n1-1];
  }
  
  init();
}

template< class T, int n>
void Histogram<T,n>::init()
{
  int hSize = 1;
  for (int i=0; i<n; i++)
    hSize = hSize*_nBins[i];
    
  for (int i=0; i<hSize; i++)
    insert(std::make_pair(i, 0.0));
}

template< class T, int n>
int Histogram<T,n>::vect2idx(const std::vector<int> & histID)
{
  int idx = 0;
  for (int i=0; i<n; i++)
  {
    int tId = 1;
    for (int j=0; j<i; j++)
      tId = tId*_nBins[j];
      
    idx += histID[i]*tId;
  }
  
  return idx;
}

template< class T, int n>
void Histogram<T,n>::idx2vect(int idx, std::vector<int>& histID)
{
  histID.clear();
  histID.resize(n);
  
  int tId = idx;
  for (int i=0; i<n; i++)
  {
    histID[i] = tId % _nBins[i];
    tId = tId / _nBins[i];
  }
}


#endif // _DTIDROID_HISTOGRM_H
