The ParaMonte Documentation Website
Current view: top level - kernel - Math_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 229 233 98.3 %
Date: 2021-01-08 13:07:16 Functions: 27 27 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !>  \brief This module contains mathematical procedures.
      44             : !>  \author Amir Shahmoradi
      45             : 
      46             : module Math_mod
      47             : 
      48             :     implicit none
      49             : 
      50             :     character(*), parameter :: MODULE_NAME = "@Math_mod"
      51             : 
      52             :     interface getCumSum
      53             :         module procedure :: getCumSum_IK, getCumSum_RK
      54             :     end interface getCumSum
      55             : 
      56             :     interface getCumSumReverse
      57             :         module procedure :: getCumSumReverse_IK, getCumSumReverse_RK
      58             :     end interface getCumSumReverse
      59             : 
      60             :     interface getLogSumExp
      61             :         module procedure :: getLogSumExp_RK, getLogSumExp_CK
      62             :     end interface getLogSumExp
      63             : 
      64             :     interface getLogSubExp
      65             :         module procedure :: getLogSubExp_RK
      66             :     end interface getLogSubExp
      67             : 
      68             :     interface getLogEggBox
      69             :         module procedure :: getLogEggBoxSD_RK, getLogEggBoxMD_RK
      70             :         module procedure :: getLogEggBoxSD_CK, getLogEggBoxMD_CK
      71             :     end interface getLogEggBox
      72             : 
      73             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      74             : 
      75             : contains
      76             : 
      77             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      78             : 
      79             :     !> \brief
      80             :     !> Return the distance squared between the two input points.
      81             :     !>
      82             :     !> @param[in]   nd      : The size of the input vectors.
      83             :     !> @param[in]   Point1  : The first point.
      84             :     !> @param[in]   Point2  : The second point.
      85             :     !>
      86             :     !> \return
      87             :     !> `distanceSq` : The distance squared.
      88           3 :     pure function getDistanceSq(nd,Point1,Point2) result(distanceSq)
      89             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      90             :         !DEC$ ATTRIBUTES DLLEXPORT :: getDistanceSq
      91             : #endif
      92             :         use Constants_mod, only: IK, RK
      93             :         integer(IK) , intent(in)    :: nd
      94             :         real(RK)    , intent(in)    :: Point1(nd),Point2(nd)
      95             :         real(RK)                    :: distanceSq
      96             :         integer(IK)                 :: i
      97           3 :         distanceSq = 0._RK
      98          18 :         do i = 1, nd
      99          18 :             distanceSq = distanceSq + (Point2(i)-Point1(i))**2
     100             :         end do
     101           3 :     end function getDistanceSq
     102             : 
     103             : 
     104             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     105             : 
     106             :     !> \brief
     107             :     !> Return the correlation coefficient (`-1 < corCoef < 1`) corresponding to the input Fisher z-transformation.
     108             :     !>
     109             :     !> @param[in]   fisherTrans : The Fisher transformation.
     110             :     !>
     111             :     !> \return
     112             :     !> `corCoef` : The correlation coefficient corresponding to the input Fisher transformation.
     113           3 :     pure elemental function getCorCeofFromFisherTrans(fisherTrans) result(corCoef)
     114             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     115             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCorCeofFromFisherTrans
     116             : #endif
     117           3 :         use Constants_mod, only: RK
     118             :         real(RK), intent(in) :: fisherTrans
     119             :         real(RK)             :: corCoef !, twiceFisherTrans
     120           3 :         corCoef = tanh(fisherTrans)
     121             :         !twiceFisherTrans = 2._RK * corCoef
     122             :         !fisherTrans = (exp(twiceFisherTrans)-1._RK) / (exp(twiceFisherTrans)+1._RK)
     123           6 :     end function getCorCeofFromFisherTrans
     124             : 
     125             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     126             : 
     127             :     !> \brief
     128             :     !> Return Fisher z-transformation of an input correlation coefficient (`-1 < corCoef < 1`).
     129             :     !>
     130             :     !> @param[in]   corCoef : The input correlation coefficient.
     131             :     !>
     132             :     !> \return
     133             :     !> `fisherTrans` : The output Fisher transformation.
     134           3 :     pure elemental function getFisherTransFromCorCoef(corCoef) result(fisherTrans)
     135             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     136             :         !DEC$ ATTRIBUTES DLLEXPORT :: getFisherTransFromCorCoef
     137             : #endif
     138           3 :         use Constants_mod, only: RK
     139             :         real(RK), intent(in) :: corCoef
     140             :         real(RK)             :: fisherTrans
     141           3 :         fisherTrans = atanh(corCoef)
     142             :         !fisherTrans = 0.5_RK * log(1._RK+corCoef) / (1._RK-corCoef)
     143           6 :     end function getFisherTransFromCorCoef
     144             : 
     145             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     146             : 
     147             :     !> \brief
     148             :     !> Return the cumulative sum of the input integer array.
     149             :     !>
     150             :     !> @param[in]   vecLen  : The length of the input vector.
     151             :     !> @param[in]   Vec     : The input vector.
     152             :     !>
     153             :     !> \return
     154             :     !> `CumSum` : An integer array of length `vecLen` representing the cumulative sum.
     155             :     !>
     156             :     !> \remark
     157             :     !> The first element of `CumSum` is the same as the first element of `Vec`.
     158        2916 :     pure function getCumSum_IK(vecLen,Vec) result(CumSum)
     159             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     160             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCumSum_IK
     161             : #endif
     162           3 :         use Constants_mod, only: IK
     163             :         integer(IK), intent(in) :: vecLen
     164             :         integer(IK), intent(in) :: Vec(vecLen)
     165             :         integer(IK)             :: CumSum(vecLen)
     166             :         integer(IK)             :: i
     167        1458 :         CumSum(1) = Vec(1)
     168      351944 :         do i = 2, vecLen
     169      351944 :             CumSum(i) = CumSum(i-1) + Vec(i)
     170             :         end do
     171        1458 :     end function getCumSum_IK
     172             : 
     173             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     174             : 
     175             :     !> \brief
     176             :     !> Return the cumulative sum of the input real array.
     177             :     !> @param[in]   vecLen  : The length of the input vector.
     178             :     !> @param[in]   Vec     : The input vector.
     179             :     !>
     180             :     !> \return
     181             :     !> `CumSum` : A real array of length `vecLen` representing the cumulative sum.
     182             :     !>
     183             :     !> \remark
     184             :     !> The first element of `CumSum` is the same as the first element of `Vec`.
     185     1175220 :     pure function getCumSum_RK(vecLen,Vec) result(CumSum)
     186             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     187             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCumSum_RK
     188             : #endif
     189        1458 :         use Constants_mod, only: IK, RK
     190             :         integer(IK), intent(in) :: vecLen
     191             :         real(RK)   , intent(in) :: Vec(vecLen)
     192             :         real(RK)                :: CumSum(vecLen)
     193             :         integer(IK)             :: i
     194         306 :         CumSum(1) = Vec(1)
     195     1174300 :         do i = 2, vecLen
     196     1174300 :             CumSum(i) = CumSum(i-1) + Vec(i)
     197             :         end do
     198         306 :     end function getCumSum_RK
     199             : 
     200             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     201             : 
     202             :     !> \brief
     203             :     !> Return the cumulative sum of the input integer array, in the backward direction.
     204             :     !> @param[in]   vecLen  : The length of the input vector.
     205             :     !> @param[in]   Vec     : The input vector.
     206             :     !>
     207             :     !> \return
     208             :     !> `CumSumReverse` : An integer array of length `vecLen` representing the cumulative sum.
     209             :     !>
     210             :     !> \remark
     211             :     !> The last element of `CumSumReverse` is the same as the last element of `Vec`.
     212           6 :     pure function getCumSumReverse_IK(vecLen,Vec) result(CumSumReverse)
     213             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     214             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCumSumReverse_IK
     215             : #endif
     216         306 :         use Constants_mod, only: IK
     217             :         integer(IK), intent(in) :: vecLen
     218             :         integer(IK), intent(in) :: Vec(vecLen)
     219             :         integer(IK)             :: CumSumReverse(vecLen)
     220             :         integer(IK)             :: i, indx
     221           3 :         CumSumReverse(1) = Vec(vecLen)
     222          30 :         do i = vecLen-1,1,-1
     223          27 :             indx = vecLen - i
     224          30 :             CumSumReverse(indx+1) = CumSumReverse(indx) + Vec(i)
     225             :         end do
     226           3 :     end function getCumSumReverse_IK
     227             : 
     228             :     !> \brief
     229             :     !> Return the cumulative sum of the input real array, in the backward direction.
     230             :     !>
     231             :     !> @param[in]   vecLen  : The length of the input vector.
     232             :     !> @param[in]   Vec     : The input vector.
     233             :     !>
     234             :     !> \return
     235             :     !> `CumSumReverse` : A real array of length `vecLen` representing the cumulative sum.
     236             :     !>
     237             :     !> \remark
     238             :     !> The last element of `CumSumReverse` is the same as the last element of `Vec`.
     239          39 :     pure function getCumSumReverse_RK(vecLen,Vec) result(CumSumReverse)
     240             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     241             :         !DEC$ ATTRIBUTES DLLEXPORT :: getCumSumReverse_RK
     242             : #endif
     243           3 :         use Constants_mod, only: IK, RK
     244             :         integer(IK), intent(in) :: vecLen
     245             :         real(RK)   , intent(in) :: Vec(vecLen)
     246             :         real(RK)                :: CumSumReverse(vecLen)
     247             :         integer(IK)             :: i, indx
     248           3 :         CumSumReverse(1) = Vec(vecLen)
     249          30 :         do i = vecLen-1,1,-1
     250          27 :             indx = vecLen - i
     251          30 :             CumSumReverse(indx+1) = CumSumReverse(indx) + Vec(i)
     252             :         end do
     253           3 :     end function getCumSumReverse_RK
     254             : 
     255             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     256             : 
     257             :     !> \brief
     258             :     !> Return `log( exp(logValueLarger) - exp(logValueSamller) )` robustly (without overflow or underflow).
     259             :     !>
     260             :     !> @param[in]   logValueLarger  : A real number.
     261             :     !> @param[in]   logValueSamller : A real number.
     262             :     !>
     263             :     !> \return
     264             :     !> `logSubExp` : A real number.
     265             :     !>
     266             :     !> \warning
     267             :     !> The onus is on the user to ensure `logValueLarger > logValueSamller`.
     268             :     !>
     269             :     !> \remark
     270             :     !> This function is very useful for situations where `exp(logValueLarger)` is likely to cause overflow.
     271        8261 :     pure function getLogSubExp_RK(logValueLarger,logValueSamller) result(logSubExp)
     272             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     273             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogSubExp_RK
     274             : #endif
     275           3 :         use Constants_mod, only: RK, LOGTINY_RK
     276             :         real(RK)   , intent(in) :: logValueLarger, logValueSamller
     277        8261 :         real(RK)                :: logSubExp, diff
     278        8261 :         logSubExp = logValueLarger
     279        8261 :         diff = logValueSamller - logValueLarger
     280        8260 :         if (diff>LOGTINY_RK) logSubExp = logSubExp + log( 1._RK - exp(diff) )
     281        8261 :     end function getLogSubExp_RK
     282             : 
     283             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     284             : 
     285             :     !> \brief
     286             :     !> Return the logarithm of the sum of the exponential of the input real vector robustly (without overflow or underflow).
     287             :     !>
     288             :     !> @param[in]   lenLogValue : The length of the input vector.
     289             :     !> @param[in]   LogValue    : The input vector of log-values whose log-sum-exp must be computed.
     290             :     !> @param[in]   maxLogValue : The maximum of the input `LogValue` argument (**optional**).
     291             :     !>
     292             :     !> \return
     293             :     !> `logSumExp` : A real number.
     294           6 :     pure function getLogSumExp_RK(lenLogValue, LogValue, maxLogValue) result(logSumExp)
     295             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     296             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogSumExp_RK
     297             : #endif
     298        8261 :         use Constants_mod, only: IK, RK, LOGTINY_RK
     299             :         integer(IK) , intent(in)            :: lenLogValue
     300             :         real(RK)    , intent(in)            :: LogValue(lenLogValue)
     301             :         real(RK)    , intent(in), optional  :: maxLogValue
     302             :         real(RK)                            :: logSumExp
     303           6 :         real(RK)                            :: maxLogValueDefault
     304          24 :         real(RK)                            :: LogValueCopy(lenLogValue)
     305             :         integer(IK)                         :: i
     306           6 :         if (present(maxLogValue)) then
     307           3 :             maxLogValueDefault = maxLogValue
     308             :         else
     309          15 :             maxLogValueDefault = maxval(LogValue)
     310             :         end if
     311          24 :         LogValueCopy = LogValue - maxLogValueDefault
     312           6 :         do concurrent(i=1:lenLogValue)
     313          24 :             if (LogValueCopy(i)<LOGTINY_RK) then
     314           0 :                 LogValueCopy(i) = 0._RK
     315             :             else
     316          18 :                 LogValueCopy(i) = exp(LogValueCopy(i))
     317             :             end if
     318             :         end do
     319          24 :         logSumExp = maxLogValueDefault + log(sum(LogValueCopy))
     320           6 :     end function getLogSumExp_RK
     321             : 
     322             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     323             : 
     324             :     !> \brief
     325             :     !> Return the logarithm of the sum of the exponential of the input complex vector robustly (without overflow or underflow).
     326             :     !>
     327             :     !> @param[in]   lenLogValue : The length of the input vector.
     328             :     !> @param[in]   LogValue    : The input vector of log-values whose log-sum-exp must be computed.
     329             :     !> @param[in]   maxLogValue : The maximum of the real component of the input `LogValue` argument (**optional**).
     330             :     !>
     331             :     !> \return
     332             :     !> `logSumExp` : A complex number.
     333           6 :     pure function getLogSumExp_CK(lenLogValue, LogValue, maxLogValue) result(logSumExp)
     334             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     335             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogSumExp_CK
     336             : #endif
     337           6 :         use Constants_mod, only: IK, RK, LOGTINY_RK
     338             :         integer(IK) , intent(in)            :: lenLogValue
     339             :         complex(RK) , intent(in)            :: LogValue(lenLogValue)
     340             :         complex(RK) , intent(in), optional  :: maxLogValue
     341             :         complex(RK)                         :: logSumExp
     342          24 :         complex(RK)                         :: LogValueCopy(lenLogValue)
     343             :         complex(RK)                         :: maxLogValueDefault
     344             :         integer(IK)                         :: i
     345           6 :         if (present(maxLogValue)) then
     346           3 :             maxLogValueDefault = maxLogValue
     347             :         else
     348          15 :             maxLogValueDefault = maxval(real(LogValue))
     349             :         end if
     350          24 :         LogValueCopy = LogValue - maxLogValueDefault
     351           6 :         do concurrent(i=1:lenLogValue)
     352          24 :             if (real(LogValueCopy(i))<LOGTINY_RK) then
     353           0 :                 LogValueCopy(i) = 0._RK
     354             :             else
     355          18 :                 LogValueCopy(i) = exp(LogValueCopy(i))
     356             :             end if
     357             :         end do
     358          24 :         logSumExp = maxLogValueDefault + log(sum(LogValueCopy))
     359           6 :     end function getLogSumExp_CK
     360             : 
     361             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     362             : 
     363             :     !> \brief
     364             :     !> Return the logarithm of the egg-box probability density function in one dimension.
     365           3 :     pure function getLogEggBoxSD_RK(constant,exponent,coef,point) result(logEggBox)
     366             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     367             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxSD_RK
     368             : #endif
     369           6 :         use Constants_mod, only: IK, RK
     370             :         implicit none
     371             :         real(RK), intent(in) :: constant
     372             :         real(RK), intent(in) :: exponent
     373             :         real(RK), intent(in) :: coef
     374             :         real(RK), intent(in) :: point
     375             :         real(RK)             :: logEggBox
     376           3 :         logEggBox = exponent * ( constant + cos(coef*point) )
     377           6 :     end function getLogEggBoxSD_RK
     378             : 
     379             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     380             : 
     381             :     !> \brief
     382             :     !> Return the logarithm of the egg-box probability density function in one dimension, as a complex number.
     383           3 :     pure function getLogEggBoxSD_CK(constant,exponent,coef,point) result(logEggBox)
     384             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     385             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxSD_CK
     386             : #endif
     387           3 :         use Constants_mod, only: IK, RK
     388             :         implicit none
     389             :         complex(RK), intent(in) :: constant
     390             :         complex(RK), intent(in) :: exponent
     391             :         complex(RK), intent(in) :: coef
     392             :         complex(RK), intent(in) :: Point
     393             :         complex(RK)             :: logEggBox
     394           3 :         logEggBox = exponent * ( constant + cos(coef*point) )
     395           6 :     end function getLogEggBoxSD_CK
     396             : 
     397             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     398             : 
     399             :     !> \brief
     400             :     !> Return the logarithm of the egg-box probability density function in multiple dimensions, as a real number.
     401           3 :     pure function getLogEggBoxMD_RK(nd,constant,exponent,Coef,Point) result(logEggBox)
     402             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     403             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxMD_RK
     404             : #endif
     405           3 :         use Constants_mod, only: IK, RK
     406             :         implicit none
     407             :         integer(IK), intent(in) :: nd
     408             :         real(RK), intent(in) :: constant
     409             :         real(RK), intent(in) :: exponent
     410             :         real(RK), intent(in) :: Coef(nd)
     411             :         real(RK), intent(in) :: Point(nd)
     412             :         real(RK)             :: logEggBox
     413             :         integer(IK)          :: i
     414           3 :         logEggBox = 1._RK
     415          12 :         do i = 1, nd
     416          12 :             logEggBox = logEggBox * cos(Coef(i)*Point(i))
     417             :         end do
     418           3 :         logEggBox = exponent * ( constant + logEggBox )
     419           3 :     end function getLogEggBoxMD_RK
     420             : 
     421             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     422             : 
     423             :     !> \brief
     424             :     !> Return the logarithm of the egg-box probability density function in multiple dimensions, as a complex number.
     425             :     !> \remark
     426             :     !> The multidimensional EggBox function follows an equation of the form,
     427             :     !> \f{equation}
     428             :     !> f({\mathbf{\theta}}) = \bigg[ \exp( \text{constant} + \prod_{i=1)^{\text{nd}} ~ \cos(\text{Coef}_i^2)
     429             :     !> \f}
     430           3 :     pure function getLogEggBoxMD_CK(nd,constant,exponent,Coef,Point) result(logEggBox)
     431             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     432             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEggBoxMD_CK
     433             : #endif
     434           3 :         use Constants_mod, only: IK, RK
     435             :         implicit none
     436             :         integer(IK), intent(in) :: nd
     437             :         complex(RK), intent(in) :: constant
     438             :         complex(RK), intent(in) :: exponent
     439             :         complex(RK), intent(in) :: Coef(nd)
     440             :         complex(RK), intent(in) :: Point(nd)
     441             :         complex(RK)             :: logEggBox
     442             :         integer(IK)          :: i
     443           3 :         logEggBox = 1._RK
     444          12 :         do i = 1, nd
     445          12 :             logEggBox = logEggBox * cos(Coef(i)*Point(i))
     446             :         end do
     447           3 :         logEggBox = exponent * ( constant + logEggBox )
     448           3 :     end function getLogEggBoxMD_CK
     449             : 
     450             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     451             : 
     452             :     !> \brief
     453             :     !> Return the log(factorial) for a whole integer input. This is basically `logGamma( positiveInteger + 1 )`.
     454             :     !>
     455             :     !> \param[in]   positiveInteger :   The input positive integer whose log factorial must be computed.
     456             :     !>
     457             :     !> \return
     458             :     !> `logFactorial`  :   The natural logarithm of the factorial.
     459             :     !>
     460             :     !> \remark
     461             :     !> This function is mostly useful for large input integers.
     462         123 :     pure function getLogFactorial(positiveInteger) result(logFactorial)
     463             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     464             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogFactorial
     465             : #endif
     466             :         use Constants_mod, only: IK, RK ! LCOV_EXCL_LINE
     467             :         implicit none
     468             :         integer(IK), intent(in) :: positiveInteger
     469             :         integer(IK)             :: i
     470             :         real(RK)                :: logFactorial
     471         123 :         logFactorial = 0._RK
     472         174 :         do i = 2, positiveInteger
     473         174 :             logFactorial = logFactorial + log(real(i,kind=RK))
     474             :         end do
     475         123 :     end function getLogFactorial
     476             : 
     477             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     478             : 
     479             :     !> \brief
     480             :     !> Return the factorial for a whole integer input. This is basically `Gamma( intNum + 1 )`.
     481             :     !>
     482             :     !> \param[in]   positiveInteger :   The input positive integer whose log factorial must be computed.
     483             :     !>
     484             :     !> \return
     485             :     !> `factorial`  :   The factorial.
     486           3 :     pure function getFactorial(positiveInteger) result(factorial)
     487             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     488             :         !DEC$ ATTRIBUTES DLLEXPORT :: getFactorial
     489             : #endif
     490         123 :         use Constants_mod, only: IK, RK
     491             :         implicit none
     492             :         integer(IK), intent(in) :: positiveInteger
     493             :         integer(IK)             :: i
     494             :         real(RK)                :: factorial
     495           3 :         factorial = 1._RK
     496          30 :         do i = 2, positiveInteger
     497          30 :             factorial = factorial * i
     498             :         end do
     499           3 :     end function getFactorial
     500             : 
     501             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     502             : 
     503             :     !> \brief
     504             :     !> Return the Gamma function for a half-integer input as real of kind `RK`.
     505             :     !>
     506             :     !> \param[in]   positiveHalfInteger :   The input half integer as a real number
     507             :     !>
     508             :     !> \remark
     509             :     !> The equation for half-integer Gamma-function is given as,
     510             :     !> \f{equation}
     511             :     !>     \Gamma \left( \frac{n}{2} \right) = \sqrt \pi \frac{ (n-2)!! }{ 2^\frac{n-1}{2} } ~,
     512             :     !> \f}
     513             :     !> 
     514             :     !> \return
     515             :     !> `gammaHalfInt`   :   The Gamma function for a half integer input.
     516           3 :     pure function getGammaHalfInt(positiveHalfInteger) result(gammaHalfInt)
     517             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     518             :         !DEC$ ATTRIBUTES DLLEXPORT :: getGammaHalfInt
     519             : #endif
     520           3 :         use Constants_mod, only: IK, RK, SQRTPI
     521             :         implicit none
     522             :         real(RK), intent(in) :: positiveHalfInteger
     523             :         real(RK)             :: gammaHalfInt
     524             :         integer(IK)          :: i,k
     525           3 :         gammaHalfInt = SQRTPI
     526           3 :         k = nint(positiveHalfInteger-0.5_RK,kind=IK) ! positiveHalfInteger = k + 1/2
     527           3 :         do i = k+1, 2*k
     528           3 :             gammaHalfInt = gammaHalfInt * 0.25_RK * i
     529             :         end do
     530           3 :     end function getGammaHalfInt
     531             : 
     532             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     533             : 
     534             :     !> \brief
     535             :     !> Return the natural logarithm of the Gamma function for a half-integer input as real of kind `RK`.
     536             :     !>
     537             :     !> \param[in]   positiveHalfInteger :   The input half integer as a real number
     538             :     !>
     539             :     !> \remark
     540             :     !> The equation for half-integer Gamma-function is given as,
     541             :     !> \f{equation}
     542             :     !>     \Gamma \left( \frac{n}{2} \right) = \sqrt \pi \frac{ (n-2)!! }{ 2^\frac{n-1}{2} } ~,
     543             :     !> \f}
     544             :     !> 
     545             :     !> \return
     546             :     !> `gammaHalfInt`   :   The Gamma function for a half integer input.
     547           3 :     pure function getLogGammaHalfInt(positiveHalfInteger) result(logGammaHalfInt)
     548             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     549             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogGammaHalfInt
     550             : #endif
     551           3 :         use Constants_mod, only: IK, RK, SQRTPI
     552             :         implicit none
     553             :         real(RK), intent(in)    :: positiveHalfInteger
     554             :         real(RK), parameter     :: COEF = log(0.25_RK)
     555             :         real(RK), parameter     :: LOG_SQRTPI = log(SQRTPI)
     556             :         real(RK)                :: logGammaHalfInt
     557             :         integer(IK)             :: i, k
     558           3 :         k = nint(positiveHalfInteger-0.5_RK,kind=IK) ! positiveHalfInteger = k + 1/2
     559           3 :         logGammaHalfInt = LOG_SQRTPI
     560           3 :         do i = k+1, 2*k
     561           3 :             logGammaHalfInt = logGammaHalfInt + COEF + log(real(i,kind=RK))
     562             :         end do
     563           3 :     end function getLogGammaHalfInt
     564             : 
     565             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     566             : 
     567             :     ! see for example: http://math.stackexchange.com/questions/606184/volume-of-n-dimensional-ellipsoid
     568             :     ! see for example: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
     569             :     ! see for example: https://en.wikipedia.org/wiki/Particular_values_of_the_Gamma_function
     570             :     ! getEllVolCoef = PI^(nd/2) / gamma(nd/2+1) where n is just a positive integer
     571             :     !
     572             :     !> \brief
     573             :     !> Return the coefficient of the volume of an `nd`-dimensional ellipsoid.
     574             :     !>
     575             :     !> param[in]    nd : The number of dimensions.
     576             :     !>
     577             :     !> \return
     578             :     !> `ellVolCoef` : The coefficient of the volume of an `nd`-dimensional ellipsoid.
     579           6 :     pure function getEllVolCoef(nd) result(ellVolCoef)
     580             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     581             :         !DEC$ ATTRIBUTES DLLEXPORT :: getEllVolCoef
     582             : #endif
     583           3 :         use Constants_mod, only: IK, RK, PI
     584             :         implicit none
     585             :         integer(IK), intent(in) :: nd
     586             :         integer(IK)             :: i,k
     587             :         real(RK), parameter     :: FOUR_PI = PI * 4._RK
     588             :         real(RK)                :: ellVolCoef
     589           6 :         if (mod(nd,2_IK)==0_IK) then ! nd is even
     590             : 
     591           3 :             ellVolCoef = PI
     592          15 :             do i = 2_IK, nd / 2_IK
     593          15 :                 ellVolCoef = ellVolCoef * PI / i    ! nd = 2k ; ellVolCoef = PI^(k) / Factorial(k)
     594             :             end do
     595             : 
     596             :         else ! nd is an odd integer
     597             : 
     598             :             ! nd = 2k-1 ; gamma(nd/2 + 1) = gamma(k + 1/2) ; gamma(k+1/2) = sqrt(PI) * (2k)! / (4^k * k!)
     599           3 :             k = (nd + 1_IK) / 2_IK
     600             : 
     601             :             ! This is to avoid an extra unnecessary division of ellVolCoef by PI
     602           3 :             ellVolCoef = 4._RK / (k + 1_IK)
     603             : 
     604          18 :             do i = k+2_IK, 2_IK*k
     605             :                 ! ellVolCoef = PI^(k-1/2) / gamma(k+1/2) = PI^(k+1) * 4^k * k! / (2k)!
     606          18 :                 ellVolCoef = ellVolCoef * FOUR_PI / i
     607             :             end do
     608             : 
     609             :         end if
     610           6 :     end function getEllVolCoef
     611             : 
     612             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     613             : 
     614             :     ! see for example: http://math.stackexchange.com/questions/606184/volume-of-n-dimensional-ellipsoid
     615             :     ! see for example: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
     616             :     ! see for example: https://en.wikipedia.org/wiki/Particular_values_of_the_Gamma_function
     617             :     ! getEllVolCoef = PI^(nd/2) / gamma(nd/2+1) where n is just a positive integer
     618             :     !
     619             :     !> \brief
     620             :     !> Return the coefficient of the volume of an `nd`-dimensional ellipsoid.
     621             :     !>
     622             :     !> param[in]    nd : The number of dimensions.
     623             :     !>
     624             :     !> \return
     625             :     !> `logEllVolCoef` : The natural logarithm of the volume of an `nd`-dimensional unit-ball.
     626           6 :     pure function getLogEllVolCoef(nd) result(logEllVolCoef)
     627             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     628             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogEllVolCoef
     629             : #endif
     630           6 :         use Constants_mod, only: IK, RK, PI
     631             :         implicit none
     632             :         integer(IK), intent(in) :: nd
     633             :         integer(IK)             :: i,k
     634             :         real(RK), parameter     :: FOUR_PI = PI * 4._RK
     635             :         real(RK)                :: logEllVolCoef
     636           6 :         if (mod(nd,2_IK)==0_IK) then ! nd is even
     637             : 
     638           3 :             logEllVolCoef = PI
     639          15 :             do i = 2_IK, nd / 2_IK
     640          15 :                 logEllVolCoef = logEllVolCoef * PI / i    ! nd = 2k ; ellVolCoef = PI^(k) / Factorial(k)
     641             :             end do
     642             : 
     643             :         else ! nd is an odd integer
     644             : 
     645             :             ! nd = 2k-1 ; gamma(nd/2 + 1) = gamma(k + 1/2) ; gamma(k+1/2) = sqrt(PI) * (2k)! / (4^k * k!)
     646           3 :             k = (nd + 1_IK) / 2_IK
     647             : 
     648             :             ! This is to avoid an extra unnecessary division of ellVolCoef by PI
     649           3 :             logEllVolCoef = 4._RK / (k + 1_IK)
     650             : 
     651          18 :             do i = k+2_IK, 2_IK*k
     652             :                 ! ellVolCoef = PI^(k-1/2) / gamma(k+1/2) = PI^(k+1) * 4^k * k! / (2k)!
     653          18 :                 logEllVolCoef = logEllVolCoef * FOUR_PI / i
     654             :             end do
     655             : 
     656             :         end if
     657           6 :         logEllVolCoef = log(logEllVolCoef)
     658           6 :     end function getLogEllVolCoef
     659             : 
     660             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     661             : 
     662             :     !> \brief
     663             :     !> Return the logarithm of the volume of an `nd`-dimensional ball of unit-radius.
     664             :     !>
     665             :     !> param[in]    nd  : The number of dimensions.
     666             :     !>
     667             :     !> \return
     668             :     !> `logVolUnitBall` : The logarithm of the volume of an `nd`-dimensional ball of unit-radius.
     669             :     !>
     670             :     !> \brief
     671             :     !> This routine is an exact replacement for [getLogEllVolCoef](@ref getlogellvolcoef).
     672             :     !> However, it is, on average, 10%-15% faster than [getLogEllVolCoef](@ref getlogellvolcoef) based
     673             :     !> on benchmarks with Intel ifort 19.4 with all optimization flags on.
     674             :     !
     675             :     ! see for example: http://math.stackexchange.com/questions/606184/volume-of-n-dimensional-ellipsoid
     676             :     ! see for example: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
     677             :     ! see for example: https://en.wikipedia.org/wiki/Particular_values_of_the_Gamma_function
     678             :     ! getEllVolCoef = PI^(nd/2) / gamma(nd/2+1) where n is just a positive integer
     679         123 :     pure function getLogVolUnitBall(nd) result(logVolUnitBall)
     680             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     681             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogVolUnitBall
     682             : #endif
     683           6 :         use Constants_mod, only: IK, RK, PI
     684             :         implicit none
     685             :         integer(IK) , intent(in)    :: nd
     686             :         real(RK)    , parameter     :: LOG_PI = log(PI)
     687             :         integer(IK)                 :: ndHalfInteger
     688             :         real(RK)                    :: logVolUnitBall
     689         123 :         real(RK)                    :: ndHalfReal
     690         123 :         if (mod(nd,2_IK)==0_IK) then ! nd is even
     691         120 :             ndHalfInteger = nd / 2_IK
     692         120 :             logVolUnitBall = ndHalfInteger*LOG_PI - getLogFactorial(ndHalfInteger) ! nd = 2k ; logVolUnitBall = PI^(k) / Factorial(k)
     693             :         else ! nd is an odd integer
     694           3 :             ndHalfReal = 0.5_RK * real(nd,kind=RK)
     695           3 :             logVolUnitBall = ndHalfReal * LOG_PI - log_gamma(ndHalfReal+1._RK)
     696             :         end if
     697         123 :     end function getLogVolUnitBall
     698             : 
     699             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     700             : 
     701             :     !> \brief
     702             :     !> Return the logarithm of the volume of an `nd`-dimensional hyper-ellipsoid.
     703             :     !>
     704             :     !> param[in]    nd                  :   The number of dimensions.
     705             :     !> param[in]    logSqrtDetCovMat    :   The logarithm of the square root of the determinant of the representative covariance matrix of the ellipsoid.
     706             :     !>
     707             :     !> \return
     708             :     !> `logVolEllipsoid` : The logarithm of the volume of an `nd`-dimensional hyper-ellipsoid.
     709             :     !
     710             :     ! see for example: https://math.stackexchange.com/questions/332391/volume-of-hyperellipsoid/332434
     711             :     ! see for example: https://math.stackexchange.com/questions/2854930/volume-of-an-n-dimensional-ellipsoid
     712           6 :     pure function getLogVolEllipsoid(nd,logSqrtDetCovMat) result(logVolEllipsoid)
     713             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     714             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogVolEllipsoid
     715             : #endif
     716         123 :         use Constants_mod, only: IK, RK
     717             :         implicit none
     718             :         integer(IK) , intent(in)    :: nd
     719             :         real(RK)    , intent(in)    :: logSqrtDetCovMat
     720             :         real(RK)                    :: logVolEllipsoid
     721           6 :         logVolEllipsoid = getLogVolUnitBall(nd) + logSqrtDetCovMat
     722          12 :     end function getLogVolEllipsoid
     723             : 
     724             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     725             : 
     726             :     !> \brief
     727             :     !> Return the volume of a list of `nEllipsoid` `nd`-dimensional hyper-ellipsoids.
     728             :     !>
     729             :     !> @param[in]   nd                  :   The number of dimensions of the ellipsoids.
     730             :     !> @param[in]   nEllipsoid          :   The number of ellipsoids.
     731             :     !> @param[in]   LogSqrtDetCovMat    :   The vector of natural logarithms of the square roots of the determinants of the covariance matrices of the ellipsoids.
     732             :     !>
     733             :     !> \return
     734             :     !> `logVolEllipsoid` : The logarithm of the volume of an `nd`-dimensional hyper-ellipsoid.
     735             :     !>
     736             :     !> \remark
     737             :     !> There is really no reason to use this function for HPC solutions, as it can be likely done more efficiently.
     738          15 :     pure function getLogVolEllipsoids(nd,nEllipsoid,LogSqrtDetCovMat) result(LogVolEllipsoids)
     739             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     740             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogVolEllipsoids
     741             : #endif
     742           6 :         use Constants_mod, only: IK, RK
     743             :         implicit none
     744             :         integer(IK)                 :: i
     745             :         integer(IK) , intent(in)    :: nd
     746             :         integer(IK) , intent(in)    :: nEllipsoid
     747             :         real(RK)    , intent(in)    :: LogSqrtDetCovMat(nEllipsoid)
     748             :         real(RK)                    :: LogVolEllipsoids(nEllipsoid)
     749           3 :         real(RK)                    :: logVolUnitBall
     750           3 :         logVolUnitBall = getLogVolUnitBall(nd)
     751           9 :         do i = 1, nEllipsoid
     752           9 :             LogVolEllipsoids(i) = logVolUnitBall + LogSqrtDetCovMat(i)
     753             :         end do
     754           3 :     end function getLogVolEllipsoids
     755             : 
     756             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     757             : 
     758             :     !> \brief
     759             :     !> Return the lower incomplete gamma function for the specified `exponent` and upper limit.
     760             :     !> `tolerance` represents the relative accuracy.
     761             :     !>
     762             :     !> \warning
     763             :     !> Do not set tolerance to a number larger than 1, or else the code will crash spectacularly.
     764             :     !>
     765             :     !> \remark
     766             :     !> `logGammaExponent = log_gamma(exponent)`
     767             :     !>
     768             :     !> \warning
     769             :     !> On input, both `exponent` and `upperLim` must be positive, otherwise, a `lowerGamma = -infinity` will be returned to signal error.
     770             :     ! This algorithm borrows from the `gammp` implementation of the Numerical Recipes by Press et al 1992.
     771          39 :     pure function getLowerGamma(exponent,logGammaExponent,upperLim,tolerance) result(lowerGamma)
     772             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     773             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLowerGamma
     774             : #endif
     775           3 :         use Constants_mod, only: IK, RK, HUGE_RK
     776             :         implicit none
     777             :         real(RK), intent(in)            :: exponent, logGammaExponent, upperLim
     778             :         real(RK), intent(in), optional  :: tolerance
     779             :         real(RK)                        :: lowerGamma
     780          39 :         if ( upperLim < 0._RK .or. exponent <= 0._RK ) then
     781           3 :             lowerGamma = -HUGE_RK
     782           3 :             return
     783          36 :         elseif ( upperLim < exponent + 1._RK ) then
     784          24 :             lowerGamma = getGammaSeries(exponent,log_gamma(exponent),upperLim,tolerance)
     785             :         else
     786          12 :             lowerGamma = 1._RK - getGammaContFrac(exponent,logGammaExponent,upperLim,tolerance)
     787             :         end if
     788          39 :     end function getLowerGamma
     789             : 
     790             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     791             : 
     792             :     !> \brief
     793             :     !> Return the lower incomplete gamma function for the specified exponent and upper limit.
     794             :     !>
     795             :     !> `tolerance` represents the relative accuracy.
     796             :     !>
     797             :     !> \warning
     798             :     !> Do not set `tolerance` to a number larger than 1, or else the code will crash spectacularly.
     799             :     !>
     800             :     !> \warning
     801             :     !> On input, both `exponent` and `lowerLim` must be positive, otherwise, a `upperGamma = -infinity` will be returned to signal error.
     802             :     ! This algorithm borrows from the gammq implementation of the Numerical Recipes by Press et al 1992.
     803          39 :     pure function getUpperGamma(exponent,logGammaExponent,lowerLim,tolerance) result(upperGamma)
     804             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     805             :         !DEC$ ATTRIBUTES DLLEXPORT :: getUpperGamma
     806             : #endif
     807          39 :         use Constants_mod, only: IK, RK, HUGE_RK
     808             :         implicit none
     809             :         real(RK), intent(in)            :: exponent, logGammaExponent, lowerLim
     810             :         real(RK), intent(in), optional  :: tolerance
     811             :         real(RK)                        :: upperGamma
     812          39 :         if ( lowerLim < 0._RK .or. exponent <= 0._RK ) then
     813           3 :             upperGamma = -HUGE_RK
     814           3 :             return
     815          36 :         elseif (lowerLim < exponent + 1._RK) then
     816          24 :             upperGamma = 1._RK - getGammaSeries(exponent,logGammaExponent,lowerLim,tolerance)
     817             :         else
     818          12 :             upperGamma = getGammaContFrac(exponent,logGammaExponent,lowerLim,tolerance)
     819             :         end if
     820          39 :     end function getUpperGamma
     821             : 
     822             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     823             : 
     824             :     !> \brief
     825             :     !> Return the lower incomplete gamma function `P(exponent, upperLim)` evaluated by its series representation as `gammaSeries`.
     826             :     !>
     827             :     !> \param[in]   exponent            :   The exponent in lower incomplete gamma function `P(exponent, upperLim)`.
     828             :     !> \param[in]   logGammaExponent    :   The input `log( gamma(exponent) )`.
     829             :     !> \param[in]   upperLim            :   The upper limit in lower incomplete gamma function `P(exponent, upperLim)`.
     830             :     !> \param[in]   tolerance           :   The input relative accuracy.
     831             :     !>
     832             :     !> \warning
     833             :     !> Do not set `tolerance` to a number larger than 1, or else the code will crash spectacularly.
     834             :     !>
     835             :     !> \warning
     836             :     !> If the algorithm fails to converge, `gammaSeries` will be set to negative infinity on output to signal error.
     837             :     ! This algorithm borrows from the gser implementation of the Numerical Recipes by Press et al 1992.
     838          51 :     pure function getGammaSeries(exponent,logGammaExponent,upperLim,tolerance) result(gammaSeries)
     839             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     840             :         !DEC$ ATTRIBUTES DLLEXPORT :: getGammaSeries
     841             : #endif
     842          39 :         use Constants_mod, only: IK, RK, HUGE_RK
     843             :         implicit none
     844             :         real(RK), intent(in)            :: exponent, logGammaExponent, upperLim
     845             :         real(RK), intent(in), optional  :: tolerance
     846             :         real(RK)                        :: gammaSeries
     847             :         integer(IK), parameter          :: ITMAX = 100
     848             :         real(RK), parameter             :: EPS_DEFAULT = epsilon(upperLim)
     849          51 :         real(RK)                        :: ap,del,summ,eps
     850             :         integer(IK)                     :: iter
     851          51 :         if (present(tolerance)) then
     852          39 :             eps = tolerance
     853             :         else
     854          12 :             eps = EPS_DEFAULT
     855             :         end if
     856          51 :         if (upperLim == 0._RK) then
     857           3 :             gammaSeries = 0._RK
     858           3 :             return
     859             :         end if
     860          48 :         ap = exponent
     861          48 :         summ = 1.0_RK / exponent
     862          48 :         del = summ
     863         606 :         do iter = 1, ITMAX
     864         606 :             ap = ap + 1.0_RK
     865         606 :             del = del*upperLim / ap
     866         606 :             summ = summ + del
     867         606 :             if (abs(del) < abs(summ)*eps) exit
     868             :         end do
     869          48 :         if (iter>ITMAX) then
     870           0 :             gammaSeries = -HUGE_RK
     871             :         else
     872          48 :             gammaSeries = summ * exp( exponent*log(upperLim) - upperLim - logGammaExponent )
     873             :         end if
     874          51 :     end function getGammaSeries
     875             : 
     876             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     877             : 
     878             :     !> \brief
     879             :     !> Return the incomplete gamma function `Q(exponent, lowerLim)` evaluated by its continued fraction representation as `gammaContFrac`.
     880             :     !>
     881             :     !> @param[in]   logGammaExponent    :   This is the `log( gamma(exponent) )`.
     882             :     !> @param[in]   tolerance           :   Represents the relative accuracy (**optional**).
     883             :     !> @param[in]   exponent            :   The exponent.
     884             :     !> @param[in]   lowerLim            :   The lower limit.
     885             :     !>
     886             :     !> \warning
     887             :     !> If the algorithm fails to converge, `gammaContFrac` will be set to negative infinity on output to signal error.
     888             :     ! This algorithm borrows from the gcf implementation of the Numerical Recipes by Press et al 1992.
     889          27 :     pure function getGammaContFrac(exponent,logGammaExponent,lowerLim,tolerance) result(gammaContFrac)
     890             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     891             :         !DEC$ ATTRIBUTES DLLEXPORT :: getGammaContFrac
     892             : #endif
     893          51 :         use Constants_mod, only: IK, RK, HUGE_RK
     894             :         implicit none
     895             :         real(RK), intent(in)            :: exponent, logGammaExponent, lowerLim
     896             :         real(RK), optional, intent(in)  :: tolerance
     897             :         real(RK)                        :: gammaContFrac
     898             :         integer(IK), parameter          :: ITMAX = 100_IK                               !< The maximum allowed number of iterations.
     899             :         real(RK), parameter             :: EPS_DEFAULT = epsilon(lowerLim)              !< The default relative accuracy.
     900             :         real(RK), parameter             :: FPMIN_DEFAULT = tiny(lowerLim) / EPS_DEFAULT !< A number near the smallest representable floating-point number.
     901          27 :         real(RK)                        :: eps, fpmin
     902          27 :         real(RK)                        :: an,b,c,d,del,h
     903             :         integer(IK)                     :: iter
     904          27 :         if (lowerLim == 0._RK) then
     905           3 :             gammaContFrac = 1._RK
     906           3 :             return
     907             :         end if
     908          24 :         if (present(tolerance)) then
     909          18 :             eps = tolerance
     910          18 :             fpmin = tiny(lowerLim) / eps
     911             :         else
     912           6 :             eps = EPS_DEFAULT
     913           6 :             fpmin = FPMIN_DEFAULT
     914             :         end if
     915          24 :         b = lowerLim + 1._RK - exponent
     916          24 :         c = 1._RK / fpmin
     917          24 :         d = 1._RK / b
     918          24 :         h = d
     919         198 :         do iter = 1, ITMAX
     920         198 :             an = -iter * (iter-exponent)
     921         198 :             b = b + 2.0_RK
     922         198 :             d = an*d + b
     923         198 :             if (abs(d) < fpmin) d = fpmin
     924         198 :             c = b + an/c
     925         198 :             if (abs(c) < fpmin) c = fpmin
     926         198 :             d = 1.0_RK / d
     927         198 :             del = d*c
     928         198 :             h = h*del
     929         198 :             if (abs(del-1._RK) <= eps) exit
     930             :         end do
     931             :        !if (iter > ITMAX) call nrerror('exponent too large, ITMAX too small in getGammaContFrac')
     932          24 :         if (iter>ITMAX) then
     933           0 :             gammaContFrac = -HUGE_RK
     934             :         else
     935          24 :             gammaContFrac = exp(exponent*log(lowerLim) - lowerLim - logGammaExponent)*h
     936             :         end if
     937          27 :     end function getGammaContFrac
     938             : 
     939             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     940             : 
     941             : end module Math_mod ! LCOV_EXCL_LINE

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021