TestVectorRealizer_gsl.C

Go to the documentation of this file.

#define BOOST_TEST_MODULE
#include <boost/test/included/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>

#include <queso/VectorSpace.h>
#include <queso/VectorRealizer.h>
#include <queso/GslMatrix.h>

using namespace QUESO;

BOOST_AUTO_TEST_CASE( test_uqGaussianVectorRealizer )
{
  // Initialize
  MPI_Init(NULL, NULL);
  uqFullEnvironment env; // Puts random number generator in known state... will generate same sequence each run
  uqVectorSpace<uqGslVector, uqGslMatrix> imageSpace(env, "test_space", 2, NULL);
  uqMap eMap(2, 0, env.comm());

  double tol = 1e-16;


  // Tests

  // Test 1: mean = 0, covMatrix = identity
  uqGslVector expectedValues(env, eMap, 0.0);
  uqGslMatrix lowerCholCovMatrix(env, eMap, 1.0); // identity

  int ierr = lowerCholCovMatrix.chol();
  BOOST_REQUIRE( ierr==0 ); // make sure cholesky succeeded

  lowerCholCovMatrix.zeroUpper(false); // zero upper triangular 

  uqGaussianVectorRealizer<uqGslVector, uqGslMatrix>* gaussianRealizer = 
    new uqGaussianVectorRealizer<uqGslVector, uqGslMatrix>("test_realizer", 
                                                                          imageSpace, 
                                                                          expectedValues, 
                                                                          lowerCholCovMatrix);

  uqGslVector myRealization(expectedValues);

  gaussianRealizer->realization(myRealization);
  //std::cout << "myRealization = " << std::setprecision(16) << std::scientific << myRealization << "\n";

  // NOTE: Data generated by running rng (i.e. this is a regression test).
  // If default seed value or rng algorithm are changed, these asserts will fail.
  BOOST_REQUIRE_CLOSE(myRealization[0],  2.2285703126720258e-01, tol);
  BOOST_REQUIRE_CLOSE(myRealization[1], -8.5156524864967331e-01, tol);

  delete gaussianRealizer;
  
  // Test 2: mean = [-2; 1], covMatrix = [2, 1; 1, 2]
  expectedValues[0] = -2.0; expectedValues[1] = 1.0;
  
  lowerCholCovMatrix(0,0) = 2.0; lowerCholCovMatrix(0,1) = 1.0;
  lowerCholCovMatrix(1,0) = 1.0; lowerCholCovMatrix(1,1) = 2.0;

  ierr = lowerCholCovMatrix.chol();
  BOOST_REQUIRE( ierr==0 ); // make sure cholesky succeeded

  lowerCholCovMatrix.zeroUpper(false); // zero upper triangular 

  gaussianRealizer = new uqGaussianVectorRealizer<uqGslVector, uqGslMatrix>("test_realizer", 
                                                                                           imageSpace, 
                                                                                           expectedValues, 
                                                                                           lowerCholCovMatrix);

  // Generate realization
  gaussianRealizer->realization(myRealization);
  //std::cout << "myRealization = " << std::setprecision(16) << std::scientific << myRealization << "\n";

  BOOST_REQUIRE_CLOSE(myRealization[0], -1.315078621127142e+00, tol);
  BOOST_REQUIRE_CLOSE(myRealization[1],  9.780380444774379e-01, tol);

  delete gaussianRealizer;


  // Clean up
  MPI_Finalize();

}