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_distGeom](@ref pm_distGeom).
      19             : !>
      20             : !>  \finmain
      21             : !>
      22             : !>  \author
      23             : !>  \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
      24             : 
      25             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      26             : 
      27             :         !%%%%%%%%%%%%%%%%%%%%
      28             : #if     getGeomLogPMF_ENABLED
      29             :         !%%%%%%%%%%%%%%%%%%%%
      30             : 
      31          50 :         CHECK_ASSERTION(__LINE__, 0._RKC < probSuccess, SK_"@setGeomLogPMF(): The condition `0 <= probSuccess` must hold. probSuccess = "//getStr(probSuccess)) ! fpp
      32          50 :         CHECK_ASSERTION(__LINE__, probSuccess <= 1._RKC, SK_"@setGeomLogPMF(): The condition `probSuccess <= 1` must hold. probSuccess = "//getStr(probSuccess)) ! fpp
      33          50 :         call setGeomLogPMF(logPMF, stepSuccess, log(probSuccess))
      34             : 
      35             :         !%%%%%%%%%%%%%%%%%%%%
      36             : #elif   setGeomLogPMF_ENABLED
      37             :         !%%%%%%%%%%%%%%%%%%%%
      38             : 
      39             :         ! Validate the input.
      40         100 :         CHECK_ASSERTION(__LINE__, 0_IK < stepSuccess, SK_"@setGeomLogPMF(): The condition `0 < stepSuccess` must hold. stepSuccess = "//getStr(stepSuccess)) ! fpp
      41         100 :         CHECK_ASSERTION(__LINE__, logProbSuccess <= 0._RKC, SK_"@setGeomLogPMF(): The condition `logProbSuccess <= 0.` must hold. logProbSuccess = "//getStr(logProbSuccess)) ! fpp
      42             : #if     Log_ENABLED
      43           9 :         CHECK_ASSERTION(__LINE__, abs(1._RKC - exp(logProbSuccess) - exp(logProbFailure)) < epsilon(0._RKC) * 100, \
      44             :         SK_"@setGeomLogPMF(): The condition `exp(logProbFailure) + exp(logProbSuccess) == 1.` must hold. logProbFailure, logProbSuccess = "\
      45             :         //getStr([logProbFailure, logProbSuccess])) ! fpp
      46             : #elif   Def_ENABLED
      47             : #define logProbFailure log(get1mexp(logProbSuccess))
      48             : #else
      49             : #error  "Unrecognized interface."
      50             : #endif
      51             :         ! Compute the PMF.
      52         100 :         if (logProbSuccess /= 0._RKC) then ! imperfect probability of success.
      53          98 :             logPMF = (real(stepSuccess, RKC) - 1_IK) * logProbFailure + logProbSuccess
      54             :         else ! 100% probability of success.
      55           2 :             if (1_IK < stepSuccess) then
      56           2 :                 logPMF = -huge(logPMF)
      57             :             else
      58           0 :                 logPMF = 0._RKC
      59             :             end if
      60             :         end if
      61             : 
      62             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      63             : #elif   getGeomCDF_ENABLED || setGeomCDF_ENABLED
      64             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      65             : 
      66             :         ! Validate the input.
      67         100 :         CHECK_ASSERTION(__LINE__, 0._RKC < probSuccess, SK_"@setGeomCDF(): The condition `0. < probSuccess` must hold. probSuccess = "//getStr(probSuccess)) ! fpp
      68         100 :         CHECK_ASSERTION(__LINE__, probSuccess <= 1._RKC, SK_"@setGeomCDF(): The condition `probSuccess <= 1.` must hold. probSuccess = "//getStr(probSuccess)) ! fpp
      69         100 :         CHECK_ASSERTION(__LINE__, 0_IK < stepSuccess, SK_"@setGeomCDF(): The condition `0 <= stepSuccess` must hold. stepSuccess = "//getStr(stepSuccess)) ! fpp
      70         100 :         cdf = 1._RKC - (1._RKC - probSuccess)**stepSuccess
      71             : 
      72             :         !%%%%%%%%%%%%%%%%%%
      73             : #elif   getGeomRand_ENABLED
      74             :         !%%%%%%%%%%%%%%%%%%
      75             : 
      76       20006 :         CHECK_ASSERTION(__LINE__, 0._RKC < probSuccess, SK_"@getGeomRand(): The condition `0 <= probSuccess` must hold. probSuccess = "//getStr(probSuccess)) ! fpp
      77       20006 :         if (probSuccess < 1._RKC) then
      78       20005 :             call setGeomRand(rand, log(1._RKC - probSuccess))
      79             :         else
      80           1 :             rand = 1_IK
      81           1 :             CHECK_ASSERTION(__LINE__, probSuccess <= 1._RKC, SK_"@getGeomRand(): The condition `probSuccess <= 1` must hold. probSuccess = "//getStr(probSuccess)) ! fpp
      82             :         end if
      83             : 
      84             :         !%%%%%%%%%%%%%%%%%%
      85             : #elif   setGeomRand_ENABLED
      86             :         !%%%%%%%%%%%%%%%%%%
      87             : 
      88             : #if     D0_ENABLED
      89             :         real(RKC) :: unifrnd
      90             : #elif   D1_ENABLED
      91           8 :         real(RKC) :: unifrnd(size(rand, 1, IK)), logProbFailureInverse
      92             : #else
      93             : #error  "Unrecognized interface."
      94             : #endif
      95       82022 :         CHECK_ASSERTION(__LINE__, logProbFailure < 0._RKC, SK_"@setGeomRand(): The condition `logProbFailure < 0.` must hold. logProbFailure = "//getStr(logProbFailure)) ! fpp
      96             : #if     RNGD_ENABLED || RNGF_ENABLED
      97       82020 :         call random_number(unifrnd)
      98             : #elif   RNGX_ENABLED
      99           2 :         call setUnifRand(rng, unifrnd)
     100             : #else
     101             : #error  "Unrecognized interface."
     102             : #endif
     103             : #if     D0_ENABLED
     104       82018 :         rand = 1_IK + floor(log(1._RKC - unifrnd) / logProbFailure)
     105             : #elif   D1_ENABLED
     106           4 :         logProbFailureInverse = 1._RKC / logProbFailure
     107          12 :         rand = 1_IK + floor(log(1._RKC - unifrnd) * logProbFailureInverse)
     108             : #endif
     109             : 
     110             : #else
     111             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     112             : #error  "Unrecognized interface."
     113             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     114             : #endif
     115             : #undef  logProbFailure