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 file contains procedure implementations of [pm_ellipsoid](@ref pm_ellipsoid).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Wednesday 00:57 AM, September 22, 2021, Dallas TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if (getLogVolUnitBall_ENABLED || setLogVolUnitBall_ENABLED) && Gamm_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 : real(RKC) :: ndimHalf
32 : real(RKC), parameter :: LOG_PI = log(acos(-1._RKC))
33 2026 : CHECK_ASSERTION(__LINE__, real(0, RKC) <= ndim, \
34 : SK_"@setLogVolUnitBall(): The condition `0 <= ndim` must hold. ndim = "//getStr(ndim))
35 2026 : ndimHalf = 0.5_RKC * ndim
36 2026 : if (0._RKC < ndim) then
37 2022 : logVolUnitBall = ndimHalf * LOG_PI - log_gamma(1._RKC + ndimHalf)
38 : else
39 2 : logVolUnitBall = 0._RKC
40 : end if
41 :
42 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 : #elif (getLogVolUnitBall_ENABLED || setLogVolUnitBall_ENABLED) && Iter_ENABLED
44 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45 :
46 : integer(IK) :: idim, hdim
47 : real(RKC) , parameter :: PI = acos(-1._RKC)
48 : real(RKC) , parameter :: FOUR_PI = PI * 4._RKC
49 6 : CHECK_ASSERTION(__LINE__, 0_IK <= ndim, \
50 : SK_"@setLogVolUnitBall(): The condition `0 <= ndim` must hold. ndim = "//getStr(ndim))
51 6 : if (0_IK < ndim) then
52 5 : hdim = ndim / 2_IK
53 5 : if (ndim == 2_IK * hdim) then
54 : ! ndim is even
55 : ! ndim = 2 * hdim
56 : ! logVolUnitBall = PI^(hdim) / Factorial(hdim)
57 3 : logVolUnitBall = PI
58 4 : do idim = 2_IK, ndim / 2_IK
59 4 : logVolUnitBall = logVolUnitBall * PI / idim
60 : end do
61 : else
62 : ! ndim is an odd integer
63 : ! ndim = 2 * hdim - 1
64 : ! gamma(ndim / 2 + 1) = gamma(hdim + 1 / 2);
65 : ! gamma(hdim + 1 / 2) = sqrt(PI) * (2 * hdim)! / (4**hdim * hdim!)
66 2 : hdim = ndim + 1_IK
67 2 : logVolUnitBall = 4._RKC / (hdim + 1_IK) ! This is to avoid an extra unnecessary division of `logVolUnitBall` by `PI`.
68 6 : do idim = hdim + 2_IK, 2_IK * hdim
69 : ! logVolUnitBall
70 : ! = PI**(hdim - 1 / 2) / gamma(hdim + 1 / 2)
71 : ! = PI**(hdim + 1) * 4**hdim * hdim! / (2 * hdim)!
72 6 : logVolUnitBall = logVolUnitBall * FOUR_PI / idim
73 : end do
74 :
75 : end if
76 5 : logVolUnitBall = log(logVolUnitBall)
77 : else
78 1 : logVolUnitBall = 0._RKC
79 : end if
80 :
81 : !%%%%%%%%%%%%%%%%%%%
82 : #elif getLogVolEll_ENABLED
83 : !%%%%%%%%%%%%%%%%%%%
84 :
85 3 : CHECK_ASSERTION(__LINE__, size(gramian, 1, IK) == size(gramian, 2, IK), SK_"@setLogVolEllipsoid(): The condition `size(gramian, 1) == size(gramian, 2)` must hold. shape(gramian) = "//getStr(shape(gramian, IK)))
86 3 : logVolEll = getLogVolUnitBall(real(size(gramian, 1, IK), RKC)) + getMatDetSqrtLog(gramian)
87 :
88 : !%%%%%%%%%%%%%%%%%%%%%%%%
89 : #elif getCountMemberEll_ENABLED
90 : !%%%%%%%%%%%%%%%%%%%%%%%%
91 :
92 : ! Define the center argument.
93 : #if Org_ENABLED
94 : #define CENTER
95 : #elif Cen_ENABLED
96 : #define CENTER center,
97 : #else
98 : #error "Unrecognized interface."
99 : #endif
100 : integer(IK) :: ndim, ipnt
101 : countMemberEll = 0_IK
102 4 : ndim = size(point, 1, IK)
103 20 : do ipnt = 1_IK, size(point, 2, IK)
104 : if ( & ! LCOV_EXCL_LINE
105 : #if Sph_ENABLED
106 : isMemberEll(radius, CENTER point(1:ndim, ipnt)) & ! LCOV_EXCL_LINE
107 : #elif Ell_ENABLED
108 : isMemberEll(invGram, CENTER point(1:ndim, ipnt)) & ! LCOV_EXCL_LINE
109 : #else
110 : #error "Unrecognized interface."
111 : #endif
112 10 : ) countMemberEll = countMemberEll + 1_IK
113 : end do
114 :
115 : !%%%%%%%%%%%%%%%%%%
116 : #elif isMemberEll_ENABLED
117 : !%%%%%%%%%%%%%%%%%%
118 :
119 : integer(IK) :: ndim
120 : #if Cen_ENABLED
121 24 : real(RKC) :: normedPoint(size(point, 1, IK))
122 36 : CHECK_ASSERTION(__LINE__, size(point, 1, IK) == size(center, 1, IK), SK_"@isMemberEll(): The condition `size(point) == size(center)` must hold. size(point), size(center) = "//getStr([size(point, 1, IK), size(center, 1, IK)]))
123 : #elif !Org_ENABLED
124 : #error "Unrecognized interface."
125 : #endif
126 : #if Sph_ENABLED
127 14 : CHECK_ASSERTION(__LINE__, 0._RKC < radius, SK_"@isMemberEll(): The condition `0. < radius` must hold. radius = "//getStr(radius))
128 : #elif Ell_ENABLED
129 96 : CHECK_ASSERTION(__LINE__, all(size(point, 1, IK) == shape(invGram, IK)), SK_"@isMemberEll(): The condition `all(size(point, 1) == shape(invGram))` must hold. size(point, 1), shape(invGram) = "//getStr([size(point, 1, IK), shape(invGram, IK)]))
130 : #endif
131 8 : ndim = size(point, 1, IK)
132 : #if Org_ENABLED
133 : #define NORMEDPOINT point(1 : ndim)
134 : #elif Cen_ENABLED
135 : #define NORMEDPOINT normedPoint
136 36 : normedPoint(1 : ndim) = point(1 : ndim) - center
137 : #else
138 : #error "Unrecognized interface."
139 : #endif
140 : #if Sph_ENABLED
141 42 : isMember = logical(dot_product(NORMEDPOINT, NORMEDPOINT) <= radius, LK)
142 : #elif Ell_ENABLED
143 150 : isMember = logical(dot_product(NORMEDPOINT, matmul(invGram, NORMEDPOINT)) <= 1._RKC, LK)
144 : #else
145 : #error "Unrecognized interface."
146 : #endif
147 :
148 : #else
149 : !%%%%%%%%%%%%%%%%%%%%%%%%
150 : #error "Unrecognized interface."
151 : !%%%%%%%%%%%%%%%%%%%%%%%%
152 : #endif
153 :
154 : #undef NORMEDPOINT
155 : #undef CENTER
|
