testRandomSamplerFromCDF.cc

Go to the documentation of this file.

#include <tst/Verify.h>
#include <cppunit/extensions/HelperMacros.h>
#include <utl/RandomSamplerFromCDF.h>
#include <utl/AugerException.h>
#include <utl/TabulatedFunction.h>
#include <CLHEP/Random/Randomize.h>
#include <boost/lexical_cast.hpp>

using namespace std;
using namespace utl;
using namespace tst;
using boost::lexical_cast;


#define EQUAL(x, y) CPPUNIT_ASSERT(Verify<Equal>(x, y))
#define CLOSE(x, y) CPPUNIT_ASSERT(Verify<CloseTo>(x, y))
#define CLOSE_AT(x, y, eps) CPPUNIT_ASSERT(Verify<CloseTo>(x, y, eps))


class TestRandomSamplerFromCDF : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(TestRandomSamplerFromCDF);
  CPPUNIT_TEST_EXCEPTION(TestTooShort, InvalidConfigurationException);
  CPPUNIT_TEST_EXCEPTION(TestStartZero, InvalidConfigurationException);
  CPPUNIT_TEST_EXCEPTION(TestMonotonous, InvalidConfigurationException);
  CPPUNIT_TEST(TestTabulatedFunction);
  CPPUNIT_TEST(TestVectors);
  CPPUNIT_TEST(TestVectorAndNormalization);
  CPPUNIT_TEST(TestMap);
  CPPUNIT_TEST(TestFunction);
  CPPUNIT_TEST(TestFunction2);
  CPPUNIT_TEST_SUITE_END();

public:
  void setUp() { }

  void tearDown() { }

  void
  TestTooShort()
  {
    vector<double> cdf(1, 0);
    RandomSamplerFromCDF sampler(cdf);
  }

  void
  TestStartZero()
  {
    vector<double> cdf(2, 1);
    RandomSamplerFromCDF sampler(cdf);
  }

  void
  TestMonotonous()
  {
    vector<double> cdf(3);
    cdf[1] = 1;
    cdf[2] = 0.5;
    RandomSamplerFromCDF sampler(cdf);
  }

  void
  TestTabulatedFunction()
  {
    TabulatedFunction cdf;
    cdf.PushBack(-1, 0);
    cdf.PushBack(1, 1);
    cdf.PushBack(2, 1);
    cdf.PushBack(3, 2);
    const RandomSamplerFromCDF sampler(cdf);
    CLHEP::MTwistEngine rand;
    for (int i = 0; i < 1000; ++i) {
      const double x = sampler.shoot(rand);
      CPPUNIT_ASSERT(-1 <= x && x <= 3 && !(1 < x && x < 2));
    }
    const double eps = 1e-12;
    EQUAL(sampler.GetInverseCDF(-1), -1.);
    EQUAL(sampler.GetInverseCDF(0), -1.);
    CLOSE(sampler.GetInverseCDF(0.25), 0.);
    CLOSE(sampler.GetInverseCDF(0.5-eps), 1.);
    CLOSE(sampler.GetInverseCDF(0.5), 1.);
    CLOSE(sampler.GetInverseCDF(0.5+eps), 2.);
    CLOSE(sampler.GetInverseCDF(0.75), 2.5);
    EQUAL(sampler.GetInverseCDF(1), 3.);
    EQUAL(sampler.GetInverseCDF(2), 3.);
  }

  void
  TestVectors()
  {
    const double x[] = { 10,  10,  11, 11 };
    const double y[] = {  0, 0.1, 0.1,  1 };
    const unsigned int n = sizeof(x) / sizeof(x[0]);
    vector<double> xx(x, x + n);
    vector<double> yy(y, y + n);
    const RandomSamplerFromCDF sampler(xx, yy);
    CLHEP::MTwistEngine rand;
    for (int i = 0; i < 1000; ++i) {
      const double x = sampler.shoot(rand);
      CPPUNIT_ASSERT(x == 10 || x == 11);
    }
  }

  void
  TestVectorAndNormalization()
  {
    const double y[] = { 0, 1, 3, 7, 10 };
    const unsigned int n = sizeof(y) / sizeof(y[0]);
    const vector<double> yy(y, y + n);
    const RandomSamplerFromCDF sampler(yy);
    CLHEP::MTwistEngine rand;
    for (int i = 0; i < 1000; ++i) {
      const double x = sampler.shoot(rand);
      CPPUNIT_ASSERT(0 <= x && x <= 1);
    }
  }

  void
  TestMap()
  {
    const double x[] = { -5,  0,  8,  9, 9.5,  10 };
    const double y[] = {  0, 25, 25, 50,  50, 100 };
    map<double, double> xy;
    const unsigned int n = sizeof(x) / sizeof(x[0]);
    for (unsigned int i = 0; i < n; ++i)
      xy[x[i]] = y[i];
    const RandomSamplerFromCDF sampler(xy);
    CLHEP::MTwistEngine rand;
    for (int i = 0; i < 1000; ++i) {
      const double x = sampler.shoot(rand);
      CPPUNIT_ASSERT((-5 <= x && x <= 0) || (8 <= x && x <= 9) || (9.5 <= x && x <= 10));
    }
  }

  void
  TestFunction()
  {
    const RandomSamplerFromCDF sampler("x*x - 1", 1, 2, 10);
    CLOSE(sampler.GetInverseCDF(0), 1.);
    CLOSE_AT(sampler.GetInverseCDF(0.5), sqrt(3*0.5+1), 1e-2);
    CLOSE(sampler.GetInverseCDF(1), 2.);
    CLHEP::MTwistEngine rand;
    for (int i = 0; i < 1000; ++i) {
      const double x = sampler.shoot(rand);
      CPPUNIT_ASSERT((1 <= x && x <= 2));
    }
  }

  void
  TestFunction2()
  {
    const double min = -1.222;
    const double max = 0.63;
    const string func = "sin(x) - sin(" + lexical_cast<string>(min) + ")";
    const RandomSamplerFromCDF sampler(func, min, max, 1000);
    CLOSE(sampler.GetInverseCDF(0), min);
    CLOSE(sampler.GetInverseCDF(1), max);
    CLHEP::MTwistEngine rand;
    for (int i = 0; i < 1000; ++i) {
      const double x = sampler.shoot(rand);
      CPPUNIT_ASSERT((min <= x && x <= max));
    }
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(TestRandomSamplerFromCDF);