EvtPdf.hh

Go to the documentation of this file.

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
*                                                                      *
* This file is part of EvtGen.                                         *
*                                                                      *
* EvtGen is free software: you can redistribute it and/or modify       *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or    *
* (at your option) any later version.                                  *
*                                                                      *
* EvtGen is distributed in the hope that it will be useful,            *
* but WITHOUT ANY WARRANTY; without even the implied warranty of       *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
* GNU General Public License for more details.                         *
*                                                                      *
* You should have received a copy of the GNU General Public License    *
* along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
***********************************************************************/

#ifndef EVT_PDF_HH
#define EVT_PDF_HH

#include "EvtGenBase/EvtMacros.hh"
#include "EvtGenBase/EvtPdfMax.hh"
#include "EvtGenBase/EvtPredGen.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtStreamInputIterator.hh"
#include "EvtGenBase/EvtValError.hh"

#include <assert.h>
#include <stdio.h>

/*
 *  All classes are templated on the point type T
 *
 * EvtPdf:
 *
 * Probability density function defined on an interval of phase-space.
 * Integral over the interval can be calculated by Monte Carlo integration.
 * Some (but not all) PDFs are analytic in the sense that they can be integrated
 * by numeric quadrature and distributions can be generated according to them.
 *
 * EvtPdfGen:
 *
 * Generator adaptor. Can be used to generate random points
 * distributed according to the PDF for analytic PDFs.
 *
 * EvtPdfPred:
 *
 * Predicate adaptor for PDFs. Can be used for generating random points distributed
 * according to the PDF for any PDF using rejection method. (See "Numerical Recipes").
 *
 * EvtPdfUnary:
 *
 * Adapter for generic algorithms. Evaluates the PDF and returns the value
 *
 * EvtPdfDiv:
 *
 * PDF obtained by division of one PDF by another. Because the two PDFs are
 * arbitrary this PDF is not analytic. When importance sampling is used the
 * original PDF is divided by the analytic comparison function. EvtPdfDiv is
 * used to represent the modified PDF.
 */

template <class T>
class EvtPdfPred;
template <class T>
class EvtPdfGen;

template <class T>
class EvtPdf {
  public:
    EvtPdf() {}
    EvtPdf( const EvtPdf& other ) : _itg( other._itg ) {}
    virtual ~EvtPdf() {}
    virtual EvtPdf<T>* clone() const = 0;

    double evaluate( const T& p ) const
    {
        if ( p.isValid() )
            return pdf( p );
        else
            return 0.;
    }

    // Find PDF maximum. Points are sampled according to pc

    EvtPdfMax<T> findMax( const EvtPdf<T>& pc, int N );

    // Find generation efficiency.

    EvtValError findGenEff( const EvtPdf<T>& pc, int N, int nFindMax );

    // Analytic integration. Calls cascade down until an overridden
    // method is called.

    void setItg( EvtValError itg ) { _itg = itg; }

    EvtValError getItg() const
    {
        if ( !_itg.valueKnown() )
            _itg = compute_integral();
        return _itg;
    }
    EvtValError getItg( int N ) const
    {
        if ( !_itg.valueKnown() )
            _itg = compute_integral( N );
        return _itg;
    }

    virtual EvtValError compute_integral() const
    {
        printf( "Analytic integration of PDF is not defined\n" );
        assert( 0 );
        return EvtValError{};
    }
    virtual EvtValError compute_integral( int ) const
    {
        return compute_integral();
    }

    //  Monte Carlo integration.

    EvtValError compute_mc_integral( const EvtPdf<T>& pc, int N );

    // Generation. Create predicate accept-reject generators.
    // nMax iterations will be used to find the maximum of the accept-reject predicate

    EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> accRejGen( const EvtPdf<T>& pc,
                                                       int nMax,
                                                       double factor = 1. );

    virtual T randomPoint();

  protected:
    virtual double pdf( const T& ) const = 0;
    mutable EvtValError _itg;
};

template <class T>
class EvtPdfGen {
  public:
    typedef T result_type;

    EvtPdfGen() : _pdf( 0 ) {}
    EvtPdfGen( const EvtPdfGen<T>& other ) :
        _pdf( other._pdf ? other._pdf->clone() : 0 )
    {
    }
    EvtPdfGen( const EvtPdf<T>& pdf ) : _pdf( pdf.clone() ) {}
    ~EvtPdfGen() { delete _pdf; }

    result_type operator()() { return _pdf->randomPoint(); }

  private:
    EvtPdf<T>* _pdf;
};

template <class T>
class EvtPdfPred {
  public:
    typedef T argument_type;
    typedef bool result_type;

    EvtPdfPred() {}
    EvtPdfPred( const EvtPdf<T>& thePdf ) : itsPdf( thePdf.clone() ) {}
    EvtPdfPred( const EvtPdfPred& other ) :
        COPY_PTR( itsPdf ), COPY_MEM( itsPdfMax )
    {
    }
    ~EvtPdfPred() { delete itsPdf; }

    result_type operator()( argument_type p )
    {
        assert( itsPdf );
        assert( itsPdfMax.valueKnown() );

        double random = EvtRandom::Flat( 0., itsPdfMax.value() );
        return ( random <= itsPdf->evaluate( p ) );
    }

    EvtPdfMax<T> getMax() const { return itsPdfMax; }
    void setMax( const EvtPdfMax<T>& max ) { itsPdfMax = max; }
    template <class InputIterator>
    void compute_max( InputIterator it, InputIterator end, double factor = 1. )
    {
        T p = *it++;
        itsPdfMax = EvtPdfMax<T>( p, itsPdf->evaluate( p ) * factor );

        while ( !( it == end ) ) {
            T p = *it++;
            double val = itsPdf->evaluate( p ) * factor;
            if ( val > itsPdfMax.value() )
                itsPdfMax = EvtPdfMax<T>( p, val );
        }
    }

  private:
    EvtPdf<T>* itsPdf;
    EvtPdfMax<T> itsPdfMax;
};

template <class T>
class EvtPdfUnary {
  public:
    typedef double result_type;
    typedef T argument_type;

    EvtPdfUnary() {}
    EvtPdfUnary( const EvtPdf<T>& thePdf ) : itsPdf( thePdf.clone() ) {}
    EvtPdfUnary( const EvtPdfUnary& other ) : COPY_PTR( itsPdf ) {}
    ~EvtPdfUnary() { delete itsPdf; }

    result_type operator()( argument_type p )
    {
        assert( itsPdf );
        double ret = itsPdf->evaluate( p );
        return ret;
    }

  private:
    EvtPdf<T>* itsPdf;
};

template <class T>
class EvtPdfDiv : public EvtPdf<T> {
  public:
    EvtPdfDiv() : itsNum( 0 ), itsDen( 0 ) {}
    EvtPdfDiv( const EvtPdf<T>& theNum, const EvtPdf<T>& theDen ) :
        EvtPdf<T>(), itsNum( theNum.clone() ), itsDen( theDen.clone() )
    {
    }
    EvtPdfDiv( const EvtPdfDiv<T>& other ) :
        EvtPdf<T>( other ), COPY_PTR( itsNum ), COPY_PTR( itsDen )
    {
    }
    virtual ~EvtPdfDiv()
    {
        delete itsNum;
        delete itsDen;
    }
    EvtPdf<T>* clone() const override { return new EvtPdfDiv( *this ); }

    double pdf( const T& p ) const override
    {
        double num = itsNum->evaluate( p );
        double den = itsDen->evaluate( p );
        assert( den != 0 );
        return num / den;
    }

  private:
    EvtPdf<T>* itsNum;    // numerator
    EvtPdf<T>* itsDen;    // denominator
};

template <class T>
EvtPdfMax<T> EvtPdf<T>::findMax( const EvtPdf<T>& pc, int N )
{
    EvtPdfPred<T> pred( *this );
    EvtPdfGen<T> gen( pc );
    pred.compute_max( iter( gen, N ), iter( gen ) );
    EvtPdfMax<T> p = pred.getMax();
    return p;
}

template <class T>
EvtValError EvtPdf<T>::findGenEff( const EvtPdf<T>& pc, int N, int nFindMax )
{
    assert( N > 0 || nFindMax > 0 );
    EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> gen = accRejGen( pc, nFindMax );
    int i;
    for ( i = 0; i < N; i++ )
        gen();
    double eff = double( gen.getPassed() ) / double( gen.getTried() );
    double err = sqrt( double( gen.getPassed() ) ) / double( gen.getTried() );
    return EvtValError( eff, err );
}

template <class T>
EvtValError EvtPdf<T>::compute_mc_integral( const EvtPdf<T>& pc, int N )
{
    assert( N > 0 );

    EvtPdfDiv<T> pdfdiv( *this, pc );
    EvtPdfUnary<T> unary( pdfdiv );

    EvtPdfGen<T> gen( pc );
    EvtStreamInputIterator<T> begin = iter( gen, N );
    EvtStreamInputIterator<T> end;

    double sum = 0.;
    double sum2 = 0.;
    while ( !( begin == end ) ) {
        double value = pdfdiv.evaluate( *begin++ );
        sum += value;
        sum2 += value * value;
    }

    EvtValError x;
    if ( N > 0 ) {
        double av = sum / ( (double)N );
        if ( N > 1 ) {
            double dev2 = ( sum2 - av * av * N ) / ( (double)( N - 1 ) );
            // Due to numerical precision dev2 may sometimes be negative
            if ( dev2 < 0. )
                dev2 = 0.;
            double error = sqrt( dev2 / ( (double)N ) );
            x = EvtValError( av, error );
        } else
            x = EvtValError( av );
    }
    _itg = x * pc.getItg();
    return _itg;
}

template <class T>
T EvtPdf<T>::randomPoint()
{
    printf( "Function defined for analytic PDFs only\n" );
    assert( 0 );
    T temp;
    return temp;
}

template <class T>
EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> EvtPdf<T>::accRejGen( const EvtPdf<T>& pc,
                                                              int nMax,
                                                              double factor )
{
    EvtPdfGen<T> gen( pc );
    EvtPdfDiv<T> pdfdiv( *this, pc );
    EvtPdfPred<T> pred( pdfdiv );
    pred.compute_max( iter( gen, nMax ), iter( gen ), factor );
    return EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>>( gen, pred );
}

#endif