TabulatedPDF.cc

Go to the documentation of this file.

#include "TabulatedPDF.h"

#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>

#include <sstream>
#include <cmath>
#include <cfloat>
#include <algorithm>

using namespace std;
using namespace utl;
using namespace tls;

typedef TabulatedPDF::Array Array;

TabulatedPDF::TabulatedPDF()
{
}

TabulatedPDF::TabulatedPDF(const TabulatedPDF& src)
{ *this = src; }

TabulatedPDF::TabulatedPDF(
  double xmin, double xstep, const std::vector< double > &yValues)
{
  FillTable(xmin, xstep, yValues);
}

TabulatedPDF&
TabulatedPDF::operator=(const TabulatedPDF& src)
{
  fXMin  = src.fXMin;
  fXStep = src.fXStep;
  fY     = src.fY;
  fYSum  = src.fYSum;
  return *this;
}

void
TabulatedPDF::FillTable(
  double xmin, double xstep, const Array& yValues)
{
  const unsigned int size = yValues.size();

  Clear();

  fXMin  = xmin;
  fXStep = xstep;
  for (unsigned int i=0;i<size;++i) PushBack(yValues[i]);
  Normalise();
}

void
TabulatedPDF::PushBack(double y)
{
  const unsigned int size = fY.size();

  if (y < 0.0) {
    ostringstream msg;
    msg << "Value " << y << " at x = " << fXMin+(size-1)*fXStep << " is not a valid probability!";
    ERROR(msg);
    y = 0.0;
  }

  if (size==0)
    fYSum.push_back(0.5*y*fXStep);
  else {
    const double ysum = fYSum.back() + 0.5*(y+fY.back())*fXStep;
    fYSum.push_back(ysum);
  }

  fY.push_back(y);
}

void
TabulatedPDF::Clear()
{
  fXMin = 0.0;
  fXStep = 0.0;
  Array().swap(fY);
  Array().swap(fYSum);
}

void
TabulatedPDF::Normalise()
{
  if (fYSum.back()==1.0) return;

  const double factor = 1./fYSum.back();

  for (unsigned int i=0;i<fY.size();++i) {
    fY[i]    *= factor;
    fYSum[i] *= factor;
  }
}

double
TabulatedPDF::Density(double x)
  const
{
  const unsigned int size = fY.size();

  if (!size) throw NonExistentComponentException("No points in table.");

  if (size == 1)
    return fY[0];

  if (x < fXMin-fXStep/2)
    return 0.0;
  if (x < fXMin)
    return fY[0];
  else if (x > fXMin + (size-1)*fXStep)
    return 0.0;

  unsigned int index = static_cast<unsigned int>((x-fXMin)/fXStep);
  const double x1 = fXMin+index*fXStep;
  const double& dx = fXStep;
  const double y1 = fY[index];
  if (x1 == x) return y1;
  ++index;
  const double y2 = fY[index];

  double prob = y1 + (y2 - y1)*(x - x1)/dx;
  return prob;
}

double
TabulatedPDF::Integral(double xmin, double xmax)
  const
{
  const unsigned int size = fY.size();

  if (!size) throw NonExistentComponentException("No points in table.");

  if (size == 1)
    return fY[0];

  // primitive out of bounds handling
  if (xmin > fXMin + (size-1)*fXStep)
    return 0.0;
  if (xmax < fXMin)
    return 0.0;
  if (xmin < fXMin)
    xmin = fXMin;
  if (xmax > fXMin + (size-1)*fXStep)
    xmax = fXMin + (size-2)*fXStep;

  const unsigned int iBegin = static_cast<unsigned int>((xmin-fXMin)/fXStep);
  const unsigned int iEnd   = static_cast<unsigned int>((xmax-fXMin)/fXStep)-1;
  const double& dx = fXStep;
  const double x1l = fXMin + iBegin*fXStep;
  const double y1l = fYSum[iBegin];
  const double y1u = fYSum[iBegin+1];

  const double x2l = fXMin + iEnd*fXStep;
  const double y2l = fYSum[iEnd];
  const double y2u = fYSum[iEnd+1];

  const double p1 = y1l + (y1u - y1l)*(xmin - x1l)/dx;
  const double p2 = y2l + (y2u - y2l)*(xmax - x2l)/dx;
  return p2-p1;
}

double
TabulatedPDF::Moment(unsigned short order, double xmin, double xmax)
  const
{
  const unsigned int size = fY.size();

  if (!size) throw NonExistentComponentException("No points in table.");

  if (size==1)
    return fY[0]*pow(fXMin,order);

  // primitive out of bounds handling
  if (xmin < fXMin)
    xmin = fXMin;
  else if (xmin > fXMin + (size-1)*fXStep)
    return 0.0;
  if (xmax < fXMin)
    return 0.0;
  else if (xmax > fXMin + (size-1)*fXStep)
    xmax = fXMin + (size-1)*fXStep;

  const unsigned int iBegin = static_cast<unsigned int>((xmin-fXMin)/fXStep);
  const unsigned int iEnd   = static_cast<unsigned int>((xmax-fXMin)/fXStep);

  // simpson's rule for integration
  // formula R = h ( I0 + 4 I1 + 2 I2 + 4 I3 + ... + 2 In-2 + 4 In-1 + In) / 3
  double mom = 0;
  double sum = 0;

  // first term
  const double& h = fXStep;
  mom = h*fY[iBegin]*pow(fXMin+iBegin*fXStep,order);
  sum = h*fY[iBegin];

  // intermediate terms
  bool iloop = true;
  for (unsigned int i=iBegin+1;i<iEnd;++i)
  {
    const int weight = iloop ? 4:2 ;
    iloop=!iloop;
    mom += weight*h*fY[i]*pow(fXMin+i*fXStep,order);
    sum += weight*h*fY[i];
  }

  // last term
  mom += h*fY[iEnd-1]*pow(fXMin+iEnd*fXStep,order);
  sum += h*fY[iEnd-1];

  return mom/sum;
}

double
TabulatedPDF::Variance(double xmin, double xmax)
  const
{
  const unsigned int size = fY.size();

  if (!size) throw NonExistentComponentException("No points in table.");

  if (size==1)
    return 0;

  // primitive out of bounds handling
  if (xmin < fXMin)
    xmin = fXMin;
  else if (xmin > fXMin + (size-1)*fXStep)
    return 0.0;
  if (xmax < fXMin)
    return 0.0;
  else if (xmax > fXMin + (size-1)*fXStep)
    xmax = fXMin + (size-1)*fXStep;

  const unsigned int iBegin = static_cast<unsigned int>((xmin-fXMin)/fXStep);
  const unsigned int iEnd   = static_cast<unsigned int>((xmax-fXMin)/fXStep);

  // simpson's rule for integration
  // formula R = h ( I0 + 4 I1 + 2 I2 + 4 I3 + ... + 2 In-2 + 4 In-1 + In) / 3
  double mom1 = 0;
  double mom2 = 0;
  double sum = 0;

  // first term
  const double& h = fXStep;
  mom1 = h*fY[iBegin]*(fXMin+iBegin*fXStep);
  mom2 = h*fY[iBegin]*(fXMin+iBegin*fXStep)*(fXMin+iBegin*fXStep);
  sum = h*fY[iBegin];

  // intermediate terms
  bool iloop = true;
  for (unsigned int i=iBegin+1;i<iEnd;++i)
  {
    const int weight = iloop ? 4:2 ;
    iloop=!iloop;
    mom1 += weight*h*fY[i]*(fXMin+i*fXStep);
    mom2 += weight*h*fY[i]*(fXMin+i*fXStep)*(fXMin+i*fXStep);
    sum += weight*h*fY[i];
  }

  // last term
  mom1 += h*fY[iEnd]*(fXMin+iEnd*fXStep);
  mom2 += h*fY[iEnd]*(fXMin+iEnd*fXStep)*(fXMin+iEnd*fXStep);
  sum += h*fY[iEnd];

  return mom2/sum-pow(mom1/sum,2);
}

double
TabulatedPDF::Quantile(const double psum)
  const
{
  const unsigned int size = fY.size();

  if (!size) throw NonExistentComponentException("No points in table.");

  int index = 0;
  for (unsigned int i=0;i<size-1;++i)
    if (fYSum[i] < psum && fYSum[i+1] > psum)
      { index = i; break; }
  const double x0 = fXMin+index*fXStep;
  const double& dx = fXStep;
  const double& y0 = fYSum[index];
  const double& y1 = fYSum[index+1];
  const double quantile = (psum - y0)/(y1-y0)*dx+x0;
  return quantile;
}