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_distUnifEll](@ref pm_distUnifEll).
      19             : !>
      20             : !>  \finmain
      21             : !>
      22             : !>  \author
      23             : !>  \AmirShahmoradi, April 23, 2017, 1:36 AM, Institute for Computational Engineering and Sciences (ICES), University of Texas at Austin
      24             : 
      25             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      26             : 
      27             :         ! Set the uniform RNG
      28             : #if     RNGD_ENABLED
      29             : #define RNG
      30             : #elif   RNGF_ENABLED || RNGX_ENABLED
      31             : #define RNG rng,
      32             : #elif   getMUR_ENABLED || setMUR_ENABLED
      33             : #error  "Unrecognized interface."
      34             : #endif
      35             :         !%%%%%%%%%%%%%%%%%%%%%%%
      36             : #if     getUnifEllLogPDF_ENABLED
      37             :         !%%%%%%%%%%%%%%%%%%%%%%%
      38             : 
      39             : #if     D0_ENABLED
      40           4 :         CHECK_ASSERTION(__LINE__, 0_IK < ndim, SK_"@getUnifEllLogPDF(): The condition `0 < ndim` must hold. ndim = "//getStr(ndim))
      41           4 :         logPDF = -ndim * logChoDia - getLogVolUnitBall(real(ndim, RKC))
      42             : #elif   D1_ENABLED
      43           5 :         logPDF = -sum(logChoDia) - getLogVolUnitBall(real(size(logChoDia, 1, IK), RKC))
      44             : #elif   D2_ENABLED && AIP_ENABLED
      45             :         integer(IK) :: info
      46           4 :         real(RKC) :: chol(size(gramian, 1, IK), size(gramian, 2, IK))
      47           2 :         CHECK_ASSERTION(__LINE__, size(gramian, 1, IK) == size(gramian, 2, IK), SK_"@getUnifEllLogPDF(): The condition `size(gramian, 1) == size(gramian, 2)` must hold. shape(gramian) = "//getStr(shape(gramian,IK)))
      48           2 :         call setMatCopy(chol, rdpack, gramian, rdpack, uppDia, doff = 0_IK)
      49           2 :         call setMatDetSqrtLog(chol, uppDia, logPDF, info, chol, transHerm)
      50           2 :         if (info /= 0_IK) error stop "The Cholesky factorization of Gramian failed. Graming is not positive definite."
      51           2 :         logPDF = -(logPDF + getLogVolUnitBall(real(size(gramian, 1, IK), RKC)))
      52             :         !logPDF = -getLogVolUnitBall(real(size(gramian, 1, IK), RKC)) - getMatMulTraceLog(chol)
      53             : #else
      54             : #error  "Unrecognized interface."
      55             : #endif
      56             : 
      57             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%
      58             : #elif   D1_ENABLED && setMUR_ENABLED
      59             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%
      60             : 
      61             :         integer(IK) :: idim, ndim
      62             :         real(RKC)   :: unifrnd, ndimInv
      63             :         real(RKC)   :: sumSqUnifBallRand
      64             : #if     AC_ENABLED
      65             : #define UNIFBALLRAND unifBallRand
      66       44008 :         real(RKC)   :: unifBallRand(size(rand, 1, IK))
      67             : #elif   DC_ENABLED
      68             : #define UNIFBALLRAND rand
      69             : #else
      70             : #error  "Unrecognized interface."
      71             : #endif
      72             :         ! Perform runtime bound checks.
      73             : #if     AM_ENABLED
      74       60012 :         CHECK_ASSERTION(__LINE__, size(rand, 1, IK) == size(mean, 1, IK), SK_"@setUnifEllRand(): The condition `size(rand) == size(mean)` must hold. size(rand), size(mean) = "//getStr([size(rand, 1, IK), size(mean, 1, IK)]))
      75             : #elif   AC_ENABLED
      76       96016 :         CHECK_ASSERTION(__LINE__, all(size(rand, 1, IK) == shape(chol, IK)), SK_"@setUnifEllRand(): The condition `all(size(rand) == shape(chol))` must hold. size(rand), shape(chol) = "//getStr([size(rand, 1, IK), shape(chol, IK)]))
      77             : #elif   !(DM_ENABLED || DC_ENABLED)
      78             : #error  "Unrecognized interface."
      79             : #endif
      80        5001 :         ndim = size(rand, kind = IK)
      81       37007 :         ndimInv = 1._RKC / real(ndim, RKC)
      82             :         do
      83           0 :             sumSqUnifBallRand = 0._RKC
      84      111021 :             do idim = 1_IK, ndim
      85       74014 :                 call setNormRand(RNG UNIFBALLRAND(idim))
      86      111021 :                 sumSqUnifBallRand = sumSqUnifBallRand + UNIFBALLRAND(idim)**2
      87             :             end do
      88             :             ! Ensure the vector is not origin. Highly unlikely but possible.
      89       37007 :             if (0._RKC < sumSqUnifBallRand) exit
      90             :         end do
      91             : #if     RNGD_ENABLED || RNGF_ENABLED
      92       37007 :         call random_number(unifrnd)
      93             : #elif   RNGX_ENABLED
      94           0 :         call setUnifRand(rng, unifrnd)
      95             : #endif
      96       37007 :         unifrnd = unifrnd**ndimInv / sqrt(sumSqUnifBallRand)
      97      111021 :         UNIFBALLRAND = UNIFBALLRAND * unifrnd ! a uniform random point from inside of nd-sphere
      98             : #if     AC_ENABLED
      99             :         ! Define the indexing rules.
     100             : #if     XLD_ENABLED
     101             : #define GET_INDEX(I,J)I,J
     102             : #elif   UXD_ENABLED
     103             : #define GET_INDEX(I,J)J,I
     104             : #elif   D1_ENABLED && setMUR_ENABLED
     105             : #error  "Unrecognized interface."
     106             : #endif
     107             :         ! Define the default mean.
     108             : #if     DM_ENABLED
     109             : #define MEAN_PLUS(I)
     110             : #elif   AM_ENABLED
     111             : #define MEAN_PLUS(I) mean(I) +
     112             : #else
     113             : #error  "Unrecognized interface."
     114             : #endif
     115             :         ! Separate the first to allow the possibility of adding an optional `mean`.
     116       66012 :         rand(1 : ndim) = MEAN_PLUS(1 : ndim) chol(GET_INDEX(1 : ndim, 1)) * unifBallRand(1)
     117       44008 :         do idim = 2_IK, ndim
     118       66012 :             rand(idim : ndim) = rand(idim : ndim) + chol(GET_INDEX(idim : ndim, idim)) * unifBallRand(idim)
     119             :         end do
     120             : #endif
     121             : #undef  UNIFBALLRAND
     122             : #undef  MEAN_PLUS
     123             : #undef  GET_INDEX
     124             : #undef  MEAN
     125             : 
     126             :         !%%%%%%%%%%%%%
     127             : #elif   getMUR_ENABLED
     128             :         !%%%%%%%%%%%%%
     129             : 
     130             : #if     AM_ENABLED && DC_ENABLED
     131        5001 :         call setUnifEllRand(RNG rand, mean)
     132             : #elif   DM_ENABLED && AC_ENABLED
     133        5003 :         call setUnifEllRand(RNG rand, chol, subset)
     134             : #elif   AM_ENABLED && AC_ENABLED
     135        5001 :         call setUnifEllRand(RNG rand, mean, chol, subset)
     136             : #else
     137             : #error  "Unrecognized interface."
     138             : #endif
     139             : 
     140             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%
     141             : #elif   D2_ENABLED && setMUR_ENABLED
     142             :         !%%%%%%%%%%%%%%%%%%%%%%%%%%%
     143             : 
     144             :         integer(IK) :: ipnt, ndim
     145           2 :         ndim = size(rand, 1, IK)
     146        2002 :         do ipnt = 1_IK, size(rand, 2, IK)
     147             : #if         DM_ENABLED && DC_ENABLED
     148           0 :             call setUnifEllRand(RNG rand(1:ndim, ipnt))
     149             : #elif       AM_ENABLED && DC_ENABLED
     150           0 :             call setUnifEllRand(RNG rand(1:ndim, ipnt), mean)
     151             : #elif       DM_ENABLED && AC_ENABLED
     152        2002 :             call setUnifEllRand(RNG rand(1:ndim, ipnt), chol, subset)
     153             : #elif       AM_ENABLED && AC_ENABLED
     154           0 :             call setUnifEllRand(RNG rand(1:ndim, ipnt), mean, chol, subset)
     155             : #else
     156             : #error  "Unrecognized interface."
     157             : #endif
     158             :         end do
     159             : 
     160             : #else
     161             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     162             : #error  "Unrecognized interface."
     163             :         !%%%%%%%%%%%%%%%%%%%%%%%%
     164             : #endif
     165             : #undef  RNG