NICE_GP_HIK_Exp/progs/testLogDetApproximation.cpp

/** 
* @file testLogDetApproximation.cpp
* @brief 
* @author Alexander Freytag
* @date 05-01-2012 (dd-mm-yyyy)
*/

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <unistd.h>

#include "core/vector/MatrixT.h"
#include "core/vector/VectorT.h"

#include "core/algebra/GenericMatrix.h"
#include "core/algebra/EigValuesTRLAN.h"
#include "core/algebra/GMStandard.h"

#include "gp-hik-core/tools.h"

#include "gp-hik-core/algebra/LogDetApproxBaiAndGolub.h"


using namespace std;
using namespace NICE;

/**
 * @brief Printing main menu.
 * @author Alexander Freytag
 * @date 12/06/2011
 * 
 * @return void
 **/
void print_main_menu()
{
  std::cerr << std::endl << " - test program for logDet Approximation" << std::endl;
  std::cerr << "Input options:" << std::endl;
  std::cerr << "   -n <number>  dimension of K"<< std::endl;
  std::cerr << "   -v 1/0  verbose mode"<< std::endl;
  return;
}

int main (int argc, char* argv[])
{
  int n (5);
  bool verbose(false);
  
  int rc;
  if (argc<2)
  {
    print_main_menu();
    return -1;
  }
  
  while ((rc=getopt(argc,argv,"n:v:h"))>=0)
  {
    switch(rc)
    {
      case 'n': n = atoi(optarg); break;
      case 'v': verbose = atoi(optarg); break;
      default: print_main_menu();
    }
  }
  
  if (verbose)
  {
    std::cerr << "Testing logDet Approximation for n = " << n << std::endl;
  }
  
  srand ( time(NULL) );
  
  NICE::Matrix ARand(generateRandomMatrix(n,n));
  NICE::Matrix A;
  // A shall be positive definite
  A.multiply(ARand, ARand, true);
  
  NICE::GMStandard genericA(A);
  
  //compute GT LogDet based on eigenvalues of A
  NICE::Vector eigenvalues;
  NICE::Matrix eigenvectors;
  try
  {
    NICE::EigValuesTRLAN eigValuesComputation;
    eigValuesComputation.getEigenvalues(genericA, eigenvalues,eigenvectors, n );
  }
  catch (...)
  {
    NICE::EVArnoldi eigValuesComputation;
    eigValuesComputation.getEigenvalues(genericA, eigenvalues,eigenvectors, n );
  }
  
  
  double logDetGT(0.0);
  for (int i = 0; i < n; i++)
  {
    logDetGT += log(eigenvalues[i]);
  }
  
  if (verbose)
  {
    std::cerr << "GT logDet: " << logDetGT << std::endl;
  }
  
  //replace this later on using only the k largest eigenvalues
  double frobNorm(A.squaredFrobeniusNorm());
  
  NICE::LogDetApproxBaiAndGolub logDetApproximator;
  double logDetApprox(logDetApproximator.getLogDetApproximation(A.trace(), frobNorm, eigenvalues.Max(), eigenvalues.Min(), n ) );
  
  if (verbose)
  {
    std::cerr << "logDetApprox: " << logDetApprox << std::endl;
  }
  
  double logDetApproxUpperBound(logDetApproximator.getLogDetApproximationUpperBound(A.trace(), frobNorm, eigenvalues.Max(), n ) );
  
  if (verbose)
  {
    std::cerr << "logDetApproxUpperBound: " << logDetApproxUpperBound << std::endl;
  }
  
  double logDetApproxLowerBound(logDetApproximator.getLogDetApproximationUpperBound(A.trace(), frobNorm, eigenvalues.Min(), n ) );
  
  if (verbose)
  {
    std::cerr << "logDetApproxLowerBound: " << logDetApproxLowerBound << std::endl;
  }
  
  return 0;
}