//
// 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 "vislearning/nice_nonvis.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)
        srand48(time(NULL));

    // 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);
     
}