#ifndef _utl_RandomSamplerFromPDF_h_
#define _utl_RandomSamplerFromPDF_h_

#include <vector>
#include <map>
#include <string>

#include <utl/VRandomSampler.h>
#include <utl/ShadowPtr.h>


namespace utl {

  class TabulatedFunction;
  class RandomSamplerFromCDF;

  /**
    \class RandomSamplerFromPDF RandomSamplerFromPDF.h "utl/RandomSamplerFromPDF.h"

    This class is based on a re-write of CLHEP's randGeneral class with a couple of
    additions and changes.

    It can be constructed using a TabulatedFunction as the probability distribution
    function. It can also be constructed giving a string formula (as in ROOT).

    If you construct the distribution with an std::vector<double> as the PDF it will
    behave the same way RandGeneral does (changing double* by std::vector<double>, of
    course). In this case it will give numbers in [0,1) dividing the interval in
    subintervals of equal length.

    In all cases, there are different InterpolationTypes that specify whether the
    distribution will be discrete (only the values specified in the TabulatedFunction),
    a step function or linerly interpolated.

    Input does not have to be normalized.

    \author Darko Veberic
    \date 29 Aug 2009
    \version $Id: RandomSamplerFromPDF.h 14726 2009-09-19 15:46:30Z rodolfo $
    \ingroup math
  */

  class RandomSamplerFromPDF : public VRandomSampler {

  public:
    /** Interpolation types
      Depending on the interpolation type chosen at construction, the distribution of
      values will conform to:
      - a discrete pdf over the values in the x axis of the TabulatedFunction, with
        relative probabilities given by the y values of the TabulatedFunction,
      - a step pdf, constant in between consecutive x values and relative probability
        densities given by the y values (note that in this case the y value in the
        last point of the TabulatedFunction is discarded), or
      - a linearly interpolated, continuous function.

      The behaviour when creating with a vector is the same except the range of the
      function is [0, 1], equally subdivided.
    */
    enum InterpolationType {
      eLinear,
      eDiscrete,
      eStep
    };

    /// Construct using a PDF defined by an X and Y vectors
    RandomSamplerFromPDF(const std::vector<double>& x, const std::vector<double>& y,
                         const InterpolationType intType = eStep)
    { Init(x, y, intType); }

    /// Construct using a PDF defined by a Tabulated Function
    RandomSamplerFromPDF(const TabulatedFunction& pdf,
                         const InterpolationType intType = eStep);

    /// Construct using a PDF defined by an std::vector<double>, the old CLHEP behaviour
    RandomSamplerFromPDF(const std::vector<double>& pdf,
                         const InterpolationType intType = eStep)
    { Init(pdf, intType); }

    RandomSamplerFromPDF(const std::map<double, double>& pdf,
                         const InterpolationType intType = eDiscrete);

    /// Construct the PDF from a function defined by a string
    RandomSamplerFromPDF(const std::string& function,
                         const std::string& independentVariableName,
                         const double min, const double max, const int bins = 100)
    { SetFunction(function, independentVariableName, min, max, bins); }

    /// Construct the PDF from a function defined by a string
    RandomSamplerFromPDF(const std::string& function,
                         const double min, const double max, const int bins = 100)
    { SetFunction(function, "x", min, max, bins); }

    virtual ~RandomSamplerFromPDF();

    double GetInverseCDF(const double y) const;

  protected:
    virtual double MapRandom(const double rand) const
    { return GetInverseCDF(rand); }

  private:
    void SetFunction(const std::string& function,
                     const std::string& independentVariableName,
                     const double min, const double max, const int bins);

    void Init(const std::vector<double>& x, const std::vector<double>& y,
              const InterpolationType intType);

    void Init(const std::vector<double>& y,
              const InterpolationType intType);

    void InitDiscrete(const std::vector<double>& x, const std::vector<double>& y);

    void InitStep(const std::vector<double>& x, const std::vector<double>& y);

    void InitLinear(const std::vector<double>& x, const std::vector<double>& y);

    RandomSamplerFromCDF* fCDF;
    ShadowPtr<std::vector<double> > fPDF;
  };

}


#endif