include/EvtGenBase/EvtPdf.hh

0001
0002 /***********************************************************************
0003 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
0004 *                                                                      *
0005 * This file is part of EvtGen.                                         *
0006 *                                                                      *
0007 * EvtGen is free software: you can redistribute it and/or modify       *
0008 * it under the terms of the GNU General Public License as published by *
0009 * the Free Software Foundation, either version 3 of the License, or    *
0010 * (at your option) any later version.                                  *
0011 *                                                                      *
0012 * EvtGen is distributed in the hope that it will be useful,            *
0013 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
0014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
0015 * GNU General Public License for more details.                         *
0016 *                                                                      *
0017 * You should have received a copy of the GNU General Public License    *
0018 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
0019 ***********************************************************************/
0020
0021 #ifndef EVT_PDF_HH
0022 #define EVT_PDF_HH
0023
0024 #include "EvtGenBase/EvtMacros.hh"
0025 #include "EvtGenBase/EvtPdfMax.hh"
0026 #include "EvtGenBase/EvtPredGen.hh"
0027 #include "EvtGenBase/EvtRandom.hh"
0028 #include "EvtGenBase/EvtStreamInputIterator.hh"
0029 #include "EvtGenBase/EvtValError.hh"
0030
0031 #include <assert.h>
0032 #include <stdio.h>
0033
0034 /*
0035  *  All classes are templated on the point type T
0036  *
0037  * EvtPdf:
0038  *
0039  * Probability density function defined on an interval of phase-space.
0040  * Integral over the interval can be calculated by Monte Carlo integration.
0041  * Some (but not all) PDFs are analytic in the sense that they can be integrated
0042  * by numeric quadrature and distributions can be generated according to them.
0043  *
0044  * EvtPdfGen:
0045  *
0046  * Generator adaptor. Can be used to generate random points
0047  * distributed according to the PDF for analytic PDFs.
0048  *
0049  * EvtPdfPred:
0050  *
0051  * Predicate adaptor for PDFs. Can be used for generating random points distributed
0052  * according to the PDF for any PDF using rejection method. (See "Numerical Recipes").
0053  *
0054  * EvtPdfUnary:
0055  *
0056  * Adapter for generic algorithms. Evaluates the PDF and returns the value
0057  *
0058  * EvtPdfDiv:
0059  *
0060  * PDF obtained by division of one PDF by another. Because the two PDFs are
0061  * arbitrary this PDF is not analytic. When importance sampling is used the
0062  * original PDF is divided by the analytic comparison function. EvtPdfDiv is
0063  * used to represent the modified PDF.
0064  */
0065
0066 template <class T>
0067 class EvtPdfPred;
0068 template <class T>
0069 class EvtPdfGen;
0070
0071 template <class T>
0072 class EvtPdf {
0073   public:
0074     EvtPdf() {}
0075     EvtPdf( const EvtPdf& other ) : _itg( other._itg ) {}
0076     virtual ~EvtPdf() {}
0077     virtual EvtPdf<T>* clone() const = 0;
0078
0079     double evaluate( const T& p ) const
0080     {
0081         if ( p.isValid() )
0082             return pdf( p );
0083         else
0084             return 0.;
0085     }
0086
0087     // Find PDF maximum. Points are sampled according to pc
0088
0089     EvtPdfMax<T> findMax( const EvtPdf<T>& pc, int N );
0090
0091     // Find generation efficiency.
0092
0093     EvtValError findGenEff( const EvtPdf<T>& pc, int N, int nFindMax );
0094
0095     // Analytic integration. Calls cascade down until an overridden
0096     // method is called.
0097
0098     void setItg( EvtValError itg ) { _itg = itg; }
0099
0100     EvtValError getItg() const
0101     {
0102         if ( !_itg.valueKnown() )
0103             _itg = compute_integral();
0104         return _itg;
0105     }
0106     EvtValError getItg( int N ) const
0107     {
0108         if ( !_itg.valueKnown() )
0109             _itg = compute_integral( N );
0110         return _itg;
0111     }
0112
0113     virtual EvtValError compute_integral() const
0114     {
0115         printf( "Analytic integration of PDF is not defined\n" );
0116         assert( 0 );
0117         return EvtValError{};
0118     }
0119     virtual EvtValError compute_integral( int ) const
0120     {
0121         return compute_integral();
0122     }
0123
0124     //  Monte Carlo integration.
0125
0126     EvtValError compute_mc_integral( const EvtPdf<T>& pc, int N );
0127
0128     // Generation. Create predicate accept-reject generators.
0129     // nMax iterations will be used to find the maximum of the accept-reject predicate
0130
0131     EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> accRejGen( const EvtPdf<T>& pc,
0132                                                        int nMax,
0133                                                        double factor = 1. );
0134
0135     virtual T randomPoint();
0136
0137   protected:
0138     virtual double pdf( const T& ) const = 0;
0139     mutable EvtValError _itg;
0140 };
0141
0142 template <class T>
0143 class EvtPdfGen {
0144   public:
0145     typedef T result_type;
0146
0147     EvtPdfGen() : _pdf( 0 ) {}
0148     EvtPdfGen( const EvtPdfGen<T>& other ) :
0149         _pdf( other._pdf ? other._pdf->clone() : 0 )
0150     {
0151     }
0152     EvtPdfGen( const EvtPdf<T>& pdf ) : _pdf( pdf.clone() ) {}
0153     ~EvtPdfGen() { delete _pdf; }
0154
0155     result_type operator()() { return _pdf->randomPoint(); }
0156
0157   private:
0158     EvtPdf<T>* _pdf;
0159 };
0160
0161 template <class T>
0162 class EvtPdfPred {
0163   public:
0164     typedef T argument_type;
0165     typedef bool result_type;
0166
0167     EvtPdfPred() {}
0168     EvtPdfPred( const EvtPdf<T>& thePdf ) : itsPdf( thePdf.clone() ) {}
0169     EvtPdfPred( const EvtPdfPred& other ) :
0170         COPY_PTR( itsPdf ), COPY_MEM( itsPdfMax )
0171     {
0172     }
0173     ~EvtPdfPred() { delete itsPdf; }
0174
0175     result_type operator()( argument_type p )
0176     {
0177         assert( itsPdf );
0178         assert( itsPdfMax.valueKnown() );
0179
0180         double random = EvtRandom::Flat( 0., itsPdfMax.value() );
0181         return ( random <= itsPdf->evaluate( p ) );
0182     }
0183
0184     EvtPdfMax<T> getMax() const { return itsPdfMax; }
0185     void setMax( const EvtPdfMax<T>& max ) { itsPdfMax = max; }
0186     template <class InputIterator>
0187     void compute_max( InputIterator it, InputIterator end, double factor = 1. )
0188     {
0189         T p = *it++;
0190         itsPdfMax = EvtPdfMax<T>( p, itsPdf->evaluate( p ) * factor );
0191
0192         while ( !( it == end ) ) {
0193             T p = *it++;
0194             double val = itsPdf->evaluate( p ) * factor;
0195             if ( val > itsPdfMax.value() )
0196                 itsPdfMax = EvtPdfMax<T>( p, val );
0197         }
0198     }
0199
0200   private:
0201     EvtPdf<T>* itsPdf;
0202     EvtPdfMax<T> itsPdfMax;
0203 };
0204
0205 template <class T>
0206 class EvtPdfUnary {
0207   public:
0208     typedef double result_type;
0209     typedef T argument_type;
0210
0211     EvtPdfUnary() {}
0212     EvtPdfUnary( const EvtPdf<T>& thePdf ) : itsPdf( thePdf.clone() ) {}
0213     EvtPdfUnary( const EvtPdfUnary& other ) : COPY_PTR( itsPdf ) {}
0214     ~EvtPdfUnary() { delete itsPdf; }
0215
0216     result_type operator()( argument_type p )
0217     {
0218         assert( itsPdf );
0219         double ret = itsPdf->evaluate( p );
0220         return ret;
0221     }
0222
0223   private:
0224     EvtPdf<T>* itsPdf;
0225 };
0226
0227 template <class T>
0228 class EvtPdfDiv : public EvtPdf<T> {
0229   public:
0230     EvtPdfDiv() : itsNum( 0 ), itsDen( 0 ) {}
0231     EvtPdfDiv( const EvtPdf<T>& theNum, const EvtPdf<T>& theDen ) :
0232         EvtPdf<T>(), itsNum( theNum.clone() ), itsDen( theDen.clone() )
0233     {
0234     }
0235     EvtPdfDiv( const EvtPdfDiv<T>& other ) :
0236         EvtPdf<T>( other ), COPY_PTR( itsNum ), COPY_PTR( itsDen )
0237     {
0238     }
0239     virtual ~EvtPdfDiv()
0240     {
0241         delete itsNum;
0242         delete itsDen;
0243     }
0244     EvtPdf<T>* clone() const override { return new EvtPdfDiv( *this ); }
0245
0246     double pdf( const T& p ) const override
0247     {
0248         double num = itsNum->evaluate( p );
0249         double den = itsDen->evaluate( p );
0250         assert( den != 0 );
0251         return num / den;
0252     }
0253
0254   private:
0255     EvtPdf<T>* itsNum;    // numerator
0256     EvtPdf<T>* itsDen;    // denominator
0257 };
0258
0259 template <class T>
0260 EvtPdfMax<T> EvtPdf<T>::findMax( const EvtPdf<T>& pc, int N )
0261 {
0262     EvtPdfPred<T> pred( *this );
0263     EvtPdfGen<T> gen( pc );
0264     pred.compute_max( iter( gen, N ), iter( gen ) );
0265     EvtPdfMax<T> p = pred.getMax();
0266     return p;
0267 }
0268
0269 template <class T>
0270 EvtValError EvtPdf<T>::findGenEff( const EvtPdf<T>& pc, int N, int nFindMax )
0271 {
0272     assert( N > 0 || nFindMax > 0 );
0273     EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> gen = accRejGen( pc, nFindMax );
0274     int i;
0275     for ( i = 0; i < N; i++ )
0276         gen();
0277     double eff = double( gen.getPassed() ) / double( gen.getTried() );
0278     double err = sqrt( double( gen.getPassed() ) ) / double( gen.getTried() );
0279     return EvtValError( eff, err );
0280 }
0281
0282 template <class T>
0283 EvtValError EvtPdf<T>::compute_mc_integral( const EvtPdf<T>& pc, int N )
0284 {
0285     assert( N > 0 );
0286
0287     EvtPdfDiv<T> pdfdiv( *this, pc );
0288     EvtPdfUnary<T> unary( pdfdiv );
0289
0290     EvtPdfGen<T> gen( pc );
0291     EvtStreamInputIterator<T> begin = iter( gen, N );
0292     EvtStreamInputIterator<T> end;
0293
0294     double sum = 0.;
0295     double sum2 = 0.;
0296     while ( !( begin == end ) ) {
0297         double value = pdfdiv.evaluate( *begin++ );
0298         sum += value;
0299         sum2 += value * value;
0300     }
0301
0302     EvtValError x;
0303     if ( N > 0 ) {
0304         double av = sum / ( (double)N );
0305         if ( N > 1 ) {
0306             double dev2 = ( sum2 - av * av * N ) / ( (double)( N - 1 ) );
0307             // Due to numerical precision dev2 may sometimes be negative
0308             if ( dev2 < 0. )
0309                 dev2 = 0.;
0310             double error = sqrt( dev2 / ( (double)N ) );
0311             x = EvtValError( av, error );
0312         } else
0313             x = EvtValError( av );
0314     }
0315     _itg = x * pc.getItg();
0316     return _itg;
0317 }
0318
0319 template <class T>
0320 T EvtPdf<T>::randomPoint()
0321 {
0322     printf( "Function defined for analytic PDFs only\n" );
0323     assert( 0 );
0324     T temp;
0325     return temp;
0326 }
0327
0328 template <class T>
0329 EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>> EvtPdf<T>::accRejGen( const EvtPdf<T>& pc,
0330                                                               int nMax,
0331                                                               double factor )
0332 {
0333     EvtPdfGen<T> gen( pc );
0334     EvtPdfDiv<T> pdfdiv( *this, pc );
0335     EvtPdfPred<T> pred( pdfdiv );
0336     pred.compute_max( iter( gen, nMax ), iter( gen ), factor );
0337     return EvtPredGen<EvtPdfGen<T>, EvtPdfPred<T>>( gen, pred );
0338 }
0339
0340 #endif