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 procedure implementation of [pm_distMultiNorm](@ref pm_distMultiNorm).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
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 getMNR_ENABLED || setMNR_ENABLED
33 : #error "Unrecognized interface."
34 : #endif
35 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 : #if getMultiNormLogPDFNF_ENABLED
37 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
38 :
39 : real(RKC), parameter :: LOG_INVERSE_SQRT_TWO_PI = -0.5_RKC * log(2 * acos(-1._RKC))
40 : #if I_ENABLED
41 : integer(IK) :: ndim
42 : real(RKC) :: logSqrtDetInvCov
43 250003 : logSqrtDetInvCov = getMatDetSqrtLog(invCov, uppDia)
44 250003 : ndim = size(invCov, 1, IK)
45 : #elif IF_ENABLED
46 : integer(IK) :: ndim
47 2 : real(RKC) :: logSqrtDetInvCov, chol(size(invCov, 1, IK), size(invCov, 1, IK))
48 1 : CHECK_ASSERTION(__LINE__, size(invCov, 1, IK) == size(invCov, 2, IK), SK_"@getMultiNormLogPDFNF(): The condition `size(invCov, 1) == size(invCov, 2)` must hold. shape(invCov) = "//getStr(shape(invCov, IK)))
49 1 : call setMatDetSqrtLog(invCov, uppDia, logSqrtDetInvCov, info, chol, nothing)
50 : ndim = size(invCov, 1, IK)
51 1 : if (info /= 0_IK) return
52 : #elif ND_ENABLED
53 613108 : CHECK_ASSERTION(__LINE__, 0_IK < ndim, SK_"@getMultiNormLogPDFNF(): The input `ndim` must be positive. ndim = "//getStr(ndim))
54 : #else
55 : #error "Unrecognized interface."
56 : #endif
57 1 : logPDFNF = ndim * LOG_INVERSE_SQRT_TWO_PI + logSqrtDetInvCov
58 :
59 : !%%%%%%%%%%%%%%%%%%%%%%%%%
60 : #elif getMultiNormLogPDF_ENABLED
61 : !%%%%%%%%%%%%%%%%%%%%%%%%%
62 :
63 : ! Set the normalization factor for different interfaces.
64 : #if DDD_ENABLED || MDD_ENABLED
65 : #define LOGNORMFAC getMultiNormLogPDFNF(size(X,1,IK), 0._RKC)
66 : #elif DID_ENABLED || MID_ENABLED
67 : #define LOGNORMFAC getMultiNormLogPDFNF(invCov)
68 : #elif DDN_ENABLED || MDN_ENABLED || DIN_ENABLED || MIN_ENABLED
69 : #define LOGNORMFAC logPDFNF
70 : #else
71 : #error "Unrecognized interface."
72 : #endif
73 : ! Set the Mahalanobis distance.
74 : #if DDD_ENABLED || DDN_ENABLED
75 : #define MAHAL_SQ sum(X**2)
76 : #elif DID_ENABLED || DIN_ENABLED
77 : #define MAHAL_SQ getMahalSq(X, invCov)
78 0 : CHECK_ASSERTION(__LINE__, all(shape(invCov, IK) == size(X, 1, IK)), SK_"@getMultiNormLogPDF(): The condition `all(shape(invCov) == size(X))` must hold. shape(invCov), size(X) = "//getStr([shape(invCov, IK), size(X, 1, IK)]))
79 : #elif MID_ENABLED || MIN_ENABLED
80 : #define MAHAL_SQ getMahalSq(X, invCov, mean)
81 2000016 : CHECK_ASSERTION(__LINE__, all(shape(invCov, IK) == size(X, 1, IK)), SK_"@getMultiNormLogPDF(): The condition `all(shape(invCov) == size(X))` must hold. shape(invCov), size(X) = "//getStr([shape(invCov, IK), size(X, 1, IK)]))
82 750006 : CHECK_ASSERTION(__LINE__, size(mean, 1, IK) == size(X, 1, IK), SK_"@getMultiNormLogPDF(): The condition `size(mean) == size(X)` must hold. size(mean), size(X) = "//getStr([size(mean, 1, IK), size(X, 1, IK)]))
83 : #elif MDD_ENABLED || MDN_ENABLED
84 : #define MAHAL_SQ sum((X - mean)**2)
85 1839312 : CHECK_ASSERTION(__LINE__, size(mean, 1, IK) == size(X, 1, IK), SK_"@getMultiNormLogPDF(): The condition `size(mean) == size(X)` must hold. size(mean), size(X) = "//getStr([size(mean, 1, IK), size(X, 1, IK)]))
86 : #else
87 : #error "Unrecognized interface."
88 : #endif
89 3315524 : logPDF = LOGNORMFAC - 0.5_RKC * MAHAL_SQ
90 :
91 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
92 : #elif D1_ENABLED && (getMNR_ENABLED || setMNR_ENABLED)
93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94 :
95 : integer(IK) :: ndim
96 1121986 : ndim = size(rand, kind = IK)
97 : #if AM_ENABLED
98 3305952 : CHECK_ASSERTION(__LINE__, size(rand, 1, IK) == size(mean, 1, IK), SK_"@setMultiNormRand(): The condition `size(rand) == size(mean)` must hold. size(rand), size(mean) = "//getStr([size(rand, 1, IK), size(mean, 1, IK)]))
99 : #elif !DM_ENABLED
100 : #error "Unrecognized interface."
101 : #endif
102 : ! Generate random MVN.
103 : #if DC_ENABLED && DM_ENABLED
104 5001 : call setNormRand(RNG rand)
105 : #elif DC_ENABLED && AM_ENABLED
106 10002 : call setNormRand(RNG rand)
107 30006 : rand = rand + mean
108 : #elif AC_ENABLED && (AM_ENABLED || DM_ENABLED)
109 8895872 : CHECK_ASSERTION(__LINE__, all(size(rand, 1, IK) == shape(chol, IK)), SK_"@setMultiNormRand(): The condition `all(size(rand) == shape(chol))` must hold. size(rand), shape(chol) = "//getStr([size(rand, 1, IK), shape(chol, IK)]))
110 : ! Define the indexing rules.
111 : #if XLD_ENABLED
112 : #define GET_INDEX(I,J)I,J
113 : #elif UXD_ENABLED
114 : #define GET_INDEX(I,J)J,I
115 : #elif D1_ENABLED && setMUR_ENABLED
116 : #error "Unrecognized interface."
117 : #endif
118 : ! Define the default mean.
119 : #if DM_ENABLED
120 : #define XPLUS(X)
121 : #elif AM_ENABLED
122 : #define XPLUS(X)X +
123 : #else
124 : #error "Unrecognized interface."
125 : #endif
126 : block
127 : integer(IK) :: idim
128 : real(RKC) :: normrnd
129 : ! Separate the first to allow the possibility of adding an optional `mean`.
130 1111984 : call setNormRand(RNG normrnd)
131 4876224 : rand(1 : ndim) = XPLUS(mean) chol(GET_INDEX(1 : ndim, 1)) * normrnd
132 3764240 : do idim = 2_IK, ndim
133 20004 : call setNormRand(RNG normrnd)
134 8835765 : rand(idim : ndim) = rand(idim : ndim) + chol(GET_INDEX(idim : ndim, idim)) * normrnd
135 : end do
136 : end block
137 : #undef GET_INDEX
138 : #undef XPLUS
139 : #endif
140 :
141 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
142 : #elif D2_ENABLED && getMNR_ENABLED
143 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
144 :
145 : #if AM_ENABLED && DC_ENABLED
146 1 : call setMultiNormRand(RNG rand, mean)
147 : #elif DM_ENABLED && AC_ENABLED
148 2 : call setMultiNormRand(RNG rand, chol, subset)
149 : #elif AM_ENABLED && AC_ENABLED
150 0 : call setMultiNormRand(RNG rand, mean, chol, subset)
151 : #else
152 : #error "Unrecognized interface."
153 : #endif
154 :
155 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
156 : #elif D2_ENABLED && setMNR_ENABLED
157 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
158 :
159 : integer(IK) :: ipnt, ndim
160 7 : ndim = size(rand, 1, IK)
161 35007 : do ipnt = 1_IK, size(rand, 2, IK)
162 : #if DM_ENABLED && DC_ENABLED
163 5001 : call setMultiNormRand(RNG rand(1:ndim, ipnt))
164 : #elif AM_ENABLED && DC_ENABLED
165 10002 : call setMultiNormRand(RNG rand(1:ndim, ipnt), mean)
166 : #elif DM_ENABLED && AC_ENABLED
167 20004 : call setMultiNormRand(RNG rand(1:ndim, ipnt), chol, subset)
168 : #elif AM_ENABLED && AC_ENABLED
169 0 : call setMultiNormRand(RNG rand(1:ndim, ipnt), mean, chol, subset)
170 : #else
171 : #error "Unrecognized interface."
172 : #endif
173 : end do
174 :
175 : #else
176 : !%%%%%%%%%%%%%%%%%%%%%%%%
177 : #error "Unrecognized interface."
178 : !%%%%%%%%%%%%%%%%%%%%%%%%
179 : #endif
180 : #undef LOGNORMFAC
181 : #undef MAHAL_SQ
182 : #undef RNG
|
