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 implementations of procedures of [pm_distBand](@ref pm_distBand).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #define CHECK_ALPHA_NEQ_TWO \
28 : CHECK_ASSERTION(__LINE__, alpha /= -2._RKC, SK_": The condition `alpha /= -2.` must hold. alpha = "//getStr(alpha))
29 : #define CHECK_ALPHA_GEQ_BETA \
30 : CHECK_ASSERTION(__LINE__, beta < alpha, SK_": The condition `beta < alpha` must hold. beta, alpha = "//getStr([beta, alpha]))
31 : #define CHECK_EBREAK_GEQ_ZERO \
32 : CHECK_ASSERTION(__LINE__, 0._RKC < ebreak, SK_": The condition `0. < ebreak` must hold. ebreak = "//getStr(ebreak))
33 :
34 : !%%%%%%%%%%%%%%%%%%%
35 : #if getBandEpeak_ENABLED
36 : !%%%%%%%%%%%%%%%%%%%
37 :
38 6 : CHECK_ALPHA_NEQ_TWO
39 18 : CHECK_ALPHA_GEQ_BETA
40 6 : CHECK_EBREAK_GEQ_ZERO
41 6 : epeak = ebreak * (alpha + 2._RKC) / (alpha - beta)
42 :
43 : !%%%%%%%%%%%%%%%%%%%%
44 : #elif getBandEbreak_ENABLED
45 : !%%%%%%%%%%%%%%%%%%%%
46 :
47 12317 : CHECK_ALPHA_NEQ_TWO
48 36951 : CHECK_ALPHA_GEQ_BETA
49 12317 : CHECK_ASSERTION(__LINE__, 0._RKC < epeak, SK_"@getBandEbreak(): The condition `0. < epeak` must hold. epeak = "//getStr(epeak))
50 12317 : ebreak = epeak * (alpha - beta) / (alpha + 2._RKC)
51 :
52 : !%%%%%%%%%%%%%%%%%%
53 : #elif getBandZeta_ENABLED
54 : !%%%%%%%%%%%%%%%%%%
55 :
56 31239 : CHECK_EBREAK_GEQ_ZERO
57 31239 : zeta = ebreak**(alpha - beta) * exp(beta - alpha)
58 : !logZeta = (alpha - beta) * (ebreak - 1._RKC)
59 :
60 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 : #elif getBandUDF_ENABLED && Any_ENABLED
62 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 :
64 3998 : CHECK_ALPHA_NEQ_TWO
65 11994 : CHECK_ALPHA_GEQ_BETA
66 3998 : CHECK_EBREAK_GEQ_ZERO
67 3998 : CHECK_ASSERTION(__LINE__, 0._RKC < energy, SK_"@getBandUDF(): The condition `0. < energy` must hold. energy = "//getStr(energy))
68 3998 : if (energy < ebreak) then
69 848 : if (present(invEfold)) then
70 7 : CHECK_ASSERTION(__LINE__, abs(invEfold - (alpha - beta) / ebreak) < 10 * epsilon(energy), \
71 : SK_"@getBandUDF(): The condition `abs(invEfold - (alpha - beta) / ebreak) < 10 * epsilon(energy)` must hold. alpha, beta, ebreak, invEfold, (alpha - beta) / ebreak, epsilon(energy) = "//\
72 : getStr([alpha, beta, ebreak, invEfold, (alpha - beta) / ebreak, epsilon(energy)]))
73 1 : udf = energy**alpha * exp(-energy * invEfold)
74 : else
75 847 : udf = energy**alpha * exp((beta - alpha) * (energy / ebreak))
76 : end if
77 3150 : elseif (present(zeta)) then
78 0 : CHECK_ASSERTION(__LINE__, abs(zeta - getBandZeta(alpha, beta, ebreak)) < 10 * epsilon(energy), \
79 : SK_"@getBandUDF(): The condition `abs(zeta - getBandZeta(alpha, beta, ebreak)) < 10 * epsilon(energy)` must hold. alpha, beta, ebreak, zeta, getBandZeta(alpha, beta, ebreak), epsilon(energy) = "//\
80 : getStr([alpha, beta, ebreak, zeta, getBandZeta(alpha, beta, ebreak), epsilon(energy)]))
81 0 : udf = zeta * energy**beta
82 : else
83 3150 : udf = getBandZeta(alpha, beta, ebreak) * energy**beta
84 : end if
85 :
86 : !%%%%%%%%%%%%%%%%%%
87 : #elif setBandUCDF_ENABLED
88 : !%%%%%%%%%%%%%%%%%%
89 :
90 : real(RKC) :: mb, tmp, invEfold
91 :
92 33239 : CHECK_ALPHA_NEQ_TWO
93 99717 : CHECK_ALPHA_GEQ_BETA
94 33239 : CHECK_EBREAK_GEQ_ZERO
95 :
96 33239 : info = 0_IK
97 33239 : ucdf = 0._RKC
98 33239 : if (lb < ub) then ! integrate the spectrum
99 33235 : CHECK_ASSERTION(__LINE__, 0._RKC < lb, SK_"@setBandUCDF(): The condition `0. < lb` must hold. lb = "//getStr(lb))
100 33235 : if (lb < ebreak) then
101 : ! First compute the contribution from the lower component via either Gamma CDF or brute-force integration.
102 28563 : mb = min(ub, ebreak)
103 28563 : invEfold = (alpha - beta) / ebreak
104 28563 : if (-1._RKC < alpha) then
105 : ! use Gamma CDF.
106 : block
107 : real(RKC) :: kappa, dumm
108 23098 : kappa = alpha + 1._RKC
109 23098 : dumm = log_gamma(kappa)
110 23098 : call setGammaCDF(ucdf, lb, dumm, kappa, invEfold, info)
111 : if (info < 0_IK) return ! LCOV_EXCL_LINE
112 23098 : call setGammaCDF(tmp, mb, dumm, kappa, invEfold, info)
113 : if (info < 0_IK) return ! LCOV_EXCL_LINE
114 23098 : ucdf = (tmp - ucdf) * gamma(kappa) / invEfold**kappa
115 : end block
116 : else
117 : block
118 : integer(IK) :: neval, nint, sindex(1000)
119 : real(RKC) :: abstol, reltol, abserr, sinfo(4, 1000)
120 5465 : abstol = 0._RKC
121 5465 : reltol = epsilon(0._RKC)**(2./3.)
122 5465 : info = getQuadErr(getBandLowUDF, lb, mb, abstol, reltol, GK21, weps, ucdf, abserr, sinfo, sindex, neval, nint)
123 : if (info /= 0_IK) info = -info ! LCOV_EXCL_LINE
124 : if (info < 0_IK) return ! LCOV_EXCL_LINE
125 : end block
126 : end if
127 28563 : if (mb == ub) return
128 : else
129 4672 : mb = lb
130 : end if
131 : ! Add the remaining part of the ucdf from the high-energy component.
132 28083 : tmp = beta + 1._RKC
133 28083 : if (tmp /= 0._RKC) then
134 23871 : ucdf = ucdf + getBandZeta(alpha, beta, ebreak) * (ub**tmp - mb**tmp) / tmp
135 : else
136 4212 : ucdf = ucdf + getBandZeta(alpha, beta, ebreak) * log(ub / mb)
137 : end if
138 : end if
139 :
140 : contains
141 :
142 1401351 : PURE function getBandLowUDF(energy) result(udf)
143 : implicit none
144 : real(RKC), intent(in) :: energy
145 : real(RKC) :: udf
146 5605404 : CHECK_ASSERTION(__LINE__, lb <= energy .and. energy <= ub, SK_"@setBandUCDF(): The condition `lb <= energy .and. energy <= ub` must hold. lb, energy, ub = "//getStr([lb, energy, ub]))
147 1401351 : udf = energy**alpha * exp(-invEfold * energy)
148 : !print *, udf, invEfold, energy
149 1401351 : end function
150 :
151 : !%%%%%%%%%%%%%%%%%%
152 : #elif setBandMean_ENABLED
153 : !%%%%%%%%%%%%%%%%%%
154 :
155 : #if Def_ENABLED
156 : #define LBNEW lb
157 : #define UBNEW ub
158 : #elif !New_ENABLED
159 : #error "Unrecognized interface."
160 : #endif
161 : real(RKC) :: denom
162 1 : call setBandUCDF(denom, lb, ub, alpha, beta, ebreak, info)
163 : if (info < 0_IK) return ! LCOV_EXCL_LINE
164 1 : call setBandUCDF(mean, LBNEW, UBNEW, alpha + 1._RKC, beta + 1._RKC, ebreak, info)
165 : if (info < 0_IK) return ! LCOV_EXCL_LINE
166 1 : mean = mean / denom
167 : #undef LBNEW
168 : #undef UBNEW
169 :
170 : !%%%%%%%%%%%%%%%%%%%%
171 : #elif setBandPhoton_ENABLED
172 : !%%%%%%%%%%%%%%%%%%%%
173 :
174 : #if FromPhoton_ENABLED && NewB_ENABLED
175 : real(RKC) :: denom, numer
176 0 : CHECK_ASSERTION(__LINE__, lb <= ub, SK_"@setBandPhoton(): The condition `lb <= ub` must hold. lb, ub = "//getStr([lb, ub]))
177 0 : CHECK_ASSERTION(__LINE__, lbnew <= ubnew, SK_"@setBandPhoton(): The condition `lbnew <= ubnew` must hold. lbnew, ubnew = "//getStr([lbnew, ubnew]))
178 0 : CHECK_ASSERTION(__LINE__, 0._RKC < photon, SK_"@setBandPhoton(): The condition `0. < photon` must hold. photon = "//getStr(photon))
179 0 : call setBandUCDF(denom, lb, ub, alpha, beta, ebreak, info)
180 : if (info < 0_IK) return ! LCOV_EXCL_LINE
181 0 : call setBandUCDF(numer, lbnew, ubnew, alpha, beta, ebreak, info)
182 : if (info < 0_IK) return ! LCOV_EXCL_LINE
183 0 : photon = numer * (photon / denom)
184 : #elif FromEnergy_ENABLED
185 : #if OldB_ENABLED
186 : #define LBNEW lb
187 : #define UBNEW ub
188 : #elif !NewB_ENABLED
189 : #error "Unrecognized interface."
190 : #endif
191 : real(RKC) :: denom, numer
192 5001 : CHECK_ASSERTION(__LINE__, 0._RKC < energy, SK_"@setBandPhoton(): The condition `0. < energy` must hold. energy = "//getStr(energy))
193 5001 : call setBandUCDF(denom, lb, ub, alpha + 1._RKC, beta + 1._RKC, ebreak, info)
194 : if (info < 0_IK) return ! LCOV_EXCL_LINE
195 5001 : call setBandUCDF(numer, LBNEW, UBNEW, alpha, beta, ebreak, info)
196 : if (info < 0_IK) return ! LCOV_EXCL_LINE
197 5001 : photon = numer * (energy / denom)
198 : #undef LBNEW
199 : #undef UBNEW
200 : #else
201 : #error "Unrecognized interface."
202 : #endif
203 :
204 : !%%%%%%%%%%%%%%%%%%%%
205 : #elif setBandEnergy_ENABLED
206 : !%%%%%%%%%%%%%%%%%%%%
207 :
208 : #if FromEnergy_ENABLED && NewB_ENABLED
209 : real(RKC) :: denom, numer
210 6924 : CHECK_ASSERTION(__LINE__, lb <= ub, SK_"@setBandPhoton(): The condition `lb <= ub` must hold. lb, ub = "//getStr([lb, ub]))
211 6924 : CHECK_ASSERTION(__LINE__, lbnew <= ubnew, SK_"@setBandPhoton(): The condition `lbnew <= ubnew` must hold. lbnew, ubnew = "//getStr([lbnew, ubnew]))
212 2308 : CHECK_ASSERTION(__LINE__, 0._RKC < energy, SK_"@setBandPhoton(): The condition `0. < energy` must hold. energy = "//getStr(energy))
213 2308 : call setBandUCDF(denom, lb, ub, alpha + 1._RKC, beta + 1._RKC, ebreak, info)
214 : if (info < 0_IK) return ! LCOV_EXCL_LINE
215 2308 : call setBandUCDF(numer, lbnew, ubnew, alpha + 1._RKC, beta + 1._RKC, ebreak, info)
216 : if (info < 0_IK) return ! LCOV_EXCL_LINE
217 2308 : energy = numer * (energy / denom)
218 : #elif FromPhoton_ENABLED
219 : #if OldB_ENABLED
220 : #define LBNEW lb
221 : #define UBNEW ub
222 : #elif !NewB_ENABLED
223 : #error "Unrecognized interface."
224 : #endif
225 : real(RKC) :: denom, numer
226 7309 : CHECK_ASSERTION(__LINE__, 0._RKC < photon, SK_"@setBandPhoton(): The condition `0. < photon` must hold. photon = "//getStr(photon))
227 7309 : call setBandUCDF(denom, lb, ub, alpha, beta, ebreak, info)
228 : if (info < 0_IK) return ! LCOV_EXCL_LINE
229 7309 : call setBandUCDF(numer, LBNEW, UBNEW, alpha + 1._RKC, beta + 1._RKC, ebreak, info)
230 : if (info < 0_IK) return ! LCOV_EXCL_LINE
231 7309 : energy = numer * (photon / denom)
232 : #undef LBNEW
233 : #undef UBNEW
234 : #else
235 : #error "Unrecognized interface."
236 : #endif
237 :
238 : #else
239 : !%%%%%%%%%%%%%%%%%%%%%%%%
240 : #error "Unrecognized interface."
241 : !%%%%%%%%%%%%%%%%%%%%%%%%
242 : #endif
243 : #undef CHECK_EBREAK_GEQ_ZERO
244 : #undef CHECK_ALPHA_GEQ_BETA
245 : #undef CHECK_ALPHA_NEQ_TWO
|
