//
// C++ Implementation: TestPDF
//
// Description:
//
//
// Author: Michael Koch <Koch.Michael@uni-jena.de>, (C) 2009
//
// Copyright: See COPYING file that comes with this distribution
//
/**
 * @file TestPDF.cpp
 * @brief TestPDF
 * @author Michael Koch
 * @date Di Aug 4 2009
 */

#include "TestPDF.h"
#include <string>

#include "core/basics/cppunitex.h"
#include "core/basics/numerictools.h"
#include "core/vector/Distance.h"

#include "core/vector/VVector.h"
#include "vislearning/math/pdf/PDFGaussian.h"

#include <time.h>

using namespace std;
using namespace NICE;
using namespace OBJREC;

Matrix computeCovariance ( const VVector & vecs )
{
  Vector mean ( vecs[0].size(), 0.0 );

  for ( unsigned int i = 0; i < vecs.size(); ++i )
    mean += vecs[i];

  mean *= 1.0 / vecs.size();

  Matrix cov ( vecs[0].size(), vecs[0].size(), 0.0 );

  for ( unsigned int i = 0; i < vecs.size(); ++i )
  {
    Vector diff = vecs[i] - mean;
    cov.addTensorProduct ( 1.0, diff, diff );
  }

  cov *= 1.0 / ( vecs.size() );

  return cov;
}

CPPUNIT_TEST_SUITE_REGISTRATION ( TestPDF );



void TestPDF::setUp()
{
}

void TestPDF::tearDown()
{
}

void TestPDF::TestPDFComputation()
{
  uint dim = 3;
  bool init_random = true ;
  uint samples = dim * 500;

  NICE::Matrix C ( dim, dim );
  NICE::Vector mean ( dim, 0.0 );
  VVector samplevectors;
  C.set ( 0.0 );
  cerr << "Sample from Gaussian" << endl;
  //init random
  if ( init_random )
#ifdef WIN32
	srand ( time ( NULL ) );
#else
    srand48 ( time ( NULL ) );
#endif
    

  // generate random symmetric matrix
  for ( uint i = 0 ; i < dim ; i++ )
    for ( uint j = i ; j < dim ; j++ )
    {
      C ( i, j ) = drand48();
      C ( j, i ) = C ( i, j );
    }
  C = C * C.transpose();
  cerr << "Ground-Truth covariance" << endl;
  cerr << C << endl;
  //initialize GaussPDF
  PDFGaussian pdf_gauss ( C, mean );

  //draw samples
  pdf_gauss.sample ( samplevectors, samples );
  Matrix Cov_test = computeCovariance ( samplevectors );
  cerr << "Estimated covariance" << endl;
  cerr << Cov_test << endl;

  NICE::Matrix diff = C - Cov_test;
  double frobNorm = diff.frobeniusNorm();

  cerr << "Frobenius norm: " << frobNorm << endl;
  CPPUNIT_ASSERT_DOUBLES_EQUAL_NOT_NAN ( 0.0, frobNorm, 0.1 );

}