#ifndef _itkKullbackLeiblerNMF_h
#define _itkKullbackLeiblerNMF_h


#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <exception>
#include <sstream>
#include <stdexcept>
#include <math.h>

#include "vnl/vnl_matrix.h"
#include "vnl/vnl_sample.h"
#include "vnl/vnl_vector.h"

#include "itkNMFBase.h"
namespace itk
{
/** \class KullbackLeiblerNMF
 * \brief Specific instantiation of the base class (NMFBase)for 
 * NMF based dimentionality reduction/clustering object using Kullback-Liebler
 * divergence as the metric.
 *
 * itkKullbackLeiblerNMF provides an implementation of non-negative matrix 
 * factorization (NMF) using Kullback-Liebler divergence as the underlying
 * metric. The implementation used the interface provided by the NMFBase
 * class. Two functions are implemented in this class Compute/NMFUpdate
 *
 * The Compute function first checks the validity of the input matrices
 * and the dimensions of the initialization of the decomposed matrices.
 * The routine then calls the NMFUpdate method iteratively and tests for
 * convergence after every few iterations set by the 
 * m_ConvergenceTestStepSize parameter. 
 *
 * The NMFUpdate function implements the equations that decompse the input
 * matrix V to W and H matrices. The equations implemented are derived from
 * the following publication: D. Lee and H. Seung 
 * ``Learning the parts of objects by non-negative matrix factorization''
 * {\em Nature}, vol. 40, pp. 788-791, 1999.
 *  
 * The NMF decompositions are interative in nature. This implementation 
 * provides a default convergence testing method. 
 * The implemented method checks for convergence in the H matrix (although W 
 * matrix could also be chosen for this purpose). The user specifies the 
 * maximum number of iterations (default set to 5000). The convergence is
 * tested every 10 iterations in the default implementation. The algorithm
 * deems the decomposition to have converged if the matrix does not change
 * after 10 * stopConv iterations. The stopConv parameter is set by the
 * user (default value is set to 5).
 *
 * \ingroup ImageFeatureExtraction */ 

template< class T >
class KullbackLeiblerNMF: public NMFBase<T>
{
public:
  typedef NMFBase<T> SuperClass;
  typedef vnl_vector<unsigned int> MLELabelVectorType;
  typedef typename SuperClass::InMatrixType InMatrixType;
  typedef typename SuperClass::OutMatrixType OutMatrixType;

  /** Number of iterations of the decomposition
    * default is set to 5000 */
  void SetNumberOfIterations( unsigned int numIter )
    {
    m_Iter = numIter;
    }

  /** Stop criterion for convergence
    * default is set to 5 */
  void SetStopConvergence( unsigned int stopConv )
    {
    m_Stop = stopConv;
    }

  /** Test e criterion for convergence
    * default is set to 10 */
  void SetConvergenceTestStepSize( unsigned int stepSize )
    {
    m_ConvergenceTestStepSize = stepSize;
    }

  /** Verbose flag to output results of intermediate 
    * convergence testing. Default is set to false */
  void SetVerbose( bool verbose )
    {
    m_Verbose = verbose;
    }

  /** The main function that calls the rest of the functions */
  virtual void Compute();

  /** The NMF update equations as stated in the referred manuscript
    * is implemented here */
  virtual void NMFUpdate();

  /** Default implementation to test convergence for matrix H */
  unsigned int TestConvergence();

  /** Generate the maximum liklihood label for the specified matrix;
    * 1 stands for W matrix; 2 stands for H matrix */
  MLELabelVectorType *GetMLELabel( OutMatrixType *m, unsigned int e);

  /** Write the MLE estimation results to a binary file */
  int WriteMLE(vnl_vector<unsigned int>  *m, std::string& str);

  KullbackLeiblerNMF();
  virtual ~KullbackLeiblerNMF() throw();

private:
  MLELabelVectorType         m_MLELabel;
  const double               m_Eps;
  const unsigned int         m_ConvergenceTestStepSize;

  unsigned int               m_Iter;
  unsigned int               m_Stop;
  bool                       m_Verbose;

  vnl_vector<int>            m_Label;
  vnl_vector<int>            m_LabelOld;


};//class itkKullbackLeiblerNMF


} // namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkKullbackLeiblerNMF.txx"
#endif

#endif