Line data    Source code

       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!                                                                                                                            !!!!
       4             : !!!!    ParaMonte: Parallel Monte Carlo and Machine Learning Library.                                                           !!!!
       5             : !!!!                                                                                                                            !!!!
       6             : !!!!    Copyright (C) 2012-present, The Computational Data Science Lab                                                          !!!!
       7             : !!!!                                                                                                                            !!!!
       8             : !!!!    This file is part of the ParaMonte library.                                                                             !!!!
       9             : !!!!                                                                                                                            !!!!
      10             : !!!!    LICENSE                                                                                                                 !!!!
      11             : !!!!                                                                                                                            !!!!
      12             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md                                                          !!!!
      13             : !!!!                                                                                                                            !!!!
      14             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      15             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      16             : 
      17             : !>  \brief
      18             : !>  This include file contains the implementation of procedures in [pm_distPois](@ref pm_distPois).
      19             : !>
      20             : !>  \finmain
      21             : !>
      22             : !>  \author
      23             : !>  \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
      24             : 
      25             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      26             : 
      27             :         !%%%%%%%%%%%%%%%%%%%%
      28             : #if     getPoisLogPMF_ENABLED
      29             :         !%%%%%%%%%%%%%%%%%%%%
      30             : 
      31          93 :         call setPoisLogPMF(logPMF, count, lambda)
      32             : 
      33             :         !%%%%%%%%%%%%%%%%%%%%
      34             : #elif   setPoisLogPMF_ENABLED
      35             :         !%%%%%%%%%%%%%%%%%%%%
      36             : 
      37             :         ! Validate the input.
      38         187 :         CHECK_ASSERTION(__LINE__, 0_IK <= count, SK_"@setPoisLogPMF(): The condition `0 < count` must hold. count = "//getStr(count)) ! fpp
      39         187 :         CHECK_ASSERTION(__LINE__, 0._RKC < lambda, SK_"@setPoisLogPMF(): The condition `0. < lambda` must hold. lambda = "//getStr(lambda)) ! fpp
      40             : #if     Log_ENABLED
      41           3 :         CHECK_ASSERTION(__LINE__, abs(logLambda - log(lambda)) < epsilon(0._RKC) * 10, \
      42             :         SK_"@setPoisLogPMF(): The condition `logLambda == log(lambda)` must hold. logLambda, log(lambda) = "//getStr([logLambda, log(lambda)])) ! fpp
      43             : #elif   Def_ENABLED
      44             : #define LOGLAMBDA log(lambda)
      45             : #else
      46             : #error  "Unrecognized interface."
      47             : #endif
      48             :         ! Compute the PMF.
      49         187 :         logPMF = real(count, RKC) * LOGLAMBDA - lambda - log_gamma(real(count, RKC) + 1._RKC)
      50             : 
      51             :         !%%%%%%%%%%%%%%%%%
      52             : #elif   getPoisCDF_ENABLED
      53             :         !%%%%%%%%%%%%%%%%%
      54             : 
      55             :         integer(IK) :: info
      56             :         real(RKC) :: countP1
      57          93 :         CHECK_ASSERTION(__LINE__, 0_IK <= count, SK_"@getPoisCDF(): The condition `0 <= count` must hold. count = "//getStr(count)) ! fpp
      58          93 :         countP1 = real(count, RKC) + 1._RKC
      59          93 :         call setPoisCDF(cdf, countP1, log_gamma(countP1), lambda, info)
      60          93 :         if (info < 0_IK) error stop MODULE_NAME//SK_"@getPoisCDF(): Call to setPoisCDF failed."
      61             : 
      62             :         !%%%%%%%%%%%%%%%%%
      63             : #elif   setPoisCDF_ENABLED
      64             :         !%%%%%%%%%%%%%%%%%
      65             : 
      66             :         ! Validate the input.
      67         186 :         CHECK_ASSERTION(__LINE__, 0._RKC < lambda, SK_"@setPoisCDF(): The condition `0. < lambda` must hold. lambda = "//getStr(lambda)) ! fpp
      68         186 :         CHECK_ASSERTION(__LINE__, 1._RKC <= countP1, SK_"@setPoisCDF(): The condition `1._RKC <= countP1` must hold. countP1 = "//getStr(countP1)) ! fpp
      69         558 :         CHECK_ASSERTION(__LINE__, abs(logGammaCountP1 - log_gamma(countP1)) <= epsilon(1._RKC) * 100, \
      70             :         SK_"@setPoisCDF(): The condition `logGammaCountP1 == log(countP1)` must hold. logGammaCountP1, log_gamma(countP1) = "\
      71             :         //getStr([logGammaCountP1, log_gamma(countP1)])) ! fpp
      72         558 :         CHECK_ASSERTION(__LINE__, mod(countP1, 1._RKC) == 0._RKC, \
      73             :         SK_"@setPoisCDF(): The condition `mod(countP1, 1._RKC) == 0._RKC` must hold. countP1, mod(countP1, 1._RKC) = "\
      74             :         //getStr([countP1, mod(countP1, 1._RKC)])) ! fpp
      75         186 :         call setGammaIncUpp(cdf, x = lambda, logGammaKappa = logGammaCountP1, kappa = countP1, info = info, tol = tol)
      76             : 
      77             :         !%%%%%%%%%%%%%%%%%%
      78             : #elif   getPoisRand_ENABLED
      79             :         !%%%%%%%%%%%%%%%%%%
      80             : 
      81       20077 :         if (lambda < 10) then
      82       15003 :             call setPoisRand(rand, exp(-lambda))
      83             :         else
      84        5074 :             call setPoisRand(rand, lambda, log(lambda), sqrt(lambda))
      85             :         end if
      86             : 
      87             :         !%%%%%%%%%%%%%%%%%%
      88             : #elif   setPoisRand_ENABLED
      89             :         !%%%%%%%%%%%%%%%%%%
      90             : 
      91             :         ! Set the URNG.
      92             : #if     RNGD_ENABLED
      93             : #define RNG
      94             : #elif   RNGF_ENABLED || RNGX_ENABLED
      95             : #define RNG rng,
      96             : #else
      97             : #error  "Unrecognized interface."
      98             : #endif
      99             :         ! Set the dimension of `rand`.
     100             : #if     D0_ENABLED
     101             : #define GET_RAND(i) rand
     102             : #elif   D1_ENABLED
     103             : #define GET_RAND(i) rand(i)
     104             :         integer(IK) :: irand
     105             : #else
     106             : #error  "Unrecognized interface."
     107             : #endif
     108             :         real(RKC), parameter :: LAMBDA_LIMIT_RKC = real(LAMBDA_LIMIT, RKC)
     109             : #if     Exp_ENABLED
     110             :         ! Use only when lambda < 10.
     111             :         real(RKC) :: prod, unifrnd
     112       15009 :         CHECK_ASSERTION(__LINE__, expNegLambda < 1._RKC, \
     113             :         SK_"@setPoisRand(): The condition `expNegLambda < 1.` must hold. expNegLambda = "//getStr(expNegLambda)) ! fpp
     114       45027 :         CHECK_ASSERTION(__LINE__, exp(-LAMBDA_LIMIT_RKC) < expNegLambda, \
     115             :         SK_"@setPoisRand(): The condition `exp(-LAMBDA_LIMIT) < expNegLambda` must hold. exp(-LAMBDA_LIMIT), expNegLambda = "\
     116             :         //getStr([exp(-LAMBDA_LIMIT_RKC), expNegLambda])) ! fpp
     117             : #if     D1_ENABLED
     118       15009 :         do irand = 1_IK, size(rand, 1, IK)
     119             : #endif
     120       30008 :             GET_RAND(irand) = 0_IK
     121           0 :             prod = 1._RKC
     122           5 :             do
     123       80568 :                 call setUnifRand(RNG unifrnd)
     124       80568 :                 prod = prod * unifrnd
     125       80568 :                 if (prod <= expNegLambda) exit
     126       50560 :                 GET_RAND(irand) = GET_RAND(irand) + 1_IK
     127             :             end do
     128             : #if     D1_ENABLED
     129             :         end do
     130             : #endif
     131             : #elif   Rej_ENABLED
     132             :         real(RKB)   , parameter :: UR = .43_RKC
     133             :         real(RKC)               :: unifrnd, rndunif, afac, bfac, invAlpha, vr, us, randr
     134       15234 :         CHECK_ASSERTION(__LINE__, LAMBDA_LIMIT_RKC <= lambda, \
     135             :         SK_"@setPoisRand(): The condition `LAMBDA_LIMIT <= lambda` must hold. lambda = "\
     136             :         //getStr([LAMBDA_LIMIT_RKC, lambda])) ! fpp
     137       15234 :         CHECK_ASSERTION(__LINE__, abs(logLambda - log(lambda)) < epsilon(lambda) * 100, \
     138             :         SK_"@setPoisRand(): The condition `logLambda == sqrt(lambda)` must hold. logLambda, log(lambda) = "\
     139             :         //getStr([logLambda, log(lambda)])) ! fpp
     140       15234 :         CHECK_ASSERTION(__LINE__, abs(sqrtLambda - sqrt(lambda)) < epsilon(lambda) * 100, \
     141             :         SK_"@setPoisRand(): The condition `sqrtLambda == sqrt(lambda)` must hold. sqrtLambda, sqrt(lambda) = "\
     142             :         //getStr([sqrtLambda, sqrt(lambda)])) ! fpp
     143        5078 :         bfac = 0.931_RKC + 2.53_RKC * sqrtLambda
     144        5078 :         afac = -0.059_RKC + 0.02483_RKC * bfac
     145        5078 :         invAlpha = 1.1239 + 1.1328_RKC / (bfac - 3.4_RKC)
     146        5077 :         vr = 0.9277_RKC - 3.6224_RKC / (bfac - 2._RKC)
     147             : #if     D1_ENABLED
     148        5001 :         do irand = 1_IK, size(rand, 1, IK)
     149             : #endif
     150           1 :             do
     151       13285 :                 call setUnifRand(RNG unifrnd)
     152       13285 :                 unifrnd = unifrnd - 0.5_RKC
     153       13285 :                 call setUnifRand(RNG rndunif)
     154       13285 :                 us = 0.5_RKC - abs(unifrnd)
     155       13285 :                 GET_RAND(irand) = floor((2 * afac / us + bfac) * unifrnd + lambda + UR, kind = IK)
     156       13285 :                 if (0.07_RKC <= us .and. rndunif <= vr) exit
     157        8274 :                 if (GET_RAND(irand) < 0_IK .or. (us < 0.013_RKC .and. us < rndunif)) cycle
     158        7805 :                 randr = real(GET_RAND(irand), RKC)
     159       12816 :                 if (log(rndunif) + log(invAlpha) - log(afac / (us * us) + bfac) <= randr * logLambda - lambda - log_gamma(randr + 1._RKC)) exit
     160             :             end do
     161             : #if     D1_ENABLED
     162             :         end do
     163             : #endif
     164             : #else
     165             : #error  "Unrecognized interface."
     166             : #endif
     167             : 
     168             : #else
     169             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     170             : #error  "Unrecognized interface."
     171             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     172             : #endif
     173             : #undef  GET_RAND
     174             : #undef  RNG