RandomSamplerFromPDF.cc

Go to the documentation of this file.

#include <utl/RandomSamplerFromPDF.h>
#include <utl/RandomSamplerFromCDF.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/Function.h>
#include <utl/TabulatedFunction.h>

#include <cmath>

using namespace std;
using namespace utl;


RandomSamplerFromPDF::RandomSamplerFromPDF(const TabulatedFunction& pdf,
                                           const InterpolationType intType)
{
  const vector<double> x(pdf.XBegin(), pdf.XEnd());
  const vector<double> y(pdf.YBegin(), pdf.YEnd());
  Init(x, y, intType);
}


RandomSamplerFromPDF::RandomSamplerFromPDF(const map<double, double>& pdf,
                                           const InterpolationType intType)
{
  vector<double> x;
  vector<double> y;
  for (map<double, double>::const_iterator it = pdf.begin();
       it != pdf.end(); ++it) {
    x.push_back(it->first);
    y.push_back(it->second);
  }
  Init(x, y, intType);
}


RandomSamplerFromPDF::~RandomSamplerFromPDF()
{
  delete fCDF;
}


double
RandomSamplerFromPDF::GetInverseCDF(const double y)
  const
{
  if (!fPDF)
    return fCDF->GetInverseCDF(y);

  if (y <= 0)
    return fCDF->GetMinX();
  if (y >= 1)
    return fCDF->GetMaxX();

  const int i2 = fCDF->GetIndex(y);
  const int i1 = i2 - 1;
  const vector<double>& cdfx = fCDF->GetCDFX();
  const vector<double>& cdfy = fCDF->GetCDFY();
  const vector<double>& pdf = *fPDF;
  const double x1 = cdfx[i1];
  const double y1 = pdf[i1];
  const double w = y - cdfy[i1];
  const double k = (pdf[i2] - y1) / (cdfx[i2] - x1);
  const double r = sqrt(y1*y1 + 2*k*w);
  return x1 + 2*w / (r + y1);
}


void
RandomSamplerFromPDF::SetFunction(const string& function,
                                  const string& independentVariableName,
                                  const double min, const double max, const int bins)
{
  vector<double> x;
  x.resize(bins);
  vector<double> y;
  y.resize(bins);
  utl::Function func(function, independentVariableName);
  const double step = (max - min) / (bins - 1);
  for (int i = 0; i < bins; ++i) {
    const double xx = min + i * step;
    x[i] = xx;
    y[i] = func(xx);
  }
  InitLinear(x, y);
}


void
RandomSamplerFromPDF::Init(const vector<double>& x, const vector<double>& y,
                           const InterpolationType intType)
{
  switch (intType) {
  case eLinear:
    InitLinear(x, y);
    break;
  case eDiscrete:
    InitDiscrete(x, y);
    break;
  case eStep:
    InitStep(x, y);
    break;
  }
}


inline
void
UniformFill(vector<double>& x, const unsigned int n, const unsigned int m)
{
  const double step = 1. / m;
  for (unsigned int i = 0; i < n; ++i)
    x[i] = step * i;
}


void
RandomSamplerFromPDF::Init(const vector<double>& y,
                           const InterpolationType intType)
{
  const unsigned int n = y.size();
  vector<double> x;
  switch (intType) {
  case eLinear:
    x.resize(n);
    UniformFill(x, n, n-1);
    Init(x, y, intType);
    break;
  case eDiscrete:
    x.resize(n);
    UniformFill(x, n, n);
    Init(x, y, intType);
    break;
  case eStep:
    {
      x.resize(n+1);
      UniformFill(x, n+1, n);
      vector<double> yy(y);
      yy.push_back(0);
      Init(x, yy, intType);
      break;
    }
  }
}


inline
void
Normalize(vector<double>& v, const double fact)
{
  for (vector<double>::iterator it = v.begin();
       it != v.end(); ++it)
    *it *= fact;
}


void
RandomSamplerFromPDF::InitDiscrete(const vector<double>& x, const vector<double>& y)
{
  const unsigned int n = x.size();
  if (n < 1 || n != y.size())
    throw InvalidConfigurationException("Not enough points in PDF");
  const unsigned int n2 = 2*n;
  vector<double> xv;
  xv.resize(n2);
  vector<double> yv;
  yv.resize(n2);
  double sum = 0;
  for (unsigned int i = 0; i < n; ++i) {
    unsigned int i2 = 2*i;
    const double xx = x[i];
    xv[i2] = xx;
    yv[i2] = sum;
    ++i2;
    sum += y[i];
    xv[i2] = xx;
    yv[i2] = sum;
  }
  if (sum != 1)
    Normalize(yv, 1/sum);
  fCDF = new RandomSamplerFromCDF(xv, yv);
}


void
RandomSamplerFromPDF::InitStep(const vector<double>& x, const vector<double>& y)
{
  const unsigned int n = x.size();
  if (n < 2 || n != y.size())
    throw InvalidConfigurationException("Not enough points in PDF");
  vector<double> xv;
  xv.resize(n);
  vector<double> yv;
  yv.resize(n);
  xv[0] = x.front();
  yv[0] = 0;
  double sum = 0;
  for (unsigned int i = 1; i < n; ++i) {
    xv[i] = x[i];
    sum += y[i-1];
    yv[i] = sum;
  }
  if (sum != 1)
    Normalize(yv, 1/sum);
  fCDF = new RandomSamplerFromCDF(xv, yv);
}


void
RandomSamplerFromPDF::InitLinear(const vector<double>& x, const vector<double>& y)
{
  const unsigned int n = x.size();
  if (n < 2 || n != y.size())
    throw InvalidConfigurationException("Not enough points in PDF");
  vector<double> xv;
  xv.resize(n);
  vector<double> yv;
  yv.resize(n);
  fPDF = new vector<double>;
  fPDF->resize(n);
  vector<double>& pdf2 = *fPDF;
  double sum = 0;
  double x1 = x.front();
  xv[0] = x1;
  yv[0] = 0;
  double y1 = y.front();
  pdf2[0] = y1;
  for (unsigned int i = 1; i < n; ++i) {
    const double x2 = x[i];
    xv[i] = x2;
    const double y2 = y[i];
    pdf2[i] = y2;
    sum += 0.5 * (x2 - x1) * (y1 + y2);
    yv[i] = sum;
    x1 = x2;
    y1 = y2;
  }
  if (sum != 1) {
    if (sum == 0)
      throw InvalidConfigurationException("Sum of samples in PDF is zero, "
                                          "impossible to normalize");
    const double fact = 1 / sum;
    Normalize(yv, fact);
    Normalize(pdf2, fact);
  }
  fCDF = new RandomSamplerFromCDF(xv, yv);
}