RandomHelpers.h

// -*- C++ -*-
//
// RandomHelpers.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_RandomHelpers_H
#define HERWIG_RandomHelpers_H
 
#include "ThePEG/Config/ThePEG.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
namespace RandomHelpers {
 
inline double sign(double x) {
  return x < 0. ? -1. : 1.;
}
 
template<class Density>
struct Generator {
 
  double lower() const;
 
  double upper() const;
 
  double value(double x) const;
 
  double normalization() const;
 
  double operator()(double r) const;
 
};
 
struct Expression {};
 
template<>
struct Generator<Expression> {
 
  virtual ~Generator() {}
 
  virtual double lower() const = 0;
 
  virtual double upper() const = 0;
 
  virtual double value(double x) const = 0;
 
  virtual double normalization() const = 0;
 
  virtual double operator()(double r) const = 0;
 
};
 
template<class Density>
struct Container {};
 
template<class Density>
class Generator<Container<Density> >
  : public Generator<Expression> {
 
  Generator<Density> generator;
 
public:
 
  Generator(const Generator<Density>& gen)
    : generator(gen) {}
 
  virtual double lower() const { return generator.lower(); }
 
  virtual double upper() const { return generator.upper(); }
 
  virtual double value(double x) const { return generator.value(x); }
 
  virtual double normalization() const { return generator.normalization(); }
 
  virtual double operator()(double r) const { return generator(r); }
 
};
 
template<class Density>
pair<double,double> generate(const Generator<Density>& gen,
                             double r) {
  double x = gen(r);
 
  if ( gen.value(x) != 0. )
    return make_pair(x,gen.normalization()/gen.value(x));
  else
    return make_pair(x,0.);
}
 
template<class Density>
struct Remap {};
 
template<class Density>
class Generator<Remap<Density> > {
 
  Generator<Density> theGenerator;
 
  double theLower;
 
  double theUpper;
 
  Generator(const Generator<Density>& gen,
            double low, double up)
    : theGenerator(gen), theLower(low), theUpper(up) {
    if ( low >= up )
      throw std::logic_error("[Generator<Remap>] Invalid boundaries.");
  }
 
  const Generator<Density>& generator() const { return theGenerator; }
 
  double lower() const { return theLower; }
 
  double upper() const { return theUpper; }
 
  double value(double y) const {
    double xm = generator().lower();
    double xp = generator().upper();
    double ym = lower();
    double yp = upper();
    double x = ((xp-xm)/(yp-ym))*y+(yp*xm-ym*xp)/(yp-ym);
    return generator().value(x);
  }
 
  double normalization() const {
    double xm = generator().lower();
    double xp = generator().upper();
    double ym = lower();
    double yp = upper();
    return ((yp-ym)/(xp-xm))*generator().normalization();
  }
 
  double operator()(double r) const {
    double xm = generator().lower();
    double xp = generator().upper();
    double ym = lower();
    double yp = upper();
    double x = ((yp-ym)/(xp-xm))*generator()(r)+(xp*ym-xm*yp)/(xp-xm);
    return x;
  }
 
};
 
struct on {
 
  double lower;
 
  double upper;
 
  on(double a, double b)
    :lower(a), upper(b) {}
};
 
template<class Density>
Generator<Remap<Density> > operator*(const Generator<Density>& gen,
                                     const on& interval) {
  return Generator<Remap<Density> >(gen,interval.lower,interval.upper);
}
 
 
template<class Density>
struct Rescale {};
 
template<class Density>
class Generator<Rescale<Density> > {
 
  Generator<Density> theGenerator;
 
  double theScale;
 
public:
 
  Generator(const Generator<Density>& gen,
            double sc)
    : theGenerator(gen), theScale(sc) {
  }
 
  const Generator<Density>& generator() const { return theGenerator; }
 
  double scale() const { return theScale; }
 
  double lower() const { return generator().lower(); }
 
  double upper() const { return generator().upper(); }
 
  double value(double x) const {
    return scale()*generator().value(x);
  }
 
  double normalization() const {
    return scale()*generator().normalization();
  }
 
  double operator()(double r) const {
    return generator()(r);
  }
 
};
 
template<class Density>
Generator<Rescale<Density> > operator*(double a, const Generator<Density>& gen) {
  return Generator<Rescale<Density> >(gen,a);
}
 
template<class Density1,
         class Density2>
struct Sum {};
 
template<class Density1,
         class Density2>
class Generator<Sum<Density1,Density2> > {
 
  Generator<Density1> theFirstGenerator;
 
  Generator<Density2> theSecondGenerator;
 
  double theLower;
 
  double theUpper;
 
  double theFraction;
 
public:
 
  Generator(const Generator<Density1>& firstGen,
            const Generator<Density2>& secondGen)
    : theFirstGenerator(firstGen), theSecondGenerator(secondGen),
      theLower(min(firstGen.lower(),secondGen.lower())), 
      theUpper(max(firstGen.upper(),secondGen.upper())), 
      theFraction(1.) {
    theFraction = 
      firstGenerator().normalization() / normalization();
  }
 
  const Generator<Density1>& firstGenerator() const { return theFirstGenerator; }
 
  const Generator<Density2>& secondGenerator() const { return theSecondGenerator; }
 
  double lower() const { return theLower; }
 
  double upper() const { return theUpper; }
 
  double fraction() const { return theFraction; }
 
  double value(double x) const {
    double res = 0.;
    if ( firstGenerator().lower() <= x &&
         x <= firstGenerator().upper() )
      res += firstGenerator().value(x);
    if ( secondGenerator().lower() <= x &&
         x <= secondGenerator().upper() )
      res += secondGenerator().value(x);
    return res;
  }
 
  double normalization() const {
    return 
      firstGenerator().normalization() + secondGenerator().normalization();
  }
 
  double operator()(double r) const {
    return
      r < fraction() ? 
      firstGenerator()(r/fraction()) : 
      secondGenerator()((r-fraction())/(1.-fraction()));
  }
 
};
 
template<class Density1,
         class Density2>
Generator<Sum<Density1,Density2> > operator+(const Generator<Density1>& first,
                                             const Generator<Density2>& second) {
  return Generator<Sum<Density1,Density2> >(first,second);
}
 
template<class Density>
struct matcher {
 
  Generator<Density> generator;
 
  matcher(const Generator<Density>& gen)
    : generator(gen) {}
 
};
 
template<class Density>
matcher<Density> match(const Generator<Density>& gen) {
  return matcher<Density>(gen);
}
 
template<class Density1,
         class Density2>
Generator<Sum<Density1,Rescale<Density2> > > operator+(const Generator<Density1>& first,
                                                       const matcher<Density2>& second) {
  double matching = 
    first.value(first.upper())/
    second.generator.value(second.generator.lower());
  return Generator<Sum<Density1,Rescale<Density2> > >(first,matching*second.generator);
}
 
template<class Density1,
         class Density2>
struct Piecewise {};
 
struct ToBeDefined {};
 
template<class Density1,
         class Density2>
class Generator<Piecewise<Density1,Density2> > {
 
  Generator<Density1> theFirstGenerator;
 
  Generator<Density2> theSecondGenerator;
 
  double theLower;
 
  double theIntermediate;
 
  double theUpper;
 
  double theFraction;
 
public:
 
  Generator(const Generator<Density1>& firstGen,
            const Generator<Density2>& secondGen)
    : theFirstGenerator(firstGen), theSecondGenerator(secondGen),
      theLower(firstGen.lower()), theIntermediate(firstGen.upper()), theUpper(secondGen.upper()), 
      theFraction(1.) {
    if ( firstGenerator().upper() != secondGenerator().lower() )
      throw std::logic_error("[Generator<Piecewise>] Invalid boundaries.");
    theFraction = 
      firstGenerator().normalization() / normalization();
  }
 
  const Generator<Density1>& firstGenerator() const { return theFirstGenerator; }
 
  const Generator<Density2>& secondGenerator() const { return theSecondGenerator; }
 
  double lower() const { return theLower; }
 
  double intermediate() const { return theIntermediate; }
 
  double upper() const { return theUpper; }
 
  double fraction() const { return theFraction; }
 
  double value(double x) const {
    return
      x < intermediate() ? 
      firstGenerator().value(x) : 
      secondGenerator().value(x);
  }
 
  double normalization() const {
    return 
      firstGenerator().normalization() + secondGenerator().normalization();
  }
 
  double operator()(double r) const {
    return
      r < fraction() ? 
      firstGenerator()(r/fraction()) : 
      secondGenerator()((r-fraction())/(1.-fraction()));
  }
 
  template<class Density3>
  Generator<Piecewise<Piecewise<Density1,Density2>,Density3> >
  operator,(const Generator<Density3>& thirdGenerator) {
    return 
      Generator<Piecewise<Piecewise<Density1,Density2>,Density3> >
      (*this,thirdGenerator);
  }
 
  template<class Density3>
  Generator<Piecewise<Piecewise<Density1,Density2>,Rescale<Density3> > >
  operator,(const matcher<Density3>& thirdGenerator) {
    double matching = 
      value(upper())/thirdGenerator.generator.value(upper());
    return
      Generator<Piecewise<Piecewise<Density1,Density2>,Rescale<Density3> > >
      (*this,matching*thirdGenerator.generator);
  }
 
};
 
template<class Density>
struct Generator<Piecewise<Density,ToBeDefined> > {
 
  Generator<Density> generator;
 
  Generator(const Generator<Density>& gen)
    : generator(gen) {}
 
  template<class Density2>
  Generator<Piecewise<Density,Density2> >
  operator,(const Generator<Density2>& secondGen) {
    return 
      Generator<Piecewise<Density,Density2> >
      (generator,secondGen);
  }
 
  template<class Density2>
  Generator<Piecewise<Density,Rescale<Density2> > >
  operator,(const matcher<Density2>& secondGen) {
    double matching = 
      generator.value(generator.upper())/secondGen.generator.value(generator.upper());
    return
      Generator<Piecewise<Density,Rescale<Density2> > >
      (generator,matching*secondGen.generator);
  }
 
};
 
template<>
struct Generator<Piecewise<ToBeDefined,ToBeDefined> > {
 
  template<class Density>
  Generator<Piecewise<Density,ToBeDefined> >
  operator,(const Generator<Density>& gen) {
    return 
      Generator<Piecewise<Density,ToBeDefined> >(gen);
  }
 
};
 
inline Generator<Piecewise<ToBeDefined,ToBeDefined> >
piecewise() {
  return Generator<Piecewise<ToBeDefined,ToBeDefined> >();
}
 
 
 
struct Flat {};
 
template<>
class Generator<Flat> {
 
  double theLower;
 
  double theUpper;
 
public:
 
  Generator(double low, double up)
    : theLower(low), theUpper(up) {}
 
  double lower() const { return theLower; }
 
  double upper() const { return theUpper; }
 
  double value(double x) const { 
    return x>=lower() && x<=upper() ? 1. : 0.;
  }
 
  double normalization() const { return upper()-lower(); }
 
  double operator()(double r) const { 
    return lower() + r*(upper()-lower());
  }
 
};
 
inline Generator<Flat> flat(double low, double up) {
  return Generator<Flat>(low,up);
}
 
struct Zero {};
 
template<>
class Generator<Zero> {
 
  double theLower;
 
  double theUpper;
 
public:
 
  Generator(double low, double up)
    : theLower(low), theUpper(up) {}
 
  double lower() const { return theLower; }
 
  double upper() const { return theUpper; }
 
  double value(double x) const { 
    return x>=lower() && x<=upper() ? Constants::epsilon : 0.;
  }
 
  double normalization() const { return 0.; }
 
  double operator()(double r) const { 
    return lower() + r*(upper()-lower());
  }
 
};
 
inline Generator<Zero> zero(double low, double up) {
  return Generator<Zero>(low,up);
}
 
struct Inverse {};
 
template<>
class Generator<Inverse> {
 
  double thePole;
 
  double theLower;
 
  double theUpper;
 
  double theScale;
 
  double theOffset;
 
public:
 
  Generator(double z,
            double l, double u)
    : thePole(z),
      theLower(l), theUpper(u),
      theScale(z < l ? log((u-z)/(l-z)) : log((z-l)/(z-u))),
      theOffset(z < l ? log(l-z) : log(z-u)) {
    if ( z >= l && z <= u )
      throw std::logic_error("[Generator<Inverse>] Pole inside sampling interval.");
  }
 
  double lower() const { return theLower; }
 
  double upper() const { return theUpper; }
 
  double pole() const { return thePole; }
 
  double scale() const { return theScale; }
 
  double offset() const { return theOffset; }
 
  double value(double x) const {
    return x>=lower() && x<=upper() ? 1/abs(x-pole()) : 0.;
  }
 
  double normalization() const { 
    return scale();
  }
 
  double operator()(double r) const { 
    return pole() + sign(upper()-pole())*exp(scale()*r+offset());
  }
 
};
 
inline Generator<Inverse> inverse(double z,
                                  double lower, double upper) {
  return Generator<Inverse>(z,lower,upper);
}
 
struct Power {};
 
template<>
class Generator<Power> {
 
  double thePole;
 
  double thePower;
 
  double theLower;
 
  double theUpper;
 
  double theScale;
 
  double theOffset;
 
public:
 
  Generator(double z, double p,
            double l, double u)
    : thePole(z), thePower(p),
      theLower(l), theUpper(u),
      theScale(z<=l ? (pow(u-z,1.+p)-pow(l-z,1.+p))/(1.+p) : (pow(z-l,1.+p)-pow(z-u,1.+p))/(1.+p)),
      theOffset(z<=l ? pow(l-z,1.+p)/(1.+p) : pow(z-u,1.+p)/(1.+p)) {
    if ( p == -1. )
      throw std::logic_error("[Generator<Power>] Unit inverse. Consider using inverse().");
    if ( z >= l && z <= u && p < 0. )
      throw std::logic_error("[Generator<Power>] Pole inside sampling interval.");
    if ( z >= l && z <= u && p > 0. )
      throw std::logic_error("[Generator<Power>] Zero inside sampling interval.");
  }
 
  double lower() const { return theLower; }
 
  double upper() const { return theUpper; }
 
  double pole() const { return thePole; }
 
  double power() const { return thePower; }
 
  double scale() const { return theScale; }
 
  double offset() const { return theOffset; }
 
  double value(double x) const {
    return x>=lower() && x<=upper() ? pow(abs(x-pole()),power()) : 0.;
  }
 
  double normalization() const { 
    return scale();
  }
 
  double operator()(double r) const { 
    return pole() + sign(upper()-pole())*pow((1.+power())*(scale()*r+offset()),1./(1.+power()));
  }
 
};
 
inline Generator<Power> power(double z, double p,
                              double lower, double upper) {
  return Generator<Power>(z,p,lower,upper);
}
 
struct BreitWigner {};
 
template<>
class Generator<BreitWigner> {
 
  double thePole;
 
  double theWidth;
 
  double theLower;
 
  double theUpper;
 
  double theScale;
 
  double theOffset;
 
  double theSqrtWidth;
 
public:
 
  Generator(double z, double w,
            double l, double u)
    : thePole(z), theWidth(w),
      theLower(l), theUpper(u),
      theScale((atan((u-z)/sqrt(abs(w*z)))-atan((l-z)/sqrt(abs(w*z))))/sqrt(abs(w*z))),
      theOffset(atan((l-z)/sqrt(abs(w*z)))/sqrt(abs(w*z))),
      theSqrtWidth(sqrt(abs(w*z))) {
    if ( w == 0. )
      throw std::logic_error("[Generator<BreitWigner>] Zero width. Consider using power().");
  }
 
  double lower() const { return theLower; }
 
  double upper() const { return theUpper; }
 
  double pole() const { return thePole; }
 
  double width() const { return theWidth; }
 
  double scale() const { return theScale; }
 
  double offset() const { return theOffset; }
 
  double sqrtWidth() const { return theSqrtWidth; }
 
  double value(double x) const {
    return 
      x>=lower() && x<=upper() ?
      1./(sqr(x-pole())+abs(width()*pole())) : 0.;
  }
 
  double normalization() const { 
    return scale();
  }
 
  double operator()(double r) const { 
    double res = pole() + sqrtWidth()*tan(sqrtWidth()*(scale()*r+offset()));
    if ( res <= lower() ) return lower()*(1+std::numeric_limits<double>::epsilon());
    else if ( res >= upper() ) return upper()*(1-std::numeric_limits<double>::epsilon());
    else return res;
  }
 
};
 
inline Generator<BreitWigner> breitWigner(double z, double w,
                                          double lower, double upper) {
  return Generator<BreitWigner>(z,w,lower,upper);
}
 
}
 
}
 
#endif // HERWIG_RandomHelpers_H