RandomSamplerFromCDF.cc

Go to the documentation of this file.

#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>
#include <sstream>

using namespace std;
using namespace utl;


RandomSamplerFromCDF::RandomSamplerFromCDF(const TabulatedFunction& cdf) :
  fX(cdf.XBegin(), cdf.XEnd()),
  fY(cdf.YBegin(), cdf.YEnd())
{
  CheckAndNormalize();
}


RandomSamplerFromCDF::RandomSamplerFromCDF(const vector<double>& cdf) :
  fY(cdf.begin(), cdf.end())
{
  const int size = cdf.size();
  fX.resize(size);
  const double step = 1. / (size - 1);
  for (int i = 0; i < size; ++i)
    fX[i] = i * step;
  CheckAndNormalize();
}


RandomSamplerFromCDF::RandomSamplerFromCDF(const map<double, double>& cdf)
{
  for (map<double, double>::const_iterator it = cdf.begin();
       it != cdf.end(); ++it) {
    fX.push_back(it->first);
    fY.push_back(it->second);
  }
  CheckAndNormalize();
}


double
RandomSamplerFromCDF::GetInverseCDF(const double y)
  const
{
  if (y <= 0)
    return fX.front();
  if (y >= 1)
    return fX.back();

  const int i2 = GetIndex(y);
  const int i1 = i2 - 1;
  const double x1 = fX[i1];
  const double y1 = fY[i1];
  return x1 + (fX[i2] - x1)*(y - y1)/(fY[i2] - y1);
}


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


void
RandomSamplerFromCDF::CheckAndNormalize()
{
  const unsigned int n = fX.size();
  if (n < 2 || n != fY.size())
    throw InvalidConfigurationException("Not enough points.");
  if (fY.front() != 0) {
    ostringstream err;
    err << "CDF starts with " << fY.front() << " instead of 0.";
    throw InvalidConfigurationException(err.str());
  }
  vector<double>::iterator it = fY.begin();
  double y1 = *it;
  for (++it; it != fY.end(); ++it) {
    const double y2 = *it;
    if (y1 > y2) {
      cerr << "CDF is not monotonously increasing: " << y1 << " > " << y2 << endl;
      Dump();
      throw InvalidConfigurationException("CDF is not monotonously increasing.");
    }
    y1 = y2;
  }
  if (y1 != 1) {
    const double fac = 1 / y1;
    for (it = fY.begin(); it != fY.end(); ++it)
      *it *= fac;
  }
}


void
RandomSamplerFromCDF::Dump()
  const
{
  cerr << "CDF (x,y) dump:" << endl;
  const size_t n = fX.size();
  for (size_t i = 0; i < n; ++i)
    cerr << fX.at(i) << ' ' << fY.at(i) << endl;
}