testRandomSamplerFromPDF.cc

Go to the documentation of this file.

#include <tst/Verify.h>
#include <cppunit/extensions/HelperMacros.h>
#include <utl/RandomSamplerFromPDF.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 TestRandomSamplerFromPDF : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(TestRandomSamplerFromPDF);
  CPPUNIT_TEST_EXCEPTION(TestLinearTooShort, InvalidConfigurationException);
  CPPUNIT_TEST_EXCEPTION(TestDiscreteTooShort, InvalidConfigurationException);
  CPPUNIT_TEST_EXCEPTION(TestStepTooShort, InvalidConfigurationException);
  CPPUNIT_TEST_EXCEPTION(TestNegative, 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
  TestLinearTooShort()
  {
    vector<double> pdf(1, 0);
    RandomSamplerFromPDF sampler(pdf, RandomSamplerFromPDF::eLinear);
  }

  void
  TestDiscreteTooShort()
  {
    vector<double> pdf;
    RandomSamplerFromPDF sampler(pdf, RandomSamplerFromPDF::eDiscrete);
  }

  void
  TestStepTooShort()
  {
    vector<double> pdf;
    RandomSamplerFromPDF sampler(pdf, RandomSamplerFromPDF::eStep);
  }

  void
  TestNegative()
  {
    vector<double> pdf(3);
    pdf[1] = 1;
    pdf[2] = -0.5;
    RandomSamplerFromPDF sampler(pdf, RandomSamplerFromPDF::eDiscrete);
  }

  void
  TestTabulatedFunction()
  {
    TabulatedFunction pdf;
    pdf.PushBack(-1, 0);
    pdf.PushBack(1, 1);
    pdf.PushBack(2, 0);
    pdf.PushBack(4, 0);
    pdf.PushBack(5, 0.5);
    pdf.PushBack(6, 0);
    pdf.PushBack(9, 0);
    pdf.PushBack(10, 3);
    CLHEP::MTwistEngine rand;
    {
      const RandomSamplerFromPDF sLin(pdf, RandomSamplerFromPDF::eLinear);
      for (int i = 0; i < 1000; ++i) {
        const double x = sLin.shoot(rand);
        CPPUNIT_ASSERT((-1 <= x && x <= 2) || (4 <= x && x <= 6) || (9 <= x && x <= 10));
      }
      const double x[] = {
        -1, -0.5,  0, 0.5, 1,   1.5, 2, 4.5, 5, 5.5,
         6,  9.5, 10, 4.75, 9.2, 4.8, 9.9
      };
      const double y[] = {
        0, 0.0625, 0.25, 0.5625, 1,   1.375, 1.5, 1.5625, 1.75, 1.9375,
        2, 2.375,  3.5, 1.64063, 2.06, 1.66, 3.215
      };
      const unsigned int n = sizeof(x) / sizeof(x[0]);
      for (unsigned int i = 0; i < n; ++i)
        CLOSE_AT(sLin.GetInverseCDF(y[i]/3.5), x[i], 1e-5);
    }
    {
      const RandomSamplerFromPDF sDis(pdf, RandomSamplerFromPDF::eDiscrete);
      for (int i = 0; i < 1000; ++i) {
        const double x = sDis.shoot(rand);
        CPPUNIT_ASSERT(x == 1 || x == 5 || x == 10);
      }
      const double x[] = { -1,   1,   1,   1, 1,   5,  10 };
      const double y[] = {  0, 0.1, 0.5, 0.9, 1, 1.1, 4.5 };
      const unsigned int n = sizeof(x) / sizeof(x[0]);
      for (unsigned int i = 0; i < n; ++i)
        EQUAL(sDis.GetInverseCDF(y[i]/4.5), x[i]);
    }
    {
      const RandomSamplerFromPDF sSt(pdf, RandomSamplerFromPDF::eStep);
      for (int i = 0; i < 1000; ++i) {
        const double x = sSt.shoot(rand);
        CPPUNIT_ASSERT((1 <= x && x <= 2) || (5 <= x && x <= 6));
      }
      const double x[] = { 1.12, 1.93, 5.58 };
      const double y[] = { 0.08, 0.62, 0.86 };
      const unsigned int n = sizeof(x) / sizeof(x[0]);
      for (unsigned int i = 0; i < n; ++i)
        CLOSE(sSt.GetInverseCDF(y[i]), x[i]);
    }
  }

  /*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(TestRandomSamplerFromPDF);