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 implementations of the tests of [pm_cosmicRate](@ref pm_cosmicRate).
19 : !>
20 : !> \fintest
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, 12:27 AM Tuesday, February 22, 2022, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if getLogRateDensity_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 : real(RKC) , parameter :: TOL = epsilon(0._RKC)*100
32 : real(RKC) :: zplus1, logzplus1
33 : real(RKC) :: logRateDensity_ref
34 : real(RKC) :: logRateDensity
35 : real(RKC) :: diff
36 : integer(IK) :: i
37 35 : assertion = .true._LK
38 17535 : do i = 1_IK, 500_IK
39 17500 : logzplus1 = getUnifRand(0._RKC, 10._RKC)
40 17500 : zplus1 = exp(logzplus1)
41 : #if H06_ENABLED
42 2500 : logRateDensity = getLogRateDensityH06(logzplus1)
43 : #elif L08_ENABLED
44 2500 : logRateDensity = getLogRateDensityL08(logzplus1)
45 : #elif B10_ENABLED
46 2500 : logRateDensity = getLogRateDensityB10(logzplus1)
47 : #elif P15_ENABLED
48 2500 : logRateDensity = getLogRateDensityP15(logzplus1)
49 : #elif M14_ENABLED
50 2500 : logRateDensity = getLogRateDensityM14(zplus1, logzplus1)
51 : #elif M17_ENABLED
52 2500 : logRateDensity = getLogRateDensityM17(zplus1, logzplus1)
53 : #elif F18_ENABLED
54 2500 : logRateDensity = getLogRateDensityF18(zplus1, logzplus1)
55 : #else
56 : #error "Unrecognized interface."
57 : #endif
58 17535 : call report(int(__LINE__, IK))
59 : end do
60 :
61 : contains
62 :
63 : #if H06_ENABLED || L08_ENABLED || B10_ENABLED || P15_ENABLED
64 :
65 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66 :
67 : pure elemental function getLogRateDensity_ref(logzplus1) result(logRateDensity)
68 : real(RKC), intent(in) :: logzplus1
69 : real(RKC) :: logRateDensity
70 : #if H06_ENABLED
71 : real(RKC), parameter :: G0 = +3.4_RKC
72 : real(RKC), parameter :: G1 = -0.3_RKC
73 : real(RKC), parameter :: G2 = -7.8_RKC
74 : real(RKC), parameter :: LOGZ0PLUS1 = log(1._RKC + 0.97_RKC)
75 : real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 4.50_RKC)
76 : real(RKC), parameter :: LOG_NORM_FAC_1 = LOGZ0PLUS1 * (G0 - G1)
77 : real(RKC), parameter :: LOG_NORM_FAC_2 = LOGZ1PLUS1 * (G1 - G2) + LOG_NORM_FAC_1
78 : #elif L08_ENABLED
79 : real(RKC), parameter :: G0 = +3.3000_RKC
80 : real(RKC), parameter :: G1 = +0.0549_RKC
81 : real(RKC), parameter :: G2 = -4.4600_RKC
82 : real(RKC), parameter :: LOGZ0PLUS1 = log(1._RKC + 0.993_RKC)
83 : real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 3.800_RKC)
84 : real(RKC), parameter :: LOG_NORM_FAC_1 = LOGZ0PLUS1 * (G0 - G1)
85 : real(RKC), parameter :: LOG_NORM_FAC_2 = LOGZ1PLUS1 * (G1 - G2) + LOG_NORM_FAC_1
86 : #elif B10_ENABLED
87 : real(RKC), parameter :: G0 = +3.14_RKC
88 : real(RKC), parameter :: G1 = +1.36_RKC
89 : real(RKC), parameter :: G2 = -2.92_RKC
90 : real(RKC), parameter :: LOGZ0PLUS1 = log(1._RKC + 0.97_RKC)
91 : real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 4.00_RKC)
92 : real(RKC), parameter :: LOG_NORM_FAC_1 = LOGZ0PLUS1 * (G0 - G1)
93 : real(RKC), parameter :: LOG_NORM_FAC_2 = LOGZ1PLUS1 * (G1 - G2) + LOG_NORM_FAC_1
94 : #elif P15_ENABLED
95 : real(RKC), parameter :: EXPONENT_HIGH_Z = -7.8_RKC
96 : real(RKC), parameter :: LOGZ1PLUS1 = log(1._RKC + 4.5_RKC)
97 : real(RKC), parameter :: LOG_NORM_FAC_2 = -EXPONENT_HIGH_Z * LOGZ1PLUS1
98 : #else
99 : #error "Unrecognized interface."
100 : #endif
101 : #if H06_ENABLED || L08_ENABLED || B10_ENABLED
102 7500 : if (logzplus1 < LOGZ0PLUS1) then
103 488 : logRateDensity = logzplus1 * G0
104 7012 : elseif (logzplus1 < LOGZ1PLUS1) then
105 698 : logRateDensity = logzplus1 * G1 + LOG_NORM_FAC_1
106 : else
107 6314 : logRateDensity = logzplus1 * G2 + LOG_NORM_FAC_2
108 : end if
109 : #elif P15_ENABLED
110 2500 : if (logzplus1 < LOGZ1PLUS1) then
111 276 : logRateDensity = 0._RKC
112 : else
113 2028 : logRateDensity = logzplus1 * EXPONENT_HIGH_Z + LOG_NORM_FAC_2
114 : end if
115 : #endif
116 6000 : end function
117 : #elif M14_ENABLED || M17_ENABLED || F18_ENABLED
118 7500 : pure elemental function getLogRateDensity_ref(zplus1, logzplus1) result(logRateDensity)
119 : real(RKC), intent(in) :: zplus1, logzplus1
120 : real(RKC) :: logRateDensity
121 : #if M14_ENABLED
122 : real(RKC), parameter :: LOG_AMPLITUDE = log(0.015_RKC)
123 : real(RKC), parameter :: EXPLONENT_LOWER = 2.7_RKC
124 : real(RKC), parameter :: EXPLONENT_UPPER = 5.6_RKC
125 : real(RKC), parameter :: ZPLUS1_BREAK = 2.9_RKC
126 : real(RKC), parameter :: ZPLUS1_COEFF = 1._RKC / (ZPLUS1_BREAK**EXPLONENT_UPPER)
127 : #elif M17_ENABLED
128 : real(RKC), parameter :: LOG_AMPLITUDE = log(0.01_RKC)
129 : real(RKC), parameter :: EXPLONENT_LOWER = 2.6_RKC
130 : real(RKC), parameter :: EXPLONENT_UPPER = 6.2_RKC
131 : real(RKC), parameter :: ZPLUS1_BREAK = 3.2_RKC
132 : real(RKC), parameter :: ZPLUS1_COEFF = 1._RKC / (ZPLUS1_BREAK**EXPLONENT_UPPER)
133 : #elif F18_ENABLED
134 : real(RKC), parameter :: LOG_AMPLITUDE = log(0.013_RKC)
135 : real(RKC), parameter :: EXPLONENT_LOWER = 2.99_RKC
136 : real(RKC), parameter :: EXPLONENT_UPPER = 6.19_RKC
137 : real(RKC), parameter :: ZPLUS1_BREAK = 2.63_RKC
138 : real(RKC), parameter :: ZPLUS1_COEFF = 1._RKC / (ZPLUS1_BREAK**EXPLONENT_UPPER)
139 : #endif
140 7500 : logRateDensity = LOG_AMPLITUDE + EXPLONENT_LOWER * logzplus1 - log(1._RKC + ZPLUS1_COEFF * zplus1**EXPLONENT_UPPER)
141 7500 : end function
142 : #else
143 : #error "Unrecognized interface."
144 : #endif
145 :
146 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147 :
148 17500 : subroutine report(line)
149 : integer(IK), intent(in) :: line
150 : #if H06_ENABLED || L08_ENABLED || B10_ENABLED || P15_ENABLED
151 10000 : logRateDensity_ref = getLogRateDensity_ref(logzplus1)
152 : #elif M14_ENABLED || M17_ENABLED || F18_ENABLED
153 7500 : logRateDensity_ref = getLogRateDensity_ref(zplus1, logzplus1)
154 : #else
155 : #error "Unrecognized interface."
156 : #endif
157 17500 : diff = abs(logRateDensity - logRateDensity_ref)
158 17500 : assertion = assertion .and. logical(diff <= TOL * abs(logRateDensity_ref), LK)
159 17500 : if (test%traceable .and. .not. assertion) then
160 : ! LCOV_EXCL_START
161 : write(test%disp%unit,"(*(g0,:,', '))")
162 : write(test%disp%unit,"(*(g0,:,', '))") "logRateDensity_ref ", logRateDensity_ref
163 : write(test%disp%unit,"(*(g0,:,', '))") "logRateDensity ", logRateDensity
164 : write(test%disp%unit,"(*(g0,:,', '))") "logzplus1 ", logzplus1
165 : write(test%disp%unit,"(*(g0,:,', '))") "zplus1 ", zplus1
166 : write(test%disp%unit,"(*(g0,:,', '))") "diff ", diff
167 : write(test%disp%unit,"(*(g0,:,', '))") "TOL ", TOL
168 : write(test%disp%unit,"(*(g0,:,', '))")
169 : ! LCOV_EXCL_STOP
170 : end if
171 17500 : call test%assert(assertion, desc = "The log rate density must be computed correctly for the specified logzplus1.", line = line)
172 17500 : end subroutine
173 :
174 : #else
175 : !%%%%%%%%%%%%%%%%%%%%%%%%
176 : #error "Unrecognized interface."
177 : !%%%%%%%%%%%%%%%%%%%%%%%%
178 : #endif
|
