<commit_before> /*****************************************************************************
  * Project: RooFit                                                           *
  *                                                                           *
  * Simple Gamma distribution
  * authors: <NAME>, <NAME>
  *                                                                           *
  *****************************************************************************/

/** \class RooGamma
    \ingroup Roofit

Implementation of the Gamma PDF for RooFit/RooStats.
\f[
f(x) = \frac{(x-\mu)^{\gamma-1} \cdot \exp^{(-(x-mu) / \beta}}{\Gamma(\gamma) \cdot \beta^{\gamma}}
\f]
defined for \f$ x \geq 0 \f$ if \f$ \mu = 0 \f$

Notes from K<NAME>:

Wikipedia and several sources refer to the Gamma distribution as

\f[
G(\mu,\alpha,\beta) = \beta^\alpha \mu^{(\alpha-1)} \frac{e^{(-\beta \mu)}}{\Gamma(\alpha)}
\f]

Below is the correspondence:

| Wikipedia       | This Implementation      |
|-----------------|--------------------------|
| \f$ \alpha \f$  | \f$ \gamma \f$           |
| \f$ \beta \f$   | \f$ \frac{1}{\beta} \f$  |
| \f$ \mu \f$     | x                        |
| 0               | \f$ \mu \f$              |


Note, that for a model Pois(N|mu), a uniform prior on mu, and a measurement N
the posterior is in the Wikipedia parameterization Gamma(mu, alpha=N+1, beta=1)
thus for this implementation it is:

`RooGamma(_x=mu,_gamma=N+1,_beta=1,_mu=0)`

Also note, that in this case it is equivalent to
RooPoison(N,mu) and treating the function as a PDF in mu.
**/

#include "RooGamma.h"

#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooRandom.h"
#include "RooMath.h"
#include "RooHelpers.h"
#include "BatchHelpers.h"
#include "RooVDTHeaders.h"

#include "TMath.h"
#include <Math/SpecFuncMathCore.h>
#include <Math/PdfFuncMathCore.h>
#include <Math/ProbFuncMathCore.h>

#include <iostream>
#include <cmath>
using namespace std;

ClassImp(RooGamma);

////////////////////////////////////////////////////////////////////////////////

RooGamma::RooGamma(const char *name, const char *title,
          RooAbsReal& _x, RooAbsReal& _gamma,
          RooAbsReal& _beta, RooAbsReal& _mu) :
  RooAbsPdf(name,title),
  x("x","Observable",this,_x),
  gamma("gamma","Mean",this,_gamma),
  beta("beta","Width",this,_beta),
  mu("mu","Para",this,_mu)
{
  RooHelpers::checkRangeOfParameters(this, {&_gamma, &_beta}, 0.);
}

////////////////////////////////////////////////////////////////////////////////

RooGamma::RooGamma(const RooGamma& other, const char* name) :
  RooAbsPdf(other,name), x("x",this,other.x), gamma("mean",this,other.gamma),
  beta("beta",this,other.beta), mu("mu",this,other.mu)
{
}

////////////////////////////////////////////////////////////////////////////////

Double_t RooGamma::evaluate() const
{
  return TMath::GammaDist(x, gamma, mu, beta) ;
}

////////////////////////////////////////////////////////////////////////////////

namespace {
//Author: <NAME>, CERN 22 August 2019

template<class Tx, class Tgamma, class Tbeta, class Tmu>
void compute(	size_t batchSize,
              double * __restrict output,
              Tx X, Tgamma G, Tbeta B, Tmu M)
{
  constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
  for (size_t i=0; i<batchSize; i++) {
    if (X[i]<M[i] || G[i] <= 0.0 || B[i] <= 0.0) {
      output[i] = NaN;
    }
    if (X[i] == M[i]) {
      output[i] = (G[i]==1.0)/B[i];
    } 
    else {
      output[i] = 0.0;
    }
  }
  
  if (G.isBatch()) {
    for (size_t i=0; i<batchSize; i++) {
      if (output[i] == 0.0) {
        output[i] = -std::lgamma(G[i]);
      }
    }
  }
  else {
    double gamma = -std::lgamma(G[2019]);
    for (size_t i=0; i<batchSize; i++) {
      if (output[i] == 0.0) {
        output[i] = gamma;
      }
    }
  }
  
  for (size_t i=0; i<batchSize; i++) {
    if (X[i] != M[i]) {
      const double invBeta = 1/B[i];
      double arg = (X[i]-M[i])*invBeta;
      output[i] -= arg;
      arg = _rf_fast_log(arg);
      output[i] += arg*(G[i]-1);
      output[i] = _rf_fast_exp(output[i]);
      output[i] *= invBeta;
    }
  }
}
};

RooSpan<double> RooGamma::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
  using namespace BatchHelpers;
<commit_msg>Fix windows warning in RooGamma<commit_after> /*****************************************************************************
  * Project: RooFit                                                           *
  *                                                                           *
  * Simple Gamma distribution
  * authors: <NAME>, <NAME>
  *                                                                           *
  *****************************************************************************/

/** \class RooGamma
    \ingroup Roofit

Implementation of the Gamma PDF for RooFit/RooStats.
\f[
f(x) = \frac{(x-\mu)^{\gamma-1} \cdot \exp^{(-(x-mu) / \beta}}{\Gamma(\gamma) \cdot \beta^{\gamma}}
\f]
defined for \f$ x \geq 0 \f$ if \f$ \mu = 0 \f$

Notes from <NAME>:

Wikipedia and several sources refer to the Gamma distribution as

\f[
G(\mu,\alpha,\beta) = \beta^\alpha \mu^{(\alpha-1)} \frac{e^{(-\beta \mu)}}{\Gamma(\alpha)}
\f]

Below is the correspondence:

| Wikipedia       | This Implementation      |
|-----------------|--------------------------|
| \f$ \alpha \f$  | \f$ \gamma \f$           |
| \f$ \beta \f$   | \f$ \frac{1}{\beta} \f$  |
| \f$ \mu \f$     | x                        |
| 0               | \f$ \mu \f$              |


Note, that for a model Pois(N|mu), a uniform prior on mu, and a measurement N
the posterior is in the Wikipedia parameterization Gamma(mu, alpha=N+1, beta=1)
thus for this implementation it is:

`RooGamma(_x=mu,_gamma=N+1,_beta=1,_mu=0)`

Also note, that in this case it is equivalent to
RooPoison(N,mu) and treating the function as a PDF in mu.
**/

#include "RooGamma.h"

#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooRandom.h"
#include "RooMath.h"
#include "RooHelpers.h"
#include "BatchHelpers.h"
#include "RooVDTHeaders.h"

#include "TMath.h"
#include <Math/SpecFuncMathCore.h>
#include <Math/PdfFuncMathCore.h>
#include <Math/ProbFuncMathCore.h>

#include <iostream>
#include <cmath>
using namespace std;

ClassImp(RooGamma);

////////////////////////////////////////////////////////////////////////////////

RooGamma::RooGamma(const char *name, const char *title,
          RooAbsReal& _x, RooAbsReal& _gamma,
          RooAbsReal& _beta, RooAbsReal& _mu) :
  RooAbsPdf(name,title),
  x("x","Observable",this,_x),
  gamma("gamma","Mean",this,_gamma),
  beta("beta","Width",this,_beta),
  mu("mu","Para",this,_mu)
{
  RooHelpers::checkRangeOfParameters(this, {&_gamma, &_beta}, 0.);
}

////////////////////////////////////////////////////////////////////////////////

RooGamma::RooGamma(const RooGamma& other, const char* name) :
  RooAbsPdf(other,name), x("x",this,other.x), gamma("mean",this,other.gamma),
  beta("beta",this,other.beta), mu("mu",this,other.mu)
{
}

////////////////////////////////////////////////////////////////////////////////

Double_t RooGamma::evaluate() const
{
  return TMath::GammaDist(x, gamma, mu, beta) ;
}

////////////////////////////////////////////////////////////////////////////////

namespace {
//Author: <NAME>, CERN 22 August 2019

template<class Tx, class Tgamma, class Tbeta, class Tmu>
void compute(	size_t batchSize,
              double * __restrict output,
              Tx X, Tgamma G, Tbeta B, Tmu M)
{
  constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
  for (size_t i=0; i<batchSize; i++) {
    if (X[i]<M[i] || G[i] <= 0.0 || B[i] <= 0.0) {
      output[i] = NaN;
    }
    if (X[i] == M[i]) {
      output[i] = ((G[i]==1.0) ? 1. : 0.)/B[i];
    }
    else {
      output[i] = 0.0;
    }
  }

  if (G.isBatch()) {
    for (size_t i=0; i<batchSize; i++) {
      if (output[i] == 0.0) {
        output[i] = -std::lgamma(G[i]);
      }
    }
  }
  else {
    double gamma = -std::lgamma(G[2019]);
    for (size_t i=0; i<batchSize; i++) {
      if (output[i] == 0.0) {
        output[i] = gamma;
      }
    }
  }

  for (size_t i=0; i<batchSize; i++) {
    if (X[i] != M[i]) {
      const double invBeta = 1/B[i];
      double arg = (X[i]-M[i])*invBeta;
      output[i] -= arg;
      arg = _rf_fast_log(arg);
      output[i] += arg*(G[i]-1);
      output[i] = _rf_fast_exp(output[i]);
      output[i] *= invBeta;
    }
  }
}
};

RooSpan<double> RooGamma::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
  using namespace BatchHelpers;

{'fraction_non_alphanumeric': 0.27563653352020556, 'fraction_numerical': 0.009577201588413922, 'mean_word_length': 1.9165531335149864, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6515362', 'n_tokens_mistral': 2976, 'n_tokens_neox': 2854, 'n_words': 664}