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 the implementation details of the routines under the generic interface [pm_sampleECDF](@ref pm_sampleECDF).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Saturday 4:40 PM, August 21, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! The accumulation method accumulates significant error for large arrays.
28 : ! Either multiplication must be used, or the accurate summation method of fabcomp().
29 : ! The current accumulation workaround is ugly and needs a revamp using fabcomp().
30 : #define ACCUMULATION_ENABLED 1
31 : #define MULTIPLICATION_ENABLED 0
32 :
33 : !%%%%%%%%%%%%%%
34 : #if setECDF_ENABLED
35 : !%%%%%%%%%%%%%%
36 :
37 : #define CHECK_LEN_LCDF \
38 : CHECK_ASSERTION(__LINE__, size(lcdf, 1, IK) == size(ecdf, 1, IK), \
39 : SK_"@setECDF(): The condition `size(lcdf) == size(ecdf)` must hold. size(lcdf), size(ecdf) = "//getStr([size(lcdf, 1, IK), size(ecdf, 1, IK)]))
40 : #define CHECK_LEN_UCDF \
41 : CHECK_ASSERTION(__LINE__, size(ucdf, 1, IK) == size(ecdf, 1, IK), \
42 : SK_"@setECDF(): The condition `size(ucdf) == size(ecdf)` must hold. size(ucdf), size(ecdf) = "//getStr([size(ucdf, 1, IK), size(ecdf, 1, IK)]))
43 :
44 : integer(IK) :: nsam, isam
45 : real(TKC), parameter :: ZERO = 0._TKC, UNIT = 1._TKC
46 : real(TKC):: nsamInv
47 171 : nsam = size(ecdf, 1, IK)
48 171 : if (nsam == 0_IK) return
49 : #if ONE_ENABLED
50 : #define GET_WEIGHTED(W,X)X
51 57 : nsamInv = UNIT / real(nsam, TKC)
52 : #elif WIK_ENABLED || WRK_ENABLED
53 : #define GET_WEIGHTED(W,X)real(W, TKC) * X
54 1329704 : CHECK_ASSERTION(__LINE__, all(0 <= weight), SK_"@setECDF(): The condition `all(0 <= weight)` must hold. weight = "//getStr(weight))
55 2659522 : CHECK_ASSERTION(__LINE__, abs(weisum - sum(weight)) < epsilon(0._TKC) * 100, SK_"@setECDF(): The condition `abs(weisum - sum(weight)) < epsilon(0._TKC) * 100` must hold. weisum, sum(weight) = "//getStr([weisum, sum(weight)]))
56 114 : nsamInv = UNIT / real(weisum, TKC)
57 : #else
58 : #error "Unrecognized interface."
59 : #endif
60 171 : if (1_IK < nsam) then
61 29 : ecdf(1) = GET_WEIGHTED(weight(1),nsamInv)
62 : #if !WRK_ENABLED
63 138 : blockExpectedLowerUpper: if (.not. (present(lcdf) .or. present(ucdf))) then
64 : #endif
65 : #if ACCUMULATION_ENABLED || WIK_ENABLED || WRK_ENABLED
66 : #define ISAM isam
67 1329644 : do isam = 2, nsam - 1
68 1329644 : ecdf(isam) = ecdf(isam - 1) + GET_WEIGHTED(weight(isam),nsamInv)
69 : end do
70 145 : if (1._TKC < ecdf(nsam - 1)) then
71 4 : nsamInv = 1._TKC / (ecdf(nsam - 1) + GET_WEIGHTED(weight(nsam),nsamInv))
72 : do concurrent(isam = 1 : nsam - 1)
73 39972 : ecdf(isam) = ecdf(isam) * nsamInv
74 : end do
75 : end if
76 : #elif MULTIPLICATION_ENABLED
77 : #define ISAM nsam
78 : do concurrent(isam = 2 : nsam - 1)
79 : ecdf(isam) = real(isam, TKC) * nsamInv
80 : end do
81 : #else
82 : #error "Unrecognized interface."
83 : #endif
84 145 : ecdf(ISAM) = UNIT
85 : #if !WRK_ENABLED
86 : else blockExpectedLowerUpper
87 : block
88 : real(TKC):: sqrt_nsamInvTwice_logTwoOverAlpha
89 22 : if (present(alpha)) then
90 4 : CHECK_ASSERTION(__LINE__, ZERO < alpha .and. alpha < UNIT, SK_"@setECDF(): The condition `0 < alpha .and. alpha < 1` must hold. alpha = "//getStr(alpha))
91 4 : sqrt_nsamInvTwice_logTwoOverAlpha = sqrt(0.5_TKC * nsamInv * log(2._TKC / alpha))
92 : else
93 18 : sqrt_nsamInvTwice_logTwoOverAlpha = sqrt(0.5_TKC * nsamInv * log(2._TKC / 0.05_TKC))
94 : end if
95 22 : if (present(lcdf) .and. present(ucdf)) then
96 66 : CHECK_LEN_LCDF ! fpp
97 66 : CHECK_LEN_UCDF ! fpp
98 22 : lcdf(1) = max(ZERO, GET_WEIGHTED(weight(1),nsamInv) - sqrt_nsamInvTwice_logTwoOverAlpha)
99 22 : ucdf(1) = min(UNIT, GET_WEIGHTED(weight(1),nsamInv) + sqrt_nsamInvTwice_logTwoOverAlpha)
100 : #if ACCUMULATION_ENABLED || WIK_ENABLED || WRK_ENABLED
101 11216 : do isam = 2, nsam - 1
102 11194 : ecdf(isam) = ecdf(isam - 1) + GET_WEIGHTED(weight(isam),nsamInv)
103 11194 : lcdf(isam) = max(ZERO, ecdf(isam) - sqrt_nsamInvTwice_logTwoOverAlpha)
104 11216 : ucdf(isam) = min(UNIT, ecdf(isam) + sqrt_nsamInvTwice_logTwoOverAlpha)
105 : end do
106 : #elif MULTIPLICATION_ENABLED
107 : do concurrent(isam = 2 : nsam - 1)
108 : ecdf(isam) = real(isam, TKC) * nsamInv
109 : lcdf(isam) = max(ZERO, ecdf(isam) - sqrt_nsamInvTwice_logTwoOverAlpha)
110 : ucdf(isam) = min(UNIT, ecdf(isam) + sqrt_nsamInvTwice_logTwoOverAlpha)
111 : end do
112 : #endif
113 22 : ecdf(ISAM) = UNIT
114 22 : lcdf(ISAM) = max(ZERO, UNIT - sqrt_nsamInvTwice_logTwoOverAlpha)
115 22 : ucdf(ISAM) = UNIT
116 0 : elseif (present(lcdf)) then
117 0 : CHECK_LEN_LCDF ! fpp
118 0 : lcdf(1) = max(ZERO, GET_WEIGHTED(weight(1),nsamInv) - sqrt_nsamInvTwice_logTwoOverAlpha)
119 : #if ACCUMULATION_ENABLED || WIK_ENABLED || WRK_ENABLED
120 0 : do isam = 2, nsam - 1
121 0 : ecdf(isam) = ecdf(isam - 1) + GET_WEIGHTED(weight(isam),nsamInv)
122 0 : lcdf(isam) = max(ZERO, ecdf(isam) - sqrt_nsamInvTwice_logTwoOverAlpha)
123 : end do
124 : #elif MULTIPLICATION_ENABLED
125 : do concurrent(isam = 2 : nsam - 1)
126 : ecdf(isam) = real(isam, TKC) * nsamInv
127 : lcdf(isam) = max(ZERO, ecdf(isam) - sqrt_nsamInvTwice_logTwoOverAlpha)
128 : end do
129 : #endif
130 0 : ecdf(ISAM) = UNIT
131 0 : lcdf(ISAM) = max(ZERO, UNIT - sqrt_nsamInvTwice_logTwoOverAlpha)
132 0 : elseif (present(ucdf)) then
133 0 : CHECK_LEN_UCDF ! fpp
134 0 : ucdf(1) = min(UNIT, GET_WEIGHTED(weight(1),nsamInv) + sqrt_nsamInvTwice_logTwoOverAlpha)
135 : #if ACCUMULATION_ENABLED || WIK_ENABLED || WRK_ENABLED
136 0 : do isam = 2, nsam - 1
137 0 : ecdf(isam) = ecdf(isam - 1) + GET_WEIGHTED(weight(isam),nsamInv)
138 0 : ucdf(isam) = min(UNIT, ecdf(isam) + sqrt_nsamInvTwice_logTwoOverAlpha)
139 : end do
140 : #elif MULTIPLICATION_ENABLED
141 : do concurrent(isam = 2 : nsam - 1)
142 : ecdf(isam) = real(isam, TKC) * nsamInv
143 : ucdf(isam) = min(UNIT, ecdf(isam) + sqrt_nsamInvTwice_logTwoOverAlpha)
144 : end do
145 : #endif
146 0 : ecdf(ISAM) = UNIT
147 0 : ucdf(ISAM) = UNIT
148 : end if
149 : end block
150 : end if blockExpectedLowerUpper
151 : #endif
152 : else
153 : ! nsam == 1
154 8 : ecdf = UNIT
155 : #if !WRK_ENABLED
156 5 : if (present(lcdf)) lcdf = UNIT
157 5 : if (present(ucdf)) ucdf = UNIT
158 : #endif
159 : end if
160 : #else
161 : !%%%%%%%%%%%%%%%%%%%%%%%%
162 : #error "Unrecognized interface."
163 : !%%%%%%%%%%%%%%%%%%%%%%%%
164 : #endif
165 : #undef CHECK_LEN_LCDF
166 : #undef CHECK_LEN_UCDF
167 : #undef GET_WEIGHTED
168 : #undef ISAM
|
