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_distanceMahal](@ref pm_distanceMahal).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Check the positive-definiteness of `invCov`.
28 : #define CHECK_POSDEF(INVCOV,ISAM) \
29 : CHECK_ASSERTION(__LINE__, isMatClass(INVCOV, posdefmat), \
30 : SK_"@setMahalSq(): The condition `isMatClass(invCov(:,:,isam), posdefmat)` must hold. isam, invCov(:, :, isam) = "//getStr(ISAM)//SK_", "//getStr(INVCOV))
31 : ! Check bounds of `invCov`.
32 : #define CHECK_LEN_INVCOV \
33 : CHECK_ASSERTION(__LINE__, all(size(point, 1, IK) == [size(invCov, 1, IK), size(invCov, 2, IK)]), \
34 : SK_"@setMahalSq(): The condition `all(size(point, 1) == [size(invCov, 1), size(invCov, 2)])` must hold. size(point), shape(invCov) = "//\
35 : getStr([size(point, 1, IK), shape(invCov, IK)]))
36 : ! Check bounds of `invCov`.
37 : #if InvDef_ENABLED
38 : #define CHECK_LEN_CENTER
39 : #define CHECK_LEN_CENTER_INVCOV
40 : #elif InvCen_ENABLED
41 : #define CHECK_LEN_CENTER \
42 : CHECK_ASSERTION(__LINE__, size(point, 1) == size(center, 1, IK), \
43 : SK_"@setMahalSq(): The condition `size(point, 1) == size(center, 1)` must hold. size(point, 1), size(center) = "//\
44 : getStr([size(point, 1, IK), size(center, 1, IK)]))
45 : ! Check bounds of `invCov`.
46 : #define CHECK_LEN_CENTER_INVCOV \
47 : CHECK_ASSERTION(__LINE__, size(center, rank(center), IK) == size(invCov, rank(invCov), IK), \
48 : SK_"@setMahalSq(): The condition `size(center, rank(center)) == size(invCov, rank(invCov))` must hold. shape(center), shape(invCov) = "//\
49 : getStr([shape(center, IK), shape(invCov, IK)]))
50 : #else
51 : #error "Unrecognized interface."
52 : #endif
53 : ! \bug
54 : ! ifort 2021.8
55 : ! error #6401: The attributes of this name conflict with those made accessible by a USE statement. [SIZE]
56 : ! It appears ifort cannot digest `size` intrinsic in the `do concurrent` declaration.
57 : !intrinsic :: size
58 : ! Set the type and kind.
59 : #if RK_ENABLED
60 : #define GET_CONJG(point) point
61 : #define TYPE_KIND real(RKC)
62 : real(RKC) , parameter :: ZERO = 0._RKC
63 : #elif CK_ENABLED
64 : #define GET_CONJG(point) conjg(point)
65 : #define TYPE_KIND complex(CKC)
66 : complex(CKC), parameter :: ZERO = (0._CKC, 0._CKC)
67 : #else
68 : #error "Unrecognized interface."
69 : #endif
70 : #if RK_ENABLED && !D0_ENABLED
71 : integer(IK) :: idim
72 : #elif !(CK_ENABLED || D0_ENABLED)
73 : #error "Unrecognized interface."
74 : #endif
75 : !%%%%%%%%%
76 : #if D0_ENABLED
77 : !%%%%%%%%%
78 :
79 : #if RK_ENABLED
80 24 : CHECK_ASSERTION(__LINE__, 0._RKC < invCov, SK_": The condition `0. < invCov` must hold. invCov = "//getStr(invCov))
81 : #elif !CK_ENABLED
82 : #error "Unrecognized interface."
83 : #endif
84 : ! Compute the distance.
85 : #if InvDef_ENABLED
86 24 : mahalSq = GET_CONJG(point) * invCov * point
87 : #elif InvCen_ENABLED
88 24 : mahalSq = GET_CONJG((point - center)) * invCov * (point - center)
89 : #endif
90 :
91 : !%%%%%%%%%%%%%%%%%%%%%%%%
92 : #elif One_ENABLED && D1_ENABLED
93 : !%%%%%%%%%%%%%%%%%%%%%%%%
94 :
95 : integer(IK) :: ndim
96 : #if InvDef_ENABLED
97 : #define PNORMED(IDIM) point(IDIM)
98 : #elif InvCen_ENABLED
99 250008 : TYPE_KIND :: pnormed(size(point, 1, IK))
100 750026 : pnormed = point - center
101 : #endif
102 4 : ndim = size(point, 1, IK)
103 :
104 250012 : CHECK_POSDEF(invCov, 1)
105 750024 : CHECK_LEN_CENTER
106 1500072 : CHECK_LEN_INVCOV
107 :
108 : ! Compute the distances.
109 : #if CK_ENABLED
110 82 : mahalSq = dot_product(PNORMED(1:ndim), matmul(invCov, PNORMED(1:ndim))) ! fpp
111 : #elif RK_ENABLED
112 : !mahalSq = dot_product(PNORMED , matmul(PNORMED, invCov)) ! fpp
113 2 : mahalSq = ZERO
114 750026 : do idim = 1_IK, ndim
115 1750072 : mahalSq = mahalSq + PNORMED(idim) * dot_product(PNORMED(1:ndim), invCov(1:ndim, idim))
116 : end do
117 : #endif
118 :
119 : !%%%%%%%%%%%%%%%%%%%%%%%%
120 : #elif One_ENABLED && D2_ENABLED
121 : !%%%%%%%%%%%%%%%%%%%%%%%%
122 :
123 : integer(IK) :: ipnt, ndim
124 : #if InvDef_ENABLED
125 : #define PNORMED(IDIM) point(IDIM, ipnt)
126 : #elif InvCen_ENABLED
127 8 : TYPE_KIND :: pnormed(size(point, 1, IK))
128 : #endif
129 4 : ndim = size(point, 1, IK)
130 :
131 8 : CHECK_POSDEF(invCov, 1)
132 12 : CHECK_LEN_CENTER
133 48 : CHECK_LEN_INVCOV
134 : #if setMahalSq_ENABLED
135 12 : CHECK_ASSERTION(__LINE__, size(point, 2, IK) == size(mahalSq, 1, IK), \
136 : SK_"@setMahalSq(): The condition `size(point, 2) == size(mahalSq, 1)` must hold. size(point, 2), shape(mahalSq, 1) = "//\
137 : getStr([size(point, 2, IK), size(mahalSq, 1, IK)]))
138 : #endif
139 :
140 : ! Compute the distances.
141 :
142 48 : do ipnt = 1_IK, size(point, 2, IK)
143 : #if InvCen_ENABLED
144 80 : pnormed(1:ndim) = point(1:ndim, ipnt) - center
145 : #endif
146 : #if CK_ENABLED
147 414 : mahalSq(ipnt) = dot_product(PNORMED(1:ndim), matmul(invCov, PNORMED(1:ndim))) ! fpp
148 : #elif RK_ENABLED
149 : !mahalSq(ipnt) = getMahalSq(PNORMED(1:ndim), invCov)
150 20 : mahalSq(ipnt) = ZERO
151 84 : do idim = 1_IK, ndim
152 260 : mahalSq(ipnt) = mahalSq(ipnt) + PNORMED(idim) * dot_product(PNORMED(1:ndim), invCov(1:ndim, idim))
153 : end do
154 : #endif
155 : end do
156 :
157 : !%%%%%%%%%%%%%%%%%%%%%%%%
158 : #elif Mix_ENABLED && D1_ENABLED
159 : !%%%%%%%%%%%%%%%%%%%%%%%%
160 :
161 : integer(IK) :: isam, ndim
162 : #if InvDef_ENABLED
163 : #define PNORMED(IDIM) point(IDIM)
164 : #elif InvCen_ENABLED
165 982008 : TYPE_KIND :: pnormed(size(point, 1, IK))
166 : #endif
167 4 : ndim = size(point, 1, IK)
168 :
169 5401044 : CHECK_LEN_CENTER_INVCOV
170 1473012 : CHECK_LEN_CENTER
171 3928064 : CHECK_LEN_INVCOV
172 : #if setMahalSq_ENABLED
173 12 : CHECK_ASSERTION(__LINE__, size(invCov, 3, IK) == size(mahalSq, 1, IK), \
174 : SK_"@setMahalSq(): The condition `size(invCov, 3) == size(mahalSq, 1)` must hold. size(invCov, 3), shape(mahalSq, 1) = "//\
175 : getStr([size(invCov, 3, IK), size(mahalSq, 1, IK)]))
176 : #endif
177 :
178 : ! Compute the distances.
179 :
180 1475024 : do isam = 1_IK, size(invCov, 3, IK)
181 984016 : CHECK_POSDEF(invCov(:, :, isam), isam)
182 : #if InvCen_ENABLED
183 2948032 : pnormed(1:ndim) = point(1:ndim) - center(1:ndim, isam)
184 : #endif
185 : #if CK_ENABLED
186 168 : mahalSq(isam) = dot_product(PNORMED(1:ndim), matmul(invCov(:, :, isam), PNORMED(1:ndim))) ! fpp
187 : #elif RK_ENABLED
188 984008 : mahalSq(isam) = ZERO
189 3439036 : do idim = 1_IK, ndim
190 6872104 : mahalSq(isam) = mahalSq(isam) + PNORMED(idim) * dot_product(PNORMED(1:ndim), invCov(1:ndim, idim, isam))
191 : end do
192 : #endif
193 : end do
194 :
195 : !%%%%%%%%%%%%%%%%%%%%%%%%
196 : #elif Mix_ENABLED && D2_ENABLED
197 : !%%%%%%%%%%%%%%%%%%%%%%%%
198 :
199 : integer(IK) :: ipnt, isam, nsam, ndim
200 : #if InvDef_ENABLED
201 : #define PNORMED(IDIM) point(IDIM, ipnt)
202 : #elif InvCen_ENABLED
203 8 : TYPE_KIND :: pnormed(size(point, 1, IK))
204 : #endif
205 4 : ndim = size(point, 1, IK)
206 4 : nsam = size(invCov, 3, IK)
207 :
208 44 : CHECK_LEN_CENTER_INVCOV
209 12 : CHECK_LEN_CENTER
210 64 : CHECK_LEN_INVCOV
211 : #if setMahalSq_ENABLED
212 36 : CHECK_ASSERTION(__LINE__, all([size(invCov, 3, IK), size(point, 2, IK)] == shape(mahalSq, IK)), \
213 : SK_"@setMahalSq(): The condition `all([size(invCov, 3), size(point, 2)] == shape(mahalSq))` must hold. size(invCov, 3), size(point, 2), shape(mahalSq) = "//\
214 : getStr([size(invCov, 3, IK), size(point, 2, IK), shape(mahalSq, IK)]))
215 : #endif
216 :
217 : ! Compute the distances.
218 :
219 48 : do ipnt = 1_IK, size(point, 2, IK)
220 128 : do isam = 1_IK, nsam
221 : #if InvCen_ENABLED
222 160 : pnormed(1:ndim) = point(1:ndim, ipnt) - center(1:ndim, isam)
223 : #endif
224 : #if CK_ENABLED
225 840 : mahalSq(isam, ipnt) = dot_product(PNORMED(1:ndim), matmul(invCov(:, :, isam), PNORMED(1:ndim))) ! fpp
226 : #elif RK_ENABLED
227 40 : mahalSq(isam, ipnt) = ZERO
228 180 : do idim = 1_IK, ndim
229 520 : mahalSq(isam, ipnt) = mahalSq(isam, ipnt) + PNORMED(idim) * dot_product(PNORMED(1:ndim), invCov(1:ndim, idim, isam))
230 : end do
231 : #endif
232 : end do
233 : end do
234 :
235 : #else
236 : !%%%%%%%%%%%%%%%%%%%%%%%%
237 : #error "Unrecognized interface."
238 : !%%%%%%%%%%%%%%%%%%%%%%%%
239 : #endif
240 : #undef CHECK_LEN_CENTER_INVCOV
241 : #undef CHECK_LEN_CENTER
242 : #undef CHECK_LEN_INVCOV
243 : #undef CHECK_POSDEF
244 : #undef TYPE_KIND
245 : #undef GET_CONJG
246 : #undef PNORMED
|
