/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: GridFeatureParser.cxx,v $
  Language:  C++
  Date:      $Date: 2003/05/23 21:21:32 $
  Version:   $Revision: 1.15 $

  Copyright (c) 2002 Insight Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.


=========================================================================*/

#include <iostream>
#include <memory>

#include <fstream>
#include "itkImageFileReader.h"
#include "vnl/vnl_matrix.h"

#include "itkKullbackLeiblerNMF.h"


//int KLNMFTest(int, char* [] )
int main( int argc, char *argv[] )
{
  //Define test matrix
  typedef unsigned int InputMatrixValueType;
  unsigned int nrows = 4;
  unsigned int ncols = 6;
  vnl_matrix< InputMatrixValueType > mat(nrows,ncols);


  mat[0][0] =1; mat[0][1] =1; mat[0][2] =1; mat[0][3] =2; mat[0][4] =2; mat[0][5] =2;
  mat[1][0] =4; mat[1][1] =4; mat[1][2] =4; mat[1][3] =3; mat[1][4] =3; mat[1][5] =3;
  mat[2][0] =4; mat[2][1] =4; mat[2][2] =4; mat[2][3] =3; mat[2][4] =3; mat[2][5] =3;
  mat[3][0] =1; mat[3][1] =1; mat[3][2] =1; mat[3][3] =2; mat[3][4] =2; mat[3][5] =2;

  itk::KullbackLeiblerNMF< InputMatrixValueType > nmf;

  nmf.SetInput( mat );
  nmf.SetNumberOfClasses(2);
  nmf.SetVerbose(false);
  nmf.Initialize();
  nmf.Compute();

  //nmf.GetHMatrix()->print(std::cout);
  vnl_vector<unsigned int> testH = *nmf.GetMLELabel( nmf.GetHMatrix(), 2);

  vnl_vector<unsigned int> gtH1(ncols), gtH2(ncols);
  gtH1.fill(1); gtH1[3]=2; gtH1[4]=2; gtH1[5]=2;
  gtH2.fill(2); gtH2[3]=1; gtH2[4]=1; gtH2[5]=1;
 
  
  if ( ((testH-gtH1).sum() != 0) || ((testH-gtH2).sum()!=0 ) )
    std::cout << "Test failed." << std::endl;
  else
    std::cout << "Test passed." << std::endl;

  return 0;
}