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_distCov](@ref pm_distCov).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Monday March 6, 2017, 3:22 pm, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin.<br>
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if getCovRand_ENABLED && GRAM_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 : type(rngf_type) :: rng
32 : #if S0_ENABLED
33 1434 : if (present(scale)) then
34 10 : call setCovRand(rng, rand, scale)
35 : else
36 1424 : call setCovRand(rng, rand)
37 : end if
38 : #elif S1_ENABLED
39 31 : call setCovRand(rng, rand, scale)
40 : #else
41 : #error "Unrecognized interface."
42 : #endif
43 :
44 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45 : #elif getCovRand_ENABLED && (DVINE_ENABLED || ONION_ENABLED)
46 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47 :
48 : type(rngf_type) :: rng
49 : ! The onion method appears unstable, use the slower but more stable method for now.
50 : !type(onion_type) :: onion
51 : #if S0_ENABLED
52 : if (present(scale)) then
53 : call setCovRand(rng, rand, method, eta, scale)
54 : else
55 : call setCovRand(rng, rand, method, eta)
56 : end if
57 : #elif S1_ENABLED
58 : call setCovRand(rng, rand, method, eta, scale)
59 : #else
60 : #error "Unrecognized interface."
61 : #endif
62 : #if ONION_ENABLED
63 : if (onion%info /= 0_IK) error stop "@getCovRand(): The Cholesky factorization of the onion method failed."
64 : #endif
65 : !call setMatCopy(rand(2:, 1:), rdpack, rand(1:, 2:), rdpack, uppDia)
66 :
67 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
68 : #elif getCovRand_ENABLED && GRAM_ENABLED
69 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70 :
71 : #error "Unrecognized interface."
72 : ! type(rngf_type) :: rng
73 : !#if S0_ENABLED
74 : ! if (present(scale)) then
75 : ! call setCovRand(rng, rand, scale)
76 : ! else
77 : ! call setCovRand(rng, rand)
78 : ! end if
79 : !#elif S1_ENABLED
80 : ! call setCovRand(rng, rand, scale)
81 : !#else
82 : !#error "Unrecognized interface."
83 : !#endif
84 :
85 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86 : #elif setCovRand_ENABLED && GRAM_ENABLED
87 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88 :
89 : real(TKC) :: normfac
90 : real(TKC), parameter :: EPS = 2 * sqrt(epsilon(0._TKC))
91 : ! Define the output kind.
92 : #if CK_ENABLED
93 : #define TYPE_OF_RAND complex(TKC)
94 : #define GET_CONJG(X)conjg(X)
95 : #define GET_RE(X)X%re
96 : complex(TKC), parameter :: LB = -cmplx(1._TKC, 1._TKC, TKC), UB = cmplx(1._TKC, 1._TKC, TKC), ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 0._TKC)
97 : #elif RK_ENABLED
98 : #define GET_RE(X)X
99 : #define GET_CONJG(X)X
100 : #define TYPE_OF_RAND real(TKC)
101 : real(TKC), parameter :: LB = -1._TKC, UB = 1._TKC, ZERO = 0._TKC, ONE = 1._TKC
102 : #else
103 : #error "Unrecognized interface."
104 : #endif
105 1536 : TYPE_OF_RAND :: upper(size(rand, 1, IK), size(rand, 2, IK))
106 : integer(IK) :: idim, jdim, ndim
107 : #if SD_ENABLED
108 : #define GET_SCALED(X)X
109 : #elif S0_ENABLED
110 : #define GET_SCALED(X)X * scaledSq
111 : real(TKC) :: scaledSq
112 20 : scaledSq = scale**2
113 20 : CHECK_ASSERTION(__LINE__, 0._TKC < scale, SK_"@setCovRand(): The condition `0. < scale` must hold. scale = "//getStr(scale))
114 : #elif S1_ENABLED
115 : #define GET_SCALED(X)X * scale(idim) * scale(jdim)
116 208 : CHECK_ASSERTION(__LINE__, all([0._TKC < scale]), SK_"@setCovRand(): The condition `all([0. < scale])` must hold. scale = "//getStr(scale))
117 246 : CHECK_ASSERTION(__LINE__, size(scale, 1, IK) == size(rand, 1, IK), SK_"@setCovRand(): The condition `size(scale) == size(rand, 1, IK)` must hold. size(scale), shape(rand) = "//getStr([size(scale, 1, IK), shape(rand, IK)]))
118 : #else
119 : #error "Unrecognized interface."
120 : #endif
121 1495 : CHECK_ASSERTION(__LINE__, size(rand, 1, IK) == size(rand, 2, IK), SK_"@setCovRand(): The condition `size(rand, 1) == size(rand, 2)` must hold. shape(rand) = "//getStr(shape(rand, IK)))
122 : ndim = size(rand, 1, IK)
123 1495 : if (ndim < 1_IK) return
124 7324 : do idim = 1, ndim
125 : do
126 34502 : call setUnifRand(rng, upper(1 : ndim, idim), LB, UB)
127 34502 : normfac = real(sqrt(dot_product(upper(1 : ndim, idim), upper(1 : ndim, idim))), TKC) ! \todo: The performance of this expression can be improved by replacing `dot_product` with `absq()`.
128 5833 : if (ZERO == normfac) cycle
129 0 : exit
130 : end do
131 4141 : normfac = 1._TKC / normfac
132 35993 : upper(1 : ndim, idim) = upper(1 : ndim, idim) * normfac
133 : end do
134 229352 : rand = matmul(transpose(GET_CONJG(upper)), upper)
135 : ! Ensure the diagonals are all pure 1, free from numerical round-off errors.
136 7324 : do jdim = 1, ndim
137 7324 : rand(jdim, jdim) = ONE
138 : end do
139 : !block
140 : ! use pm_io
141 : ! call disp%show("rand")
142 : ! call disp%show( rand )
143 : !end block
144 : !call setCor(rand, uppDia, rand, uppDia)
145 : ! Rescale.
146 : #if S0_ENABLED || S1_ENABLED
147 310 : do jdim = 1, ndim
148 1117 : do idim = 1, jdim
149 1060 : rand(idim, jdim) = GET_SCALED(rand(idim, jdim))
150 : end do
151 : end do
152 : #endif
153 : !block
154 : !use pm_io, only: disp
155 : !call disp%show("rand")
156 : !call disp%show( rand )
157 : !end block
158 : ! copy to the lower triangle.
159 1491 : call setMatCopy(rand, rdpack, rand, rdpack, upp, transHerm)
160 : #undef TYPE_OF_RAND
161 : #undef GET_SCALED
162 : #undef GET_CONJG
163 : #undef GET_RE
164 :
165 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
166 : #elif setCovRand_ENABLED && (DVINE_ENABLED || ONION_ENABLED)
167 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
168 :
169 : #define GET_INDEX(i,j) i, j
170 : !#if UXD_ENABLED
171 : !#define GET_INDEX(i,j) j, i
172 : !#elif XLD_ENABLED
173 : !#define GET_INDEX(i,j) i, j
174 : !#else
175 : !#error "Unrecognized interface."
176 : !#endif
177 : real(TKC) :: beta
178 : integer(IK) :: ndim, idim, jdim, kdim
179 : ! Set the scale
180 : #if SD_ENABLED
181 : real(TKC), parameter :: scaleSq = 1._TKC
182 : #define GET_SCALED(x, i, j)
183 : #else
184 : #if S0_ENABLED
185 : real(TKC) :: scaleSq
186 20 : scaleSq = scale * scale
187 : #define GET_SCALED(x, i, j) x = x * scaleSq
188 : #elif S1_ENABLED
189 12 : CHECK_ASSERTION(__LINE__, size(scale, 1, IK) == size(rand, 1, IK), SK_"@setCovRand(): The condition `size(scale) == size(rand, 1, IK)` must hold. size(scale), shape(rand) = "//getStr([size(scale, 1, IK), shape(rand, IK)]))
190 : #define GET_SCALED(x, i, j) x = x * scale(i) * scale(j)
191 : #else
192 : #error "Unrecognized interface."
193 : #endif
194 38 : CHECK_ASSERTION(__LINE__, all([0._TKC < scale]), SK_"@setCovRand(): The condition `all([0. < scale])` must hold. scale = "//getStr(scale))
195 : #endif
196 42 : CHECK_ASSERTION(__LINE__, 0._TKC <= eta, SK_"@setCovRand(): The condition `0. <= eta` must hold. eta = "//getStr(eta))
197 : !CHECK_ASSERTION(__LINE__, 0_IK < size(rand, 1, IK), SK_"@setCovRand(): The condition `0 < size(rand, 1)` must hold. shape(rand) = "//getStr(shape(rand)))
198 42 : CHECK_ASSERTION(__LINE__, size(rand, 1, IK) == size(rand, 2, IK), SK_"@setCovRand(): The condition `size(rand, 1) <= size(rand, 2)` must hold. shape(rand) = "//getStr(shape(rand)))
199 : ndim = size(rand, 1, IK)
200 42 : if (1_IK < ndim) then
201 : #if ONION_ENABLED
202 21 : beta = eta + 0.5_TKC * ndim
203 : block
204 : integer(IK) :: kdimp1
205 42 : real(TKC) :: brand, ssnrandInv, chol(ndim - 1, ndim - 1), wrand(ndim - 1) ! sphere rv.
206 0 : do
207 21 : call setBetaRand(rng, brand, beta, beta)
208 21 : if (0._TKC < brand .and. brand < 1._TKC) exit
209 : end do
210 21 : rand(1, 1) = 1._TKC
211 21 : rand(2, 2) = 1._TKC
212 21 : rand(1, 2) = 2._TKC * brand - 1._TKC
213 57 : do kdim = 2_IK, ndim - 1_IK
214 41 : call setNormRand(rng, wrand(1 : kdim))
215 150 : ssnrandInv = 1._TKC / norm2(wrand(1 : kdim)) ! `wrand(1 : kdim) * ssnrandInv` would be a uniform random point on the surface of the kdim-sphere.
216 41 : beta = beta - 0.5_TKC
217 : do
218 41 : call setBetaRand(rng, brand, 0.5_TKC * kdim, beta)
219 41 : if (0._TKC < brand .and. brand < 1._TKC) exit
220 : end do
221 41 : ssnrandInv = ssnrandInv * sqrt(brand)
222 41 : call setMatChol(rand, uppDia, method%info, chol, nothing, kdim, 0_IK, 0_IK, 0_IK, 0_IK)
223 41 : if (method%info /= 0_IK) return
224 : ! Compute z = cholow * wrand.
225 36 : kdimp1 = kdim + 1
226 128 : rand(1 : kdim, kdimp1) = 0._TKC
227 128 : do jdim = 1, kdim
228 92 : rand(jdim, kdimp1) = rand(jdim, kdimp1) + chol(jdim, jdim) * wrand(jdim) * ssnrandInv
229 207 : do idim = jdim + 1, kdim
230 171 : rand(idim, kdimp1) = rand(idim, kdimp1) + rand(idim, jdim) * rand(jdim, kdimp1)
231 : end do
232 : end do
233 52 : rand(kdimp1, kdimp1) = 1._TKC
234 : !block
235 : !use pm_io, only: display_type
236 : !type(display_type) :: disp
237 : !call disp%skip
238 : !call disp%show(rand(1 : kdimp1, 1 : kdimp1))
239 : !call disp%skip
240 : !end block
241 : end do
242 77 : do kdim = ndim, 1, -1
243 : #if SD_ENABLED || S0_ENABLED
244 61 : rand(kdim, kdim) = scaleSq
245 : #elif S1_ENABLED
246 0 : rand(kdim, kdim) = rand(kdim, kdim) * scale(kdim) * scale(kdim)
247 : #endif
248 169 : do idim = kdim + 1, ndim
249 : #if SD_ENABLED || S0_ENABLED
250 92 : rand(kdim, idim) = rand(kdim, idim) * scaleSq
251 : #elif S1_ENABLED
252 0 : rand(kdim, idim) = rand(kdim, idim) * scale(idim) * scale(kdim)
253 : #endif
254 153 : rand(idim, kdim) = rand(kdim, idim)
255 : end do
256 : end do
257 : !call setMatCopy(rand(2 : ndim, 1 : ndim), rdpack, rand(1 : ndim, 2 : ndim), rdpack, uppDia)
258 : end block
259 : !block
260 : !use pm_io, only: display_type
261 : !type(display_type) :: disp
262 : !call disp%skip
263 : !call disp%show(rand)
264 : !call disp%skip
265 : !end block
266 : #elif DVINE_ENABLED
267 21 : beta = eta + 0.5_TKC * (ndim + 1_IK)
268 : block
269 21 : real(TKC) :: parCorMat(size(rand, 1, IK), size(rand, 1, IK))
270 509 : parCorMat = 0._TKC
271 88 : do kdim = 1_IK, ndim - 1_IK
272 67 : beta = beta - 0.5_TKC
273 : #if SD_ENABLED || S0_ENABLED
274 60 : rand(kdim, kdim) = scaleSq
275 : #elif S1_ENABLED
276 7 : rand(kdim, kdim) = scale(kdim) * scale(kdim)
277 : #endif
278 244 : do idim = kdim + 1_IK, ndim
279 67 : loopSensibleBeta: do
280 156 : call setBetaRand(rng, parCorMat(kdim, idim), beta, beta)
281 156 : if (0._TKC < parCorMat(kdim, idim) .and. parCorMat(kdim, idim) < 1._TKC) then
282 156 : parCorMat(kdim, idim) = 2 * parCorMat(kdim, idim) - 1._TKC ! Linearly shift to the range (-1, 1).
283 : ! Convert the partial correlations to full correlation.
284 156 : rand(GET_INDEX(kdim, idim)) = parCorMat(kdim, idim)
285 314 : do jdim = kdim - 1_IK, 1_IK, -1_IK
286 : rand(GET_INDEX(kdim, idim)) = parCorMat(jdim, idim) * parCorMat(jdim, kdim) + & ! LCOV_EXCL_LINE
287 314 : rand(GET_INDEX(kdim, idim)) * sqrt((1._TKC - parCorMat(jdim, idim)**2) * (1._TKC - parCorMat(jdim, kdim)**2))
288 : end do
289 92 : GET_SCALED(rand(GET_INDEX(kdim, idim)), idim, kdim)
290 156 : rand(GET_INDEX(idim, kdim)) = rand(GET_INDEX(kdim, idim))
291 : exit loopSensibleBeta
292 : end if
293 : end do loopSensibleBeta
294 : end do
295 : end do
296 : #if SD_ENABLED || S0_ENABLED
297 20 : rand(kdim, kdim) = scaleSq
298 : #elif S1_ENABLED
299 1 : rand(kdim, kdim) = scale(kdim) * scale(kdim)
300 : #endif
301 : end block
302 : #else
303 : #error "Unrecognized interface."
304 : #endif
305 0 : elseif (1_IK == ndim) then
306 : #if SD_ENABLED || S0_ENABLED
307 0 : rand(1,1) = scaleSq
308 : #elif S1_ENABLED
309 0 : rand(1,1) = scale(1) * scale(1)
310 : #endif
311 : end if
312 : #if ONION_ENABLED
313 16 : method%info = 0_IK
314 : #endif
315 :
316 : #else
317 : !%%%%%%%%%%%%%%%%%%%%%%%%
318 : #error "Unrecognized interface."
319 : !%%%%%%%%%%%%%%%%%%%%%%%%
320 : #endif
321 : #undef GET_SCALED
322 : #undef GET_INDEX
|
