Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!!
4 : !!!! MIT License
5 : !!!!
6 : !!!! ParaMonte: plain powerful parallel Monte Carlo library.
7 : !!!!
8 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab
9 : !!!!
10 : !!!! This file is part of the ParaMonte library.
11 : !!!!
12 : !!!! Permission is hereby granted, free of charge, to any person obtaining a
13 : !!!! copy of this software and associated documentation files (the "Software"),
14 : !!!! to deal in the Software without restriction, including without limitation
15 : !!!! the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : !!!! and/or sell copies of the Software, and to permit persons to whom the
17 : !!!! Software is furnished to do so, subject to the following conditions:
18 : !!!!
19 : !!!! The above copyright notice and this permission notice shall be
20 : !!!! included in all copies or substantial portions of the Software.
21 : !!!!
22 : !!!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23 : !!!! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24 : !!!! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
25 : !!!! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
26 : !!!! DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
27 : !!!! OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
28 : !!!! OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
29 : !!!!
30 : !!!! ACKNOWLEDGMENT
31 : !!!!
32 : !!!! ParaMonte is an honor-ware and its currency is acknowledgment and citations.
33 : !!!! As per the ParaMonte library license agreement terms, if you use any parts of
34 : !!!! this library for any purposes, kindly acknowledge the use of ParaMonte in your
35 : !!!! work (education/research/industry/development/...) by citing the ParaMonte
36 : !!!! library as described on this page:
37 : !!!!
38 : !!!! https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
39 : !!!!
40 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42 :
43 : !> \brief This module contains tests of the module [Statistics_mod](@ref statistics_mod).
44 : !> \author Amir Shahmoradi
45 :
46 : module Test_Statistics_mod
47 :
48 : use Statistics_mod
49 : use Err_mod, only: Err_type
50 : use Test_mod, only: Test_type
51 : implicit none
52 :
53 : private
54 : public :: test_Statistics
55 :
56 : type(Test_type) :: Test
57 :
58 : integer(IK) , parameter :: lenRnd = 50_IK
59 :
60 : real(RK) :: UnifRnd(lenRnd) = [ 0.0759666916908419_RK &
61 : , 0.2399161535536580_RK &
62 : , 0.1233189348351660_RK &
63 : , 0.1839077882824170_RK &
64 : , 0.2399525256649030_RK &
65 : , 0.4172670690843700_RK &
66 : , 0.0496544303257421_RK &
67 : , 0.9027161099152810_RK &
68 : , 0.9447871897216460_RK &
69 : , 0.4908640924680800_RK &
70 : , 0.4892526384000190_RK &
71 : , 0.3377194098213770_RK &
72 : , 0.9000538464176620_RK &
73 : , 0.3692467811202150_RK &
74 : , 0.1112027552937870_RK &
75 : , 0.7802520683211380_RK &
76 : , 0.3897388369612530_RK &
77 : , 0.2416912859138330_RK &
78 : , 0.4039121455881150_RK &
79 : , 0.0964545251683886_RK &
80 : , 0.1319732926063350_RK &
81 : , 0.9420505907754850_RK &
82 : , 0.9561345402298020_RK &
83 : , 0.5752085950784660_RK &
84 : , 0.0597795429471558_RK &
85 : , 0.2347799133724060_RK &
86 : , 0.3531585712220710_RK &
87 : , 0.8211940401979590_RK &
88 : , 0.0154034376515551_RK &
89 : , 0.0430238016578078_RK &
90 : , 0.1689900294627040_RK &
91 : , 0.6491154749564520_RK &
92 : , 0.7317223856586700_RK &
93 : , 0.6477459631363070_RK &
94 : , 0.4509237064309450_RK &
95 : , 0.5470088922863450_RK &
96 : , 0.2963208056077730_RK &
97 : , 0.7446928070741560_RK &
98 : , 0.1889550150325450_RK &
99 : , 0.6867754333653150_RK &
100 : , 0.1835111557372700_RK &
101 : , 0.3684845964903370_RK &
102 : , 0.6256185607296900_RK &
103 : , 0.7802274351513770_RK &
104 : , 0.0811257688657853_RK &
105 : , 0.9293859709687300_RK &
106 : , 0.7757126786084020_RK &
107 : , 0.4867916324031720_RK &
108 : , 0.4358585885809190_RK &
109 : , 0.4467837494298060_RK ]
110 :
111 : real(RK) :: StdNormRnd1(lenRnd) = [ 0.537667139546100_RK &
112 : , 1.83388501459509_RK &
113 : , -2.25884686100365_RK &
114 : , 0.862173320368121_RK &
115 : , 0.318765239858981_RK &
116 : , -1.30768829630527_RK &
117 : , -0.433592022305684_RK &
118 : , 0.342624466538650_RK &
119 : , 3.57839693972576_RK &
120 : , 2.76943702988488_RK &
121 : , -1.34988694015652_RK &
122 : , 3.03492346633185_RK &
123 : , 0.725404224946106_RK &
124 : , -0.0630548731896562_RK &
125 : , 0.714742903826096_RK &
126 : , -0.204966058299775_RK &
127 : , -0.124144348216312_RK &
128 : , 1.48969760778547_RK &
129 : , 1.40903448980048_RK &
130 : , 1.41719241342961_RK &
131 : , 0.671497133608081_RK &
132 : , -1.20748692268504_RK &
133 : , 0.717238651328839_RK &
134 : , 1.63023528916473_RK &
135 : , 0.488893770311789_RK &
136 : , 1.03469300991786_RK &
137 : , 0.726885133383238_RK &
138 : , -0.303440924786016_RK &
139 : , 0.293871467096658_RK &
140 : , -0.787282803758638_RK &
141 : , 0.888395631757642_RK &
142 : , -1.14707010696915_RK &
143 : , -1.06887045816803_RK &
144 : , -0.809498694424876_RK &
145 : , -2.94428416199490_RK &
146 : , 1.43838029281510_RK &
147 : , 0.325190539456198_RK &
148 : , -0.754928319169703_RK &
149 : , 1.37029854009523_RK &
150 : , -1.71151641885370_RK &
151 : , -0.102242446085491_RK &
152 : , -0.241447041607358_RK &
153 : , 0.319206739165502_RK &
154 : , 0.312858596637428_RK &
155 : , -0.864879917324457_RK &
156 : , -0.0300512961962686_RK &
157 : , -0.164879019209038_RK &
158 : , 0.627707287528727_RK &
159 : , 1.09326566903948_RK &
160 : , 1.10927329761440_RK ]
161 :
162 : real(RK) :: StdNormRnd2(lenRnd) = [ -0.863652821988714_RK &
163 : , 0.0773590911304249_RK &
164 : , -1.21411704361541_RK &
165 : , -1.11350074148676_RK &
166 : , -0.00684932810334806_RK &
167 : , 1.53263030828475_RK &
168 : , -0.769665913753682_RK &
169 : , 0.371378812760058_RK &
170 : , -0.225584402271252_RK &
171 : , 1.11735613881447_RK &
172 : , -1.08906429505224_RK &
173 : , 0.0325574641649735_RK &
174 : , 0.552527021112224_RK &
175 : , 1.10061021788087_RK &
176 : , 1.54421189550395_RK &
177 : , 0.0859311331754255_RK &
178 : , -1.49159031063761_RK &
179 : , -0.742301837259857_RK &
180 : , -1.06158173331999_RK &
181 : , 2.35045722400204_RK &
182 : , -0.615601881466894_RK &
183 : , 0.748076783703985_RK &
184 : , -0.192418510588264_RK &
185 : , 0.888610425420721_RK &
186 : , -0.764849236567874_RK &
187 : , -1.40226896933876_RK &
188 : , -1.42237592509150_RK &
189 : , 0.488193909859941_RK &
190 : , -0.177375156618825_RK &
191 : , -0.196053487807333_RK &
192 : , 1.41931015064255_RK &
193 : , 0.291584373984183_RK &
194 : , 0.197811053464361_RK &
195 : , 1.58769908997406_RK &
196 : , -0.804465956349547_RK &
197 : , 0.696624415849607_RK &
198 : , 0.835088165072682_RK &
199 : , -0.243715140377952_RK &
200 : , 0.215670086403744_RK &
201 : , -1.16584393148205_RK &
202 : , -1.14795277889859_RK &
203 : , 0.104874716016494_RK &
204 : , 0.722254032225002_RK &
205 : , 2.58549125261624_RK &
206 : , -0.666890670701386_RK &
207 : , 0.187331024578940_RK &
208 : , -0.0824944253709554_RK &
209 : , -1.93302291785099_RK &
210 : , -0.438966153934773_RK &
211 : , -1.79467884145512_RK ]
212 :
213 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
214 :
215 : contains
216 :
217 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
218 :
219 3 : subroutine test_Statistics()
220 : implicit none
221 3 : Test = Test_type(moduleName=MODULE_NAME)
222 3 : call Test%run(test_doKS1_1, "test_doKS1_1")
223 3 : call Test%run(test_flatten_1, "test_flatten_1")
224 3 : call Test%run(test_getMVNDev_1, "test_getMVNDev_1")
225 3 : call Test%run(test_getMVUDev_1, "test_getMVUDev_1")
226 3 : call Test%run(test_getHist1D_1, "test_getHist1D_1")
227 3 : call Test%run(test_getHist1D_2, "test_getHist1D_2")
228 3 : call Test%run(test_getHist1D_3, "test_getHist1D_3")
229 3 : call Test%run(test_getHist2D_1, "test_getHist2D_1")
230 3 : call Test%run(test_getHist2D_2, "test_getHist2D_2")
231 3 : call Test%run(test_getHist2D_3, "test_getHist2D_3")
232 3 : call Test%run(test_getHist2D_4, "test_getHist2D_4")
233 3 : call Test%run(test_getHist2D_5, "test_getHist2D_5")
234 3 : call Test%run(test_getRandMVN_1, "test_getRandMVN_1")
235 3 : call Test%run(test_getRandMVU_1, "test_getRandMVU_1")
236 3 : call Test%run(test_getMean_2D_1, "test_getMean_2D_1")
237 3 : call Test%run(test_getMean_2D_2, "test_getMean_2D_2")
238 3 : call Test%run(test_getRandExp_1, "test_getRandExp_1")
239 3 : call Test%run(test_getNormPDF_1, "test_getNormPDF_1")
240 3 : call Test%run(test_getNormCDF_1, "test_getNormCDF_1")
241 3 : call Test%run(test_getQuantile_1, "test_getQuantile_1")
242 3 : call Test%run(test_getQuantile_2, "test_getQuantile_2")
243 3 : call Test%run(test_getRandGaus_1, "test_getRandGaus_1")
244 3 : call Test%run(test_getRandNorm_1, "test_getRandNorm_1")
245 3 : call Test%run(test_getRandLogn_1, "test_getRandLogn_1")
246 3 : call Test%run(test_getRandBeta_1, "test_getRandBeta_1")
247 3 : call Test%run(test_getSNormPDF_1, "test_getSNormPDF_1")
248 3 : call Test%run(test_doSortedKS2_1, "test_doSortedKS2_1")
249 3 : call Test%run(test_doUniformKS1_1, "test_doUniformKS1_1")
250 3 : call Test%run(test_mergeMeanCov_1, "test_mergeMeanCov_1")
251 3 : call Test%run(test_getRandGamma_1, "test_getRandGamma_1")
252 3 : call Test%run(test_getLogProbMVU_1, "test_getLogProbMVU_1")
253 3 : call Test%run(test_getSamCholFac_1, "test_getSamCholFac_1")
254 3 : call Test%run(test_getSamCovMean_1, "test_getSamCovMean_1")
255 3 : call Test%run(test_getRandCorMat_1, "test_getRandCorMat_1")
256 3 : call Test%run(test_getLogProbGeo_1, "test_getLogProbGeo_1")
257 3 : call Test%run(test_getUniformCDF_1, "test_getUniformCDF_1")
258 3 : call Test%run(test_getRandUniform_1, "test_getRandUniform_1")
259 3 : call Test%run(test_getNormData_2D_1, "test_getNormData_2D_1")
260 3 : call Test%run(test_getVariance_1D_1, "test_getVariance_1D_1")
261 3 : call Test%run(test_getVariance_1D_2, "test_getVariance_1D_2")
262 3 : call Test%run(test_getVariance_2D_1, "test_getVariance_1D_1")
263 3 : call Test%run(test_getVariance_2D_2, "test_getVariance_2D_2")
264 3 : call Test%run(test_getBetaCDF_SPR_1, "test_getBetaCDF_SPR_1")
265 3 : call Test%run(test_getBetaCDF_DPR_1, "test_getBetaCDF_DPR_1")
266 3 : call Test%run(test_getSNormCDF_SPR_1, "test_getSNormCDF_SPR_1")
267 3 : call Test%run(test_getSNormCDF_DPR_1, "test_getSNormCDF_DPR_1")
268 3 : call Test%run(test_getMahalSqSP_RK_1, "test_getMahalSqSP_RK_1")
269 3 : call Test%run(test_getMahalSqMP_RK_1, "test_getMahalSqMP_RK_1")
270 3 : call Test%run(test_getMahalSqSP_CK_1, "test_getMahalSqSP_CK_1")
271 3 : call Test%run(test_getMahalSqMP_CK_1, "test_getMahalSqMP_CK_1")
272 3 : call Test%run(test_getLogProbLognSP_1, "test_getLogProbLognSP_1")
273 3 : call Test%run(test_getLogProbLognMP_1, "test_getLogProbLognMP_1")
274 3 : call Test%run(test_isInsideEllipsoid_1, "test_isInsideEllipsoid_1")
275 3 : call Test%run(test_mergeMeanCovUpper_1, "test_mergeMeanCovUpper_1")
276 3 : call Test%run(test_getRandIntLecuyer_1, "test_getRandIntLecuyer_1")
277 3 : call Test%run(test_getRandRealLecuyer_1, "test_getRandRealLecuyer_1")
278 3 : call Test%run(test_getSamCovMeanTrans_1, "test_getSamCovMeanTrans_1")
279 3 : call Test%run(test_getLogProbMVNSP_RK_1, "test_getLogProbMVNSP_RK_1")
280 3 : call Test%run(test_getLogProbMVNMP_RK_1, "test_getLogProbMVNMP_RK_1")
281 3 : call Test%run(test_getLogProbMVNSP_CK_1, "test_getLogProbMVNSP_CK_1")
282 3 : call Test%run(test_getLogProbMVNMP_CK_1, "test_getLogProbMVNMP_CK_1")
283 3 : call Test%run(test_getLogProbNormSP_RK_1, "test_getLogProbNormSP_RK_1")
284 3 : call Test%run(test_getLogProbNormMP_RK_1, "test_getLogProbNormMP_RK_1")
285 3 : call Test%run(test_getLogProbNormSP_CK_1, "test_getLogProbNormSP_CK_1")
286 3 : call Test%run(test_getLogProbNormMP_CK_1, "test_getLogProbNormMP_CK_1")
287 3 : call Test%run(test_getLogProbGeoCyclic_1, "test_getLogProbGeoCyclic_1")
288 3 : call Test%run(test_getRandGammaIntShape_1, "test_getRandGammaIntShape_1")
289 3 : call Test%run(test_getRandExpWithInvMean_1, "test_getRandExpWithInvMean_1")
290 : !call Test%run(test_mergeMeanCovUpperSlow_1, "test_mergeMeanCovUpperSlow_1")
291 3 : call Test%run(test_getLogProbMixMVNSP_RK_1, "test_getLogProbMixMVNSP_RK_1")
292 3 : call Test%run(test_getLogProbMixMVNMP_RK_1, "test_getLogProbMixMVNMP_RK_1")
293 3 : call Test%run(test_getLogProbMixMVNSP_CK_1, "test_getLogProbMixMVNSP_CK_1")
294 3 : call Test%run(test_getLogProbMixMVNMP_CK_1, "test_getLogProbMixMVNMP_CK_1")
295 3 : call Test%run(test_getRandCorMatRejection_1, "test_getRandCorMatRejection_1")
296 3 : call Test%run(test_mergeMeanCovUpperDense_1, "test_mergeMeanCovUpperDense_1")
297 3 : call Test%run(test_getLogProbMixNormSP_RK_1, "test_getLogProbMixNormSP_RK_1")
298 3 : call Test%run(test_getLogProbMixNormMP_RK_1, "test_getLogProbMixNormMP_RK_1")
299 3 : call Test%run(test_getLogProbMixNormSP_CK_1, "test_getLogProbMixNormSP_CK_1")
300 3 : call Test%run(test_getLogProbMixNormMP_CK_1, "test_getLogProbMixNormMP_CK_1")
301 3 : call Test%run(test_getRandPointOnEllipsoid_1, "test_getRandPointOnEllipsoid_1")
302 3 : call Test%run(test_getSamCovUpperMeanTrans_1, "test_getSamCovUpperMeanTrans_1")
303 3 : call Test%run(test_getWeiSamCovUppMeanTrans_1, "test_getWeiSamCovUppMeanTrans_1")
304 3 : call Test%run(test_getCovMatFromCorMatUpper_1, "test_getCovMatFromCorMatUpper_1")
305 3 : call Test%run(test_getCorMatUpperFromCovMatUpper_1, "test_getCorMatUpperFromCovMatUpper_1")
306 3 : call Test%run(test_getCovMatUpperFromCorMatUpper_1, "test_getCovMatUpperFromCorMatUpper_1")
307 3 : call Test%run(test_getCovMatUpperFromCorMatLower_1, "test_getCovMatUpperFromCorMatLower_1")
308 3 : call Test%run(test_getCovMatLowerFromCorMatUpper_1, "test_getCovMatLowerFromCorMatUpper_1")
309 3 : call Test%run(test_getCovMatLowerFromCorMatLower_1, "test_getCovMatLowerFromCorMatLower_1")
310 3 : call Test%finalize()
311 3 : end subroutine test_Statistics
312 :
313 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
314 :
315 3 : function test_getMahalSqSP_RK_1() result(assertion)
316 3 : use Constants_mod, only: IK, RK
317 : implicit none
318 : integer(IK) :: i
319 : logical :: assertion
320 : integer(IK) , parameter :: nd = 3_IK
321 : real(RK) , parameter :: tolerance = 1.e-12_RK
322 : real(RK) , parameter :: mahalSq_ref = 180._RK
323 : real(RK) , parameter :: Point(nd) = [(real(i,RK),i=1,nd)]
324 : real(RK) , parameter :: MeanVec(nd) = [(real(i**2+1._RK,RK),i=1,nd)]
325 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
326 : , 0._RK, 2._RK, 0._RK &
327 : , 1._RK, 0._RK, 3._RK ], shape = shape(InvCovMat) )
328 3 : real(RK) :: mahalSq
329 3 : real(RK) :: difference
330 3 : mahalSq = getMahalSqSP_RK(nd = nd, MeanVec = MeanVec, InvCovMat = InvCovMat, Point = Point)
331 3 : difference = abs(mahalSq - mahalSq_ref) / mahalSq_ref
332 3 : assertion = difference <= tolerance
333 :
334 : ! LCOV_EXCL_START
335 : if (Test%isDebugMode .and. .not. assertion) then
336 : write(Test%outputUnit,"(*(g0,:,', '))")
337 : write(Test%outputUnit,"(*(g0,:,', '))") "mahalSq_ref ", mahalSq_ref
338 : write(Test%outputUnit,"(*(g0,:,', '))") "mahalSq ", mahalSq
339 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
340 : write(Test%outputUnit,"(*(g0,:,', '))")
341 : end if
342 : ! LCOV_EXCL_STOP
343 :
344 3 : end function test_getMahalSqSP_RK_1
345 :
346 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
347 :
348 3 : function test_getMahalSqMP_RK_1() result(assertion)
349 3 : use Constants_mod, only: IK, RK
350 : implicit none
351 : integer(IK) :: i
352 : logical :: assertion
353 : integer(IK) , parameter :: nd = 3_IK, np = 2_IK
354 : real(RK) , parameter :: tolerance = 1.e-12_RK
355 : real(RK) , parameter :: MahalSq_ref(np) = [180._RK, 36._RK]
356 : real(RK) , parameter :: Point(nd,np) = reshape([(real(i,RK),i=1,nd*np)], shape = shape(Point))
357 : real(RK) , parameter :: MeanVec(nd) = [(real(i**2+1._RK,RK),i=1,nd)]
358 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ 1._RK, 0._RK, 1._RK &
359 : , 0._RK, 2._RK, 0._RK &
360 : , 1._RK, 0._RK, 3._RK ], shape = shape(InvCovMat) )
361 : real(RK) :: MahalSq(np)
362 : real(RK) :: Difference(np)
363 3 : MahalSq = getMahalSqMP_RK(nd = nd, np = np, MeanVec = MeanVec, InvCovMat = InvCovMat, Point = Point)
364 9 : Difference = abs(MahalSq - MahalSq_ref) / MahalSq_ref
365 9 : assertion = all(Difference <= tolerance)
366 :
367 : ! LCOV_EXCL_START
368 : if (Test%isDebugMode .and. .not. assertion) then
369 : write(Test%outputUnit,"(*(g0,:,', '))")
370 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq_ref ", MahalSq_ref
371 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq ", MahalSq
372 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", Difference
373 : write(Test%outputUnit,"(*(g0,:,', '))")
374 : end if
375 : ! LCOV_EXCL_STOP
376 :
377 3 : end function test_getMahalSqMP_RK_1
378 :
379 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
380 :
381 3 : function test_getMahalSqSP_CK_1() result(assertion)
382 :
383 3 : use Constants_mod, only: IK, RK, CK
384 : implicit none
385 : integer(IK) :: i
386 : logical :: assertion
387 : integer(IK) , parameter :: nd = 3_IK
388 : real(RK) , parameter :: tolerance = 1.e-12_RK
389 : complex(CK) , parameter :: mahalSq_ref = 180._CK
390 : complex(CK) , parameter :: Point(nd) = [(cmplx(i,0.,kind=RK),i=1,nd)]
391 : complex(CK) , parameter :: MeanVec(nd) = [(cmplx(i**2+1._RK,0.,kind=RK),i=1,nd)]
392 : complex(CK) , parameter :: InvCovMat(nd,nd) = cmplx(reshape([ 1._RK, 0._RK, 1._RK &
393 : , 0._RK, 2._RK, 0._RK &
394 : , 1._RK, 0._RK, 3._RK ], shape = shape(InvCovMat) ), kind = RK )
395 : complex(CK) :: mahalSq
396 3 : real(RK) :: difference
397 3 : mahalSq = getMahalSqSP_CK(nd = nd, MeanVec = MeanVec, InvCovMat = InvCovMat, Point = Point)
398 3 : difference = abs(real(mahalSq - mahalSq_ref,RK)) / real(mahalSq_ref,RK)
399 3 : assertion = difference <= tolerance
400 :
401 : ! LCOV_EXCL_START
402 : if (Test%isDebugMode .and. .not. assertion) then
403 : write(Test%outputUnit,"(*(g0,:,', '))")
404 : write(Test%outputUnit,"(*(g0,:,', '))") "mahalSq_ref ", mahalSq_ref
405 : write(Test%outputUnit,"(*(g0,:,', '))") "mahalSq ", mahalSq
406 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
407 : write(Test%outputUnit,"(*(g0,:,', '))")
408 : end if
409 : ! LCOV_EXCL_STOP
410 :
411 3 : end function test_getMahalSqSP_CK_1
412 :
413 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
414 :
415 3 : function test_getMahalSqMP_CK_1() result(assertion)
416 3 : use Constants_mod, only: IK, RK, CK
417 : implicit none
418 : integer(IK) :: i
419 : logical :: assertion
420 : real(RK) , parameter :: tolerance = 1.e-12_RK
421 : integer(IK) , parameter :: nd = 3_IK, np = 2_IK
422 : complex(CK) , parameter :: MahalSq_ref(np) = [180._CK, 36._CK]
423 : complex(CK) , parameter :: Point(nd,np) = reshape([(cmplx(i,0.,kind=RK),i=1,nd*np)], shape = shape(Point))
424 : complex(CK) , parameter :: MeanVec(nd) = [(cmplx(i**2+1._RK,0.,kind=RK),i=1,nd)]
425 : complex(CK) , parameter :: InvCovMat(nd,nd) = cmplx(reshape([ 1._RK, 0._RK, 1._RK &
426 : , 0._RK, 2._RK, 0._RK &
427 : , 1._RK, 0._RK, 3._RK ], shape = shape(InvCovMat) ), kind=RK )
428 :
429 : complex(CK) :: MahalSq(np)
430 : real(RK) :: Difference(np)
431 3 : MahalSq = getMahalSqMP_CK(nd = nd, np = np, MeanVec = MeanVec, InvCovMat = InvCovMat, Point = Point)
432 9 : Difference = abs(real(MahalSq - MahalSq_ref,RK) / real(MahalSq_ref,RK))
433 9 : assertion = all(Difference <= tolerance)
434 :
435 : ! LCOV_EXCL_START
436 : if (Test%isDebugMode .and. .not. assertion) then
437 : write(Test%outputUnit,"(*(g0,:,', '))")
438 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq_ref ", MahalSq_ref
439 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq ", MahalSq
440 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", Difference
441 : write(Test%outputUnit,"(*(g0,:,', '))")
442 : end if
443 : ! LCOV_EXCL_STOP
444 :
445 3 : end function test_getMahalSqMP_CK_1
446 :
447 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
448 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
449 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
450 :
451 3 : function test_getLogProbNormSP_RK_1() result(assertion)
452 :
453 3 : use Constants_mod, only: IK, RK
454 : implicit none
455 :
456 : logical :: assertion
457 : real(RK) , parameter :: mean = 3._RK
458 : real(RK) , parameter :: point = 2._RK
459 : real(RK) , parameter :: inverseVariance = 1._RK / 16._RK
460 : real(RK) , parameter :: logProbNorm_ref = -2.336482894324563_RK
461 : real(RK) , parameter :: logSqrtInverseVariance = log(sqrt(inverseVariance))
462 : real(RK) , parameter :: tolerance = 1.e-12_RK
463 3 : real(RK) :: logProbNorm
464 3 : real(RK) :: difference
465 :
466 : logProbNorm = getLogProbNormSP_RK ( mean = mean &
467 : , inverseVariance = inverseVariance &
468 : , logSqrtInverseVariance = logSqrtInverseVariance &
469 : , point = point &
470 3 : )
471 :
472 3 : difference = abs( (logProbNorm - logProbNorm_ref) / logProbNorm_ref )
473 3 : assertion = difference <= tolerance
474 :
475 : ! LCOV_EXCL_START
476 : if (Test%isDebugMode .and. .not. assertion) then
477 : write(Test%outputUnit,"(*(g0,:,', '))")
478 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", logProbNorm_ref
479 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", logProbNorm
480 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
481 : write(Test%outputUnit,"(*(g0,:,', '))")
482 : end if
483 : ! LCOV_EXCL_STOP
484 :
485 3 : end function test_getLogProbNormSP_RK_1
486 :
487 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
488 :
489 3 : function test_getLogProbNormMP_RK_1() result(assertion)
490 :
491 3 : use Constants_mod, only: IK, RK
492 : implicit none
493 :
494 : logical :: assertion
495 : integer(IK) , parameter :: np = 2_IK
496 : real(RK) , parameter :: mean = 3._RK
497 : real(RK) , parameter :: point(np) = [2._RK, 3._RK]
498 : real(RK) , parameter :: inverseVariance = 1._RK / 16._RK
499 : real(RK) , parameter :: logProbNorm_ref(np) = [ -2.336482894324563_RK, -2.305232894324563_RK ]
500 : real(RK) , parameter :: logSqrtInverseVariance = log(sqrt(inverseVariance))
501 : real(RK) , parameter :: tolerance = 1.e-12_RK
502 : real(RK) :: logProbNorm(np)
503 : real(RK) :: difference(np)
504 :
505 : logProbNorm = getLogProbNormMP_RK ( np = np &
506 : , mean = mean &
507 : , inverseVariance = inverseVariance &
508 : , logSqrtInverseVariance = logSqrtInverseVariance &
509 : , point = point &
510 3 : )
511 :
512 9 : difference = abs( (logProbNorm - logProbNorm_ref) / logProbNorm_ref )
513 9 : assertion = all(difference <= tolerance)
514 :
515 : ! LCOV_EXCL_START
516 : if (Test%isDebugMode .and. .not. assertion) then
517 : write(Test%outputUnit,"(*(g0,:,', '))")
518 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", logProbNorm_ref
519 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", logProbNorm
520 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
521 : write(Test%outputUnit,"(*(g0,:,', '))")
522 : end if
523 : ! LCOV_EXCL_STOP
524 :
525 3 : end function test_getLogProbNormMP_RK_1
526 :
527 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
528 :
529 3 : function test_getLogProbMVNSP_RK_1() result(assertion)
530 :
531 3 : use Constants_mod, only: IK, RK
532 : implicit none
533 :
534 : logical :: assertion
535 : integer(IK) :: i
536 : integer(IK) , parameter :: nd = 3_IK
537 : real(RK) , parameter :: tolerance = 1.e-12_RK
538 : real(RK) , parameter :: Point(nd) = [(real(i,RK),i=1,nd)]
539 : real(RK) , parameter :: MeanVec(nd) = [(real(i**2+1._RK,RK),i=1,nd)]
540 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
541 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
542 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
543 : , shape = shape(InvCovMat) )
544 : real(RK) , parameter :: logSqrtDetInvCovMat = -0.693147180559945_RK
545 : real(RK) , parameter :: logProbNorm_ref = -15.19996278017396_RK
546 3 : real(RK) :: logProbNorm
547 3 : real(RK) :: difference
548 :
549 : logProbNorm = getLogProbMVNSP_RK( nd = nd &
550 : , MeanVec = MeanVec &
551 : , InvCovMat = InvCovMat &
552 : , logSqrtDetInvCovMat = logSqrtDetInvCovMat &
553 : , point = point &
554 3 : )
555 :
556 3 : difference = abs( (logProbNorm - logProbNorm_ref) / logProbNorm_ref )
557 3 : assertion = difference <= tolerance
558 :
559 : ! LCOV_EXCL_START
560 : if (Test%isDebugMode .and. .not. assertion) then
561 : write(Test%outputUnit,"(*(g0,:,', '))")
562 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", logProbNorm_ref
563 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", logProbNorm
564 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
565 : write(Test%outputUnit,"(*(g0,:,', '))")
566 : end if
567 : ! LCOV_EXCL_STOP
568 :
569 3 : end function test_getLogProbMVNSP_RK_1
570 :
571 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
572 :
573 3 : function test_getLogProbMVNMP_RK_1() result(assertion)
574 :
575 3 : use Constants_mod, only: IK, RK
576 : implicit none
577 :
578 : logical :: assertion
579 : integer(IK) :: i
580 : integer(IK) , parameter :: np = 2_IK
581 : integer(IK) , parameter :: nd = 3_IK
582 : real(RK) , parameter :: tolerance = 1.e-12_RK
583 : real(RK) , parameter :: Point(nd,np) = reshape( [(real(i,RK),i=1,nd*np)], shape = shape(Point) )
584 : real(RK) , parameter :: MeanVec(nd) = [(real(i**2+1._RK,RK),i=1,nd)]
585 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
586 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
587 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
588 : , shape = shape(InvCovMat) )
589 : real(RK) , parameter :: logSqrtDetInvCovMat = -0.693147180559945_RK
590 : real(RK) , parameter :: logProbNorm_ref(np) = [ -15.19996278017396_RK, -14.44996278017396_RK ]
591 : real(RK) :: logProbNorm(np)
592 : real(RK) :: difference(np)
593 :
594 : logProbNorm = getLogProbMVNMP_RK( nd = nd &
595 : , np = np &
596 : , MeanVec = MeanVec &
597 : , InvCovMat = InvCovMat &
598 : , logSqrtDetInvCovMat = logSqrtDetInvCovMat &
599 : , point = point &
600 3 : )
601 :
602 9 : difference = abs( (logProbNorm - logProbNorm_ref) / logProbNorm_ref )
603 9 : assertion = all(difference <= tolerance)
604 :
605 : ! LCOV_EXCL_START
606 : if (Test%isDebugMode .and. .not. assertion) then
607 : write(Test%outputUnit,"(*(g0,:,', '))")
608 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", logProbNorm_ref
609 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", logProbNorm
610 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
611 : write(Test%outputUnit,"(*(g0,:,', '))")
612 : end if
613 : ! LCOV_EXCL_STOP
614 :
615 3 : end function test_getLogProbMVNMP_RK_1
616 :
617 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
618 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
619 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
620 :
621 3 : function test_getLogProbNormSP_CK_1() result(assertion)
622 :
623 3 : use Constants_mod, only: IK, RK, CK
624 : implicit none
625 :
626 : logical :: assertion
627 : real(RK) , parameter :: tolerance = 1.e-12_RK
628 : complex(CK) , parameter :: mean = 3._RK
629 : complex(CK) , parameter :: point = cmplx(2._RK,kind=RK)
630 : complex(CK) , parameter :: inverseVariance = cmplx(1._RK / 16._RK,kind=RK)
631 : complex(CK) , parameter :: logProbNorm_ref = cmplx(-2.336482894324563_RK,kind=RK)
632 : complex(CK) , parameter :: logSqrtInverseVariance = cmplx(log(sqrt(inverseVariance)),kind=RK)
633 : complex(CK) :: logProbNorm
634 3 : real(RK) :: difference
635 :
636 : logProbNorm = getLogProbNormSP_CK ( mean = mean &
637 : , inverseVariance = inverseVariance &
638 : , logSqrtInverseVariance = logSqrtInverseVariance &
639 : , point = point &
640 3 : )
641 :
642 3 : difference = abs( real(logProbNorm - logProbNorm_ref,kind=RK) / real(logProbNorm_ref,kind=RK) )
643 3 : assertion = difference <= tolerance
644 :
645 : ! LCOV_EXCL_START
646 : if (Test%isDebugMode .and. .not. assertion) then
647 : write(Test%outputUnit,"(*(g0,:,', '))")
648 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", real(logProbNorm_ref, kind = RK)
649 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", real(logProbNorm, kind = RK)
650 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
651 : write(Test%outputUnit,"(*(g0,:,', '))")
652 : end if
653 : ! LCOV_EXCL_STOP
654 :
655 3 : end function test_getLogProbNormSP_CK_1
656 :
657 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
658 :
659 3 : function test_getLogProbNormMP_CK_1() result(assertion)
660 :
661 3 : use Constants_mod, only: IK, RK, CK
662 : implicit none
663 :
664 : logical :: assertion
665 : real(RK) , parameter :: tolerance = 1.e-12_RK
666 : integer(IK) , parameter :: np = 2_IK
667 : complex(CK) , parameter :: mean = 3._RK
668 : complex(CK) , parameter :: point(np) = cmplx([2._RK, 3._RK], kind=RK)
669 : complex(CK) , parameter :: inverseVariance = cmplx(1._RK / 16._RK, kind=RK)
670 : complex(CK) , parameter :: logProbNorm_ref(np) = cmplx([ -2.336482894324563_RK, -2.305232894324563_RK ], kind=RK)
671 : complex(CK) , parameter :: logSqrtInverseVariance = cmplx(log(sqrt(inverseVariance)), kind=RK)
672 : complex(CK) :: logProbNorm(np)
673 : real(RK) :: difference(np)
674 :
675 : logProbNorm = getLogProbNormMP_CK ( np = np &
676 : , mean = mean &
677 : , inverseVariance = inverseVariance &
678 : , logSqrtInverseVariance = logSqrtInverseVariance &
679 : , point = point &
680 3 : )
681 :
682 9 : difference = abs( (real(logProbNorm, kind=RK) - real(logProbNorm_ref, kind=RK)) / real(logProbNorm_ref, kind=RK) )
683 9 : assertion = all(difference <= tolerance)
684 :
685 : ! LCOV_EXCL_START
686 : if (Test%isDebugMode .and. .not. assertion) then
687 : write(Test%outputUnit,"(*(g0,:,', '))")
688 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", real(logProbNorm_ref, kind = RK)
689 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", real(logProbNorm, kind = RK)
690 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
691 : write(Test%outputUnit,"(*(g0,:,', '))")
692 : end if
693 : ! LCOV_EXCL_STOP
694 :
695 3 : end function test_getLogProbNormMP_CK_1
696 :
697 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
698 :
699 3 : function test_getLogProbMVNSP_CK_1() result(assertion)
700 :
701 3 : use Constants_mod, only: IK, RK, CK
702 : implicit none
703 :
704 : logical :: assertion
705 : integer(IK) :: i
706 : integer(IK) , parameter :: nd = 3_IK
707 : real(RK) , parameter :: tolerance = 1.e-12_RK
708 : complex(CK) , parameter :: Point(nd) = cmplx([(real(i,RK),i=1,nd)], kind=RK)
709 : complex(CK) , parameter :: MeanVec(nd) = cmplx([(real(i**2+1._RK,RK),i=1,nd)], kind=RK)
710 : complex(CK) , parameter :: InvCovMat(nd,nd) = cmplx( reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
711 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
712 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
713 : , shape = shape(InvCovMat) ), kind=RK)
714 : complex(CK) , parameter :: logSqrtDetInvCovMat = cmplx(-0.693147180559945_RK, kind=RK)
715 : complex(CK) , parameter :: logProbNorm_ref = cmplx(-15.19996278017396_RK, kind=RK)
716 : complex(CK) :: logProbNorm
717 3 : real(RK) :: difference
718 :
719 : logProbNorm = getLogProbMVNSP_CK( nd = nd &
720 : , MeanVec = MeanVec &
721 : , InvCovMat = InvCovMat &
722 : , logSqrtDetInvCovMat = logSqrtDetInvCovMat &
723 : , point = point &
724 6 : )
725 :
726 3 : difference = abs( (real(logProbNorm, kind=RK) - real(logProbNorm_ref, kind=RK)) / real(logProbNorm_ref, kind=RK) )
727 3 : assertion = difference <= tolerance
728 :
729 : ! LCOV_EXCL_START
730 : if (Test%isDebugMode .and. .not. assertion) then
731 : write(Test%outputUnit,"(*(g0,:,', '))")
732 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", real(logProbNorm_ref, kind = RK)
733 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", real(logProbNorm, kind = RK)
734 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
735 : write(Test%outputUnit,"(*(g0,:,', '))")
736 : end if
737 : ! LCOV_EXCL_STOP
738 :
739 3 : end function test_getLogProbMVNSP_CK_1
740 :
741 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
742 :
743 3 : function test_getLogProbMVNMP_CK_1() result(assertion)
744 :
745 3 : use Constants_mod, only: IK, RK, CK
746 : implicit none
747 :
748 : logical :: assertion
749 : integer(IK) :: i
750 : real(RK) , parameter :: tolerance = 1.e-12_RK
751 : integer(IK) , parameter :: np = 2_IK
752 : integer(IK) , parameter :: nd = 3_IK
753 : complex(CK) , parameter :: Point(nd,np) = reshape( cmplx([(real(i,RK),i=1,nd*np)], kind=RK), shape = shape(Point) )
754 : complex(CK) , parameter :: MeanVec(nd) = cmplx([(real(i**2+1._RK,RK),i=1,nd)], kind=RK)
755 : complex(CK) , parameter :: InvCovMat(nd,nd) = cmplx( reshape( [ 1.500000000000000_RK, 0.000000000000000_RK, -0.50000000000000_RK &
756 : , 0.000000000000000_RK, 0.500000000000000_RK, 0.000000000000000_RK &
757 : , -0.50000000000000_RK, 0.000000000000000_RK, 0.500000000000000_RK ] &
758 : , shape = shape(InvCovMat) ), kind=RK)
759 : complex(CK) , parameter :: logSqrtDetInvCovMat = cmplx(-0.693147180559945_RK, kind=RK)
760 : complex(CK) , parameter :: logProbNorm_ref(np) = cmplx([ -15.19996278017396_RK, -14.44996278017396_RK ], kind=RK)
761 : complex(CK) :: logProbNorm(np)
762 : real(RK) :: difference(np)
763 :
764 : logProbNorm = getLogProbMVNMP_CK( nd = nd &
765 : , np = np &
766 : , MeanVec = MeanVec &
767 : , InvCovMat = InvCovMat &
768 : , logSqrtDetInvCovMat = logSqrtDetInvCovMat &
769 : , point = point &
770 3 : )
771 :
772 9 : difference = abs( (real(logProbNorm, kind=RK) - real(logProbNorm_ref, kind=RK)) / real(logProbNorm_ref, kind=RK) )
773 9 : assertion = all(difference <= tolerance)
774 :
775 : ! LCOV_EXCL_START
776 : if (Test%isDebugMode .and. .not. assertion) then
777 : write(Test%outputUnit,"(*(g0,:,', '))")
778 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm_ref ", real(logProbNorm_ref, kind = RK)
779 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbNorm ", real(logProbNorm, kind = RK)
780 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
781 : write(Test%outputUnit,"(*(g0,:,', '))")
782 : end if
783 : ! LCOV_EXCL_STOP
784 :
785 3 : end function test_getLogProbMVNMP_CK_1
786 :
787 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
788 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
789 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
790 :
791 3 : function test_getLogProbMixNormSP_RK_1() result(assertion)
792 :
793 3 : use Constants_mod, only: IK, RK
794 : implicit none
795 :
796 : logical :: assertion
797 : integer(IK) , parameter :: nmode = 2_IK
798 : real(RK) , parameter :: point = 2._RK
799 : real(RK) , parameter :: tolerance = 1.e-12_RK
800 : real(RK) , parameter :: MeanVec(nmode) = [ 0.25_RK, 0.75_RK ]
801 : real(RK) , parameter :: InvCovMat(nmode) = [ 1._RK / 16._RK, 1._RK / 32._RK ]
802 : real(RK) , parameter :: LogAmplitude(nmode) = [ 3._RK, 4._RK ]
803 : real(RK) , parameter :: LogSqrtDetInvCovMat(nmode) = log(sqrt(InvCovMat))
804 : real(RK) , parameter :: logProbMixNorm_ref = 1.718832134253714_RK
805 3 : real(RK) :: logProbMixNorm
806 3 : real(RK) :: difference
807 :
808 : logProbMixNorm = getLogProbMixNorm ( nmode = nmode &
809 : , LogAmplitude = LogAmplitude &
810 : , MeanVec = MeanVec &
811 : , InvCovMat = InvCovMat &
812 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
813 : , point = point &
814 3 : )
815 :
816 3 : difference = abs( (logProbMixNorm - logProbMixNorm_ref) / logProbMixNorm_ref )
817 3 : assertion = difference <= tolerance
818 :
819 : ! LCOV_EXCL_START
820 : if (Test%isDebugMode .and. .not. assertion) then
821 : write(Test%outputUnit,"(*(g0,:,', '))")
822 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm_ref ", logProbMixNorm_ref
823 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm ", logProbMixNorm
824 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
825 : write(Test%outputUnit,"(*(g0,:,', '))")
826 : end if
827 : ! LCOV_EXCL_STOP
828 :
829 3 : end function test_getLogProbMixNormSP_RK_1
830 :
831 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
832 :
833 3 : function test_getLogProbMixNormMP_RK_1() result(assertion)
834 :
835 3 : use Constants_mod, only: IK, RK
836 : implicit none
837 :
838 : logical :: assertion
839 : integer(IK) , parameter :: np = 2_IK
840 : integer(IK) , parameter :: nmode = 2_IK
841 : real(RK) , parameter :: point(np) = [ 2._RK, 3._RK ]
842 : real(RK) , parameter :: tolerance = 1.e-12_RK
843 : real(RK) , parameter :: MeanVec(nmode) = [ 0.25_RK, 0.75_RK ]
844 : real(RK) , parameter :: InvCovMat(nmode) = [ 1._RK / 16._RK, 1._RK / 32._RK ]
845 : real(RK) , parameter :: LogAmplitude(nmode) = [ 3._RK, 4._RK ]
846 : real(RK) , parameter :: LogSqrtDetInvCovMat(nmode) = log(sqrt(InvCovMat))
847 : real(RK) , parameter :: logProbMixNorm_ref(np) = [ 1.718832134253714_RK, 1.636902047052812_RK ]
848 : real(RK) :: logProbMixNorm(np)
849 : real(RK) :: difference(np)
850 :
851 : logProbMixNorm = getLogProbMixNorm ( nmode = nmode &
852 : , np = np &
853 : , LogAmplitude = LogAmplitude &
854 : , MeanVec = MeanVec &
855 : , InvCovMat = InvCovMat &
856 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
857 : , point = point &
858 3 : )
859 :
860 9 : difference = abs( (logProbMixNorm - logProbMixNorm_ref) / logProbMixNorm_ref )
861 9 : assertion = all(difference <= tolerance)
862 :
863 : ! LCOV_EXCL_START
864 : if (Test%isDebugMode .and. .not. assertion) then
865 : write(Test%outputUnit,"(*(g0,:,', '))")
866 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm_ref ", logProbMixNorm_ref
867 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm ", logProbMixNorm
868 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
869 : write(Test%outputUnit,"(*(g0,:,', '))")
870 : end if
871 : ! LCOV_EXCL_STOP
872 :
873 3 : end function test_getLogProbMixNormMP_RK_1
874 :
875 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
876 :
877 3 : function test_getLogProbMixMVNSP_RK_1() result(assertion)
878 :
879 3 : use Constants_mod, only: IK, RK
880 : implicit none
881 :
882 : logical :: assertion
883 : integer(IK) :: i
884 : integer(IK) , parameter :: nd = 3_IK
885 : integer(IK) , parameter :: nmode = 2_IK
886 : real(RK) , parameter :: tolerance = 1.e-12_RK
887 : real(RK) , parameter :: Point(nd) = [(real(i,RK),i=1,nd)]
888 : real(RK) , parameter :: MeanVec(nd,nmode) = reshape([(real(i**2+1._RK,RK),i=1,nd*nmode)], shape = shape(MeanVec))
889 : real(RK) , parameter :: InvCovMat(nd,nd,nmode) = reshape([ 1._RK, 0._RK, 1._RK &
890 : , 0._RK, 2._RK, 0._RK &
891 : , 1._RK, 0._RK, 3._RK &
892 : , 2._RK, 0._RK, 0._RK &
893 : , 0._RK, 2._RK, 0._RK &
894 : , 0._RK, 0._RK, 2._RK &
895 : ], shape = shape(InvCovMat) )
896 : real(RK) , parameter :: LogAmplitude(nmode) = [ 3._RK, 4._RK ]
897 : real(RK) , parameter :: LogSqrtDetInvCovMat(nmode) = [-0.693147180559945_RK, -1.039720770839918_RK]
898 : real(RK) , parameter :: logProbMixMVN_ref = -90.44996278017396_RK
899 3 : real(RK) :: logProbMixMVN
900 3 : real(RK) :: difference
901 :
902 : logProbMixMVN = getlogProbMixMVN( nmode = nmode &
903 : , nd = nd &
904 : , LogAmplitude = LogAmplitude &
905 : , MeanVec = MeanVec &
906 : , InvCovMat = InvCovMat &
907 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
908 : , point = point &
909 3 : )
910 :
911 3 : difference = abs( (logProbMixMVN - logProbMixMVN_ref) / logProbMixMVN_ref )
912 3 : assertion = difference <= tolerance
913 :
914 : ! LCOV_EXCL_START
915 : if (Test%isDebugMode .and. .not. assertion) then
916 : write(Test%outputUnit,"(*(g0,:,', '))")
917 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN_ref ", logProbMixMVN_ref
918 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN ", logProbMixMVN
919 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
920 : write(Test%outputUnit,"(*(g0,:,', '))")
921 : end if
922 : ! LCOV_EXCL_STOP
923 :
924 3 : end function test_getLogProbMixMVNSP_RK_1
925 :
926 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
927 :
928 3 : function test_getLogProbMixMVNMP_RK_1() result(assertion)
929 :
930 3 : use Constants_mod, only: IK, RK
931 : implicit none
932 :
933 : logical :: assertion
934 : integer(IK) :: i
935 : integer(IK) , parameter :: nd = 3_IK
936 : integer(IK) , parameter :: np = 2_IK
937 : integer(IK) , parameter :: nmode = 2_IK
938 : real(RK) , parameter :: tolerance = 1.e-12_RK
939 : real(RK) , parameter :: Point(nd,np) = reshape([(real(i,RK),i=1,nd*np)], shape = shape(Point))
940 : real(RK) , parameter :: MeanVec(nd,nmode) = reshape([(real(i**2+1._RK,RK),i=1,nd*nmode)], shape = shape(MeanVec))
941 : real(RK) , parameter :: InvCovMat(nd,nd,nmode) = reshape([ 1._RK, 0._RK, 1._RK &
942 : , 0._RK, 2._RK, 0._RK &
943 : , 1._RK, 0._RK, 3._RK &
944 : , 2._RK, 0._RK, 0._RK &
945 : , 0._RK, 2._RK, 0._RK &
946 : , 0._RK, 0._RK, 2._RK &
947 : ], shape = shape(InvCovMat) )
948 : real(RK) , parameter :: LogAmplitude(nmode) = [ 3._RK, 4._RK ]
949 : real(RK) , parameter :: LogSqrtDetInvCovMat(nmode) = [-0.693147180559945_RK, -1.039720770839918_RK]
950 : real(RK) , parameter :: logProbMixMVN_ref(np) = [ -90.44996278017396_RK, -18.44996278017396_RK ]
951 : real(RK) :: logProbMixMVN(np)
952 : real(RK) :: difference(np)
953 :
954 : logProbMixMVN = getlogProbMixMVN( nmode = nmode &
955 : , nd = nd &
956 : , np = np &
957 : , LogAmplitude = LogAmplitude &
958 : , MeanVec = MeanVec &
959 : , InvCovMat = InvCovMat &
960 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
961 : , point = point &
962 3 : )
963 :
964 9 : difference = abs( (logProbMixMVN - logProbMixMVN_ref) / logProbMixMVN_ref )
965 9 : assertion = all(difference <= tolerance)
966 :
967 : ! LCOV_EXCL_START
968 : if (Test%isDebugMode .and. .not. assertion) then
969 : write(Test%outputUnit,"(*(g0,:,', '))")
970 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN_ref ", logProbMixMVN_ref
971 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN ", logProbMixMVN
972 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
973 : write(Test%outputUnit,"(*(g0,:,', '))")
974 : end if
975 : ! LCOV_EXCL_STOP
976 :
977 3 : end function test_getLogProbMixMVNMP_RK_1
978 :
979 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
980 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
981 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
982 :
983 3 : function test_getLogProbMixNormSP_CK_1() result(assertion)
984 :
985 3 : use Constants_mod, only: IK, RK, CK
986 : implicit none
987 :
988 : logical :: assertion
989 : real(RK) , parameter :: tolerance = 1.e-12_RK
990 : integer(IK) , parameter :: nmode = 2_IK
991 : complex(CK) , parameter :: point = 2._RK
992 : complex(CK) , parameter :: MeanVec(nmode) = cmplx([ 0.25_RK, 0.75_RK ], kind = RK)
993 : complex(CK) , parameter :: InvCovMat(nmode) = cmplx([ 1._RK / 16._RK, 1._RK / 32._RK ], kind = RK)
994 : complex(CK) , parameter :: LogAmplitude(nmode) = cmplx([ 3._RK, 4._RK ], kind = RK)
995 : complex(CK) , parameter :: LogSqrtDetInvCovMat(nmode) = cmplx(log(sqrt(InvCovMat)), kind = RK)
996 : complex(CK) , parameter :: logProbMixNorm_ref = cmplx(1.718832134253714_RK, kind = RK)
997 : complex(CK) :: logProbMixNorm
998 3 : real(RK) :: difference
999 :
1000 : logProbMixNorm = getLogProbMixNorm ( nmode = nmode &
1001 : , LogAmplitude = LogAmplitude &
1002 : , MeanVec = MeanVec &
1003 : , InvCovMat = InvCovMat &
1004 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
1005 : , point = point &
1006 6 : )
1007 :
1008 3 : difference = abs( (real(logProbMixNorm, kind = RK) - real(logProbMixNorm_ref, kind = RK)) / real(logProbMixNorm_ref, kind = RK) )
1009 3 : assertion = difference <= tolerance
1010 :
1011 : ! LCOV_EXCL_START
1012 : if (Test%isDebugMode .and. .not. assertion) then
1013 : write(Test%outputUnit,"(*(g0,:,', '))")
1014 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm_ref ", real(logProbMixNorm_ref, kind = RK)
1015 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm ", real(logProbMixNorm, kind = RK)
1016 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
1017 : write(Test%outputUnit,"(*(g0,:,', '))")
1018 : end if
1019 : ! LCOV_EXCL_STOP
1020 :
1021 3 : end function test_getLogProbMixNormSP_CK_1
1022 :
1023 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1024 :
1025 3 : function test_getLogProbMixNormMP_CK_1() result(assertion)
1026 :
1027 3 : use Constants_mod, only: IK, RK, CK
1028 : implicit none
1029 :
1030 : logical :: assertion
1031 : real(RK) , parameter :: tolerance = 1.e-12_RK
1032 : integer(IK) , parameter :: np = 2_IK
1033 : integer(IK) , parameter :: nmode = 2_IK
1034 : complex(CK) , parameter :: point(np) = cmplx([ 2._RK, 3._RK ], kind = RK)
1035 : complex(CK) , parameter :: MeanVec(nmode) = cmplx([ 0.25_RK, 0.75_RK ], kind = RK)
1036 : complex(CK) , parameter :: InvCovMat(nmode) = cmplx([ 1._RK / 16._RK, 1._RK / 32._RK ], kind = RK)
1037 : complex(CK) , parameter :: LogAmplitude(nmode) = cmplx([ 3._RK, 4._RK ], kind = RK)
1038 : complex(CK) , parameter :: LogSqrtDetInvCovMat(nmode) = cmplx(log(sqrt(InvCovMat)), kind = RK)
1039 : complex(CK) , parameter :: logProbMixNorm_ref(np) = cmplx([ 1.718832134253714_RK, 1.636902047052812_RK ], kind = RK)
1040 : complex(CK) :: logProbMixNorm(np)
1041 : real(RK) :: difference(np)
1042 :
1043 : logProbMixNorm = getLogProbMixNorm ( nmode = nmode &
1044 : , np = np &
1045 : , LogAmplitude = LogAmplitude &
1046 : , MeanVec = MeanVec &
1047 : , InvCovMat = InvCovMat &
1048 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
1049 : , point = point &
1050 3 : )
1051 :
1052 9 : difference = abs( (real(logProbMixNorm, kind = RK) - real(logProbMixNorm_ref, kind = RK)) / real(logProbMixNorm_ref, kind = RK) )
1053 9 : assertion = all(difference <= tolerance)
1054 :
1055 : ! LCOV_EXCL_START
1056 : if (Test%isDebugMode .and. .not. assertion) then
1057 : write(Test%outputUnit,"(*(g0,:,', '))")
1058 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm_ref ", real(logProbMixNorm_ref, kind = RK)
1059 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixNorm ", real(logProbMixNorm, kind = RK)
1060 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
1061 : write(Test%outputUnit,"(*(g0,:,', '))")
1062 : end if
1063 : ! LCOV_EXCL_STOP
1064 :
1065 3 : end function test_getLogProbMixNormMP_CK_1
1066 :
1067 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1068 :
1069 3 : function test_getLogProbMixMVNSP_CK_1() result(assertion)
1070 :
1071 3 : use Constants_mod, only: IK, RK, CK
1072 : implicit none
1073 :
1074 : logical :: assertion
1075 : real(RK) , parameter :: tolerance = 1.e-12_RK
1076 : integer(IK) :: i
1077 : integer(IK) , parameter :: nd = 3_IK
1078 : integer(IK) , parameter :: nmode = 2_IK
1079 : complex(CK) , parameter :: Point(nd) = cmplx([(real(i,RK),i=1,nd)], kind = RK)
1080 : complex(CK) , parameter :: MeanVec(nd,nmode) = cmplx(reshape([(real(i**2+1._RK,RK),i=1,nd*nmode)], shape = shape(MeanVec)), kind = RK)
1081 : complex(CK) , parameter :: InvCovMat(nd,nd,nmode) = cmplx( reshape( [ 1._RK, 0._RK, 1._RK &
1082 : , 0._RK, 2._RK, 0._RK &
1083 : , 1._RK, 0._RK, 3._RK &
1084 : , 2._RK, 0._RK, 0._RK &
1085 : , 0._RK, 2._RK, 0._RK &
1086 : , 0._RK, 0._RK, 2._RK &
1087 : ], shape = shape(InvCovMat) ), kind = RK)
1088 : complex(CK) , parameter :: LogAmplitude(nmode) = cmplx([ 3._RK, 4._RK ], kind = RK)
1089 : complex(CK) , parameter :: LogSqrtDetInvCovMat(nmode) = cmplx([-0.693147180559945_RK, -1.039720770839918_RK], kind = RK)
1090 : complex(CK) , parameter :: logProbMixMVN_ref = cmplx(-90.44996278017396_RK, kind = RK)
1091 : complex(CK) :: logProbMixMVN
1092 3 : real(RK) :: difference
1093 :
1094 : logProbMixMVN = getlogProbMixMVN( nmode = nmode &
1095 : , nd = nd &
1096 : , LogAmplitude = LogAmplitude &
1097 : , MeanVec = MeanVec &
1098 : , InvCovMat = InvCovMat &
1099 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
1100 : , point = point &
1101 6 : )
1102 :
1103 3 : difference = abs( (real(logProbMixMVN, kind = RK) - real(logProbMixMVN_ref, kind = RK)) / real(logProbMixMVN_ref, kind = RK) )
1104 3 : assertion = difference <= tolerance
1105 :
1106 : ! LCOV_EXCL_START
1107 : if (Test%isDebugMode .and. .not. assertion) then
1108 : write(Test%outputUnit,"(*(g0,:,', '))")
1109 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN_ref ", real(logProbMixMVN_ref, kind = RK)
1110 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN ", real(logProbMixMVN, kind = RK)
1111 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
1112 : write(Test%outputUnit,"(*(g0,:,', '))")
1113 : end if
1114 : ! LCOV_EXCL_STOP
1115 :
1116 3 : end function test_getLogProbMixMVNSP_CK_1
1117 :
1118 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1119 :
1120 3 : function test_getLogProbMixMVNMP_CK_1() result(assertion)
1121 :
1122 3 : use Constants_mod, only: IK, RK, CK
1123 : implicit none
1124 :
1125 : logical :: assertion
1126 : integer(IK) :: i
1127 : real(RK) , parameter :: tolerance = 1.e-12_RK
1128 : integer(IK) , parameter :: nd = 3_IK
1129 : integer(IK) , parameter :: np = 2_IK
1130 : integer(IK) , parameter :: nmode = 2_IK
1131 : complex(CK) , parameter :: Point(nd,np) = cmplx(reshape([(real(i,RK),i=1,nd*np)], shape = shape(Point)), kind = RK)
1132 : complex(CK) , parameter :: MeanVec(nd,nmode) = cmplx(reshape([(real(i**2+1._RK,RK),i=1,nd*nmode)], shape = shape(MeanVec)), kind = RK)
1133 : complex(CK) , parameter :: InvCovMat(nd,nd,nmode) = cmplx(reshape( [ 1._RK, 0._RK, 1._RK &
1134 : , 0._RK, 2._RK, 0._RK &
1135 : , 1._RK, 0._RK, 3._RK &
1136 : , 2._RK, 0._RK, 0._RK &
1137 : , 0._RK, 2._RK, 0._RK &
1138 : , 0._RK, 0._RK, 2._RK &
1139 : ], shape = shape(InvCovMat) ), kind = RK)
1140 : complex(CK) , parameter :: LogAmplitude(nmode) = cmplx([ 3._RK, 4._RK ], kind = RK)
1141 : complex(CK) , parameter :: LogSqrtDetInvCovMat(nmode) = cmplx([-0.693147180559945_RK, -1.039720770839918_RK], kind = RK)
1142 : complex(CK) , parameter :: logProbMixMVN_ref(np) = cmplx([ -90.44996278017396_RK, -18.44996278017396_RK ], kind = RK)
1143 : complex(CK) :: logProbMixMVN(np)
1144 : real(RK) :: difference(np)
1145 :
1146 : logProbMixMVN = getlogProbMixMVN( nmode = nmode &
1147 : , nd = nd &
1148 : , np = np &
1149 : , LogAmplitude = LogAmplitude &
1150 : , MeanVec = MeanVec &
1151 : , InvCovMat = InvCovMat &
1152 : , LogSqrtDetInvCovMat = LogSqrtDetInvCovMat &
1153 : , point = point &
1154 3 : )
1155 :
1156 9 : difference = abs( (real(logProbMixMVN, kind = RK) - real(logProbMixMVN_ref, kind = RK)) / real(logProbMixMVN_ref, kind = RK) )
1157 9 : assertion = all(difference <= tolerance)
1158 :
1159 : ! LCOV_EXCL_START
1160 : if (Test%isDebugMode .and. .not. assertion) then
1161 : write(Test%outputUnit,"(*(g0,:,', '))")
1162 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN_ref ", real(logProbMixMVN_ref, kind = RK)
1163 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMixMVN ", real(logProbMixMVN, kind = RK)
1164 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
1165 : write(Test%outputUnit,"(*(g0,:,', '))")
1166 : end if
1167 : ! LCOV_EXCL_STOP
1168 :
1169 3 : end function test_getLogProbMixMVNMP_CK_1
1170 :
1171 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1172 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1173 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1174 :
1175 3 : function test_getMean_2D_1() result(assertion)
1176 :
1177 3 : use Constants_mod, only: IK, RK
1178 : implicit none
1179 :
1180 : logical :: assertion
1181 : real(RK) , parameter :: tolerance = 1.e-12_RK
1182 : integer(IK) , parameter :: nd = 3_IK
1183 : integer(IK) , parameter :: np = 5_IK
1184 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1185 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1186 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1187 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1188 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1189 : ], shape = shape(Point) )
1190 :
1191 : real(RK) , parameter :: Mean_ref(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1192 : real(RK), allocatable :: Mean(:)
1193 : real(RK), allocatable :: Difference(:)
1194 3 : Mean = getMean(nd,np,Point)
1195 :
1196 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
1197 3 : if (allocated(Difference)) deallocate(Difference); allocate(Difference, mold = Mean)
1198 :
1199 15 : Difference = abs( (Mean - Mean_ref) / Mean_ref )
1200 12 : assertion = all(Difference < tolerance)
1201 :
1202 : ! LCOV_EXCL_START
1203 : if (Test%isDebugMode .and. .not. assertion) then
1204 : write(Test%outputUnit,"(*(g0,:,', '))")
1205 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_ref ", real(Mean_ref, kind = RK)
1206 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean ", real(Mean, kind = RK)
1207 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
1208 : write(Test%outputUnit,"(*(g0,:,', '))")
1209 : end if
1210 : ! LCOV_EXCL_STOP
1211 :
1212 3 : end function test_getMean_2D_1
1213 :
1214 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1215 :
1216 3 : function test_getMean_2D_2() result(assertion)
1217 :
1218 3 : use Constants_mod, only: IK, RK
1219 : implicit none
1220 :
1221 : logical :: assertion
1222 : integer(IK) :: i
1223 : real(RK) , parameter :: tolerance = 1.e-12_RK
1224 : integer(IK) , parameter :: nd = 3_IK
1225 : integer(IK) , parameter :: np = 5_IK
1226 : integer(IK) , parameter :: Weight(np) = [(i,i=0,np-1)]
1227 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1228 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1229 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1230 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1231 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1232 : ], shape = shape(Point) )
1233 :
1234 : real(RK) , parameter :: Mean_ref(nd) = [.4601234030687523_RK, .5228363424875015_RK, .4616067715500543_RK]
1235 : real(RK), allocatable :: Mean(:)
1236 : real(RK), allocatable :: Difference(:)
1237 3 : Mean = getMean(nd, np, Point, Weight)
1238 :
1239 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
1240 3 : if (allocated(Difference)) deallocate(Difference); allocate(Difference, mold = Mean)
1241 :
1242 15 : Difference = abs( (Mean - Mean_ref) / Mean_ref )
1243 12 : assertion = all(Difference < tolerance)
1244 :
1245 : ! LCOV_EXCL_START
1246 : if (Test%isDebugMode .and. .not. assertion) then
1247 : write(Test%outputUnit,"(*(g0,:,', '))")
1248 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_ref ", real(Mean_ref, kind = RK)
1249 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean ", real(Mean, kind = RK)
1250 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
1251 : write(Test%outputUnit,"(*(g0,:,', '))")
1252 : end if
1253 : ! LCOV_EXCL_STOP
1254 :
1255 3 : end function test_getMean_2D_2
1256 :
1257 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1258 :
1259 3 : function test_getNormData_2D_1() result(assertion)
1260 :
1261 3 : use Constants_mod, only: IK, RK
1262 : implicit none
1263 :
1264 : logical :: assertion
1265 : real(RK) , parameter :: tolerance = 1.e-12_RK
1266 : integer(IK) , parameter :: nd = 3_IK
1267 : integer(IK) , parameter :: np = 5_IK
1268 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1269 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1270 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1271 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1272 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1273 : ], shape = shape(Point) )
1274 :
1275 : real(RK) , parameter :: Mean(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1276 : real(RK) , parameter :: NormData_ref(nd,np) = reshape( [ .2566442939639110_RK &
1277 : , -.3041688273160970_RK &
1278 : , -.2468837997938950_RK &
1279 : , -.4032304034245440_RK &
1280 : , -.2388698924576700_RK &
1281 : , .2996510435725079_RK &
1282 : , .2454268289201190_RK &
1283 : , -.1890219363265699E-01_RK &
1284 : , .4264152640835700_RK &
1285 : , -.4149557135527890_RK &
1286 : , .1027426859628800_RK &
1287 : , -.1422483276617770_RK &
1288 : , .3161149940933040_RK &
1289 : , .4591982274435449_RK &
1290 : , -.3369341802004060_RK &
1291 : ], shape = shape(NormData_ref) )
1292 : real(RK), allocatable :: NormData(:,:)
1293 : real(RK), allocatable :: Difference(:,:)
1294 3 : NormData = getNormData(nd,np,Mean,Point)
1295 :
1296 : ! Gfortran 7.1 fails to automatically reallocate this array. This is not implemented in Gfortran 7.0.0
1297 3 : if (allocated(Difference)) deallocate(Difference); allocate(Difference, mold = NormData_ref)
1298 :
1299 66 : Difference = abs( (transpose(NormData) - NormData_ref) / NormData_ref )
1300 63 : assertion = all(Difference <= tolerance)
1301 :
1302 : ! LCOV_EXCL_START
1303 : if (Test%isDebugMode .and. .not. assertion) then
1304 : write(Test%outputUnit,"(*(g0,:,', '))")
1305 : write(Test%outputUnit,"(*(g0,:,', '))") "NormData_ref ", real(NormData_ref, kind = RK)
1306 : write(Test%outputUnit,"(*(g0,:,', '))") "NormData ", real(transpose(NormData), kind = RK)
1307 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", real(difference, kind = RK)
1308 : write(Test%outputUnit,"(*(g0,:,', '))")
1309 : end if
1310 : ! LCOV_EXCL_STOP
1311 :
1312 3 : end function test_getNormData_2D_1
1313 :
1314 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1315 :
1316 3 : function test_getVariance_1D_1() result(assertion)
1317 :
1318 3 : use Constants_mod, only: IK, RK
1319 : implicit none
1320 :
1321 : logical :: assertion
1322 : real(RK) , parameter :: tolerance = 1.e-12_RK
1323 : real(RK) , parameter :: Point(*) = [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1324 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1325 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1326 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1327 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK ]
1328 : integer(IK) , parameter :: np = size(Point)
1329 : real(RK) , parameter :: mean = 0.436403417501334_RK
1330 : real(RK) , parameter :: variance_ref = 0.106281559054025_RK
1331 3 : real(RK) :: variance
1332 3 : real(RK) :: difference
1333 :
1334 3 : variance = getVariance(np,mean,Point)
1335 3 : difference = abs(variance - variance_ref) / variance_ref
1336 3 : assertion = difference <= tolerance
1337 :
1338 : ! LCOV_EXCL_START
1339 : if (Test%isDebugMode .and. .not. assertion) then
1340 : write(Test%outputUnit,"(*(g0,:,', '))")
1341 : write(Test%outputUnit,"(*(g0,:,', '))") "variance_ref ", variance_ref
1342 : write(Test%outputUnit,"(*(g0,:,', '))") "variance ", variance
1343 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
1344 : write(Test%outputUnit,"(*(g0,:,', '))")
1345 : end if
1346 : ! LCOV_EXCL_STOP
1347 :
1348 3 : end function test_getVariance_1D_1
1349 :
1350 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1351 :
1352 3 : function test_getVariance_1D_2() result(assertion)
1353 :
1354 3 : use Constants_mod, only: IK, RK
1355 : implicit none
1356 :
1357 : logical :: assertion
1358 : integer(IK) :: i
1359 : real(RK) , parameter :: tolerance = 1.e-12_RK
1360 : real(RK) , parameter :: Point(*) = [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1361 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1362 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1363 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1364 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK ]
1365 : integer(IK) , parameter :: np = size(Point)
1366 : integer(IK) , parameter :: Weight(np) = [(i,i=0,np-1)]
1367 : real(RK) , parameter :: mean = 0.436403417501334_RK
1368 : real(RK) , parameter :: variance_ref = .9447699069025800E-01_RK
1369 3 : real(RK) :: variance
1370 3 : real(RK) :: difference
1371 :
1372 3 : variance = getVariance(np,mean,Point,Weight,sum(Weight))
1373 3 : difference = abs(variance - variance_ref) / variance_ref
1374 3 : assertion = difference <= tolerance
1375 :
1376 : ! LCOV_EXCL_START
1377 : if (Test%isDebugMode .and. .not. assertion) then
1378 : write(Test%outputUnit,"(*(g0,:,', '))")
1379 : write(Test%outputUnit,"(*(g0,:,', '))") "variance_ref ", variance_ref
1380 : write(Test%outputUnit,"(*(g0,:,', '))") "variance ", variance
1381 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", difference
1382 : write(Test%outputUnit,"(*(g0,:,', '))")
1383 : end if
1384 : ! LCOV_EXCL_STOP
1385 :
1386 3 : end function test_getVariance_1D_2
1387 :
1388 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1389 :
1390 3 : function test_getVariance_2D_1() result(assertion)
1391 :
1392 3 : use Constants_mod, only: IK, RK
1393 : implicit none
1394 :
1395 : logical :: assertion
1396 : real(RK) , parameter :: tolerance = 1.e-12_RK
1397 : integer(IK) , parameter :: nd = 3_IK
1398 : integer(IK) , parameter :: np = 5_IK
1399 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1400 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1401 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1402 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1403 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1404 : ], shape = shape(Point) )
1405 :
1406 : real(RK) , parameter :: Mean(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1407 : real(RK), parameter :: Variance_ref(nd) = [.1402030784811636_RK, .9283846639096889E-01_RK, .1165828911170294_RK]
1408 : real(RK), allocatable :: Variance(:)
1409 : real(RK), allocatable :: Difference(:)
1410 3 : Variance = getVariance(nd,np,Mean,Point)
1411 12 : Difference = abs(Variance - Variance_ref) / Variance_ref
1412 12 : assertion = all(Difference <= tolerance)
1413 :
1414 : ! LCOV_EXCL_START
1415 : if (Test%isDebugMode .and. .not. assertion) then
1416 : write(Test%outputUnit,"(*(g0,:,', '))")
1417 : write(Test%outputUnit,"(*(g0,:,', '))") "Variance_ref ", Variance_ref
1418 : write(Test%outputUnit,"(*(g0,:,', '))") "Variance ", Variance
1419 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", Difference
1420 : write(Test%outputUnit,"(*(g0,:,', '))")
1421 : end if
1422 : ! LCOV_EXCL_STOP
1423 :
1424 3 : end function test_getVariance_2D_1
1425 :
1426 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1427 :
1428 3 : function test_getVariance_2D_2() result(assertion)
1429 :
1430 3 : use Constants_mod, only: IK, RK
1431 : implicit none
1432 :
1433 : logical :: assertion
1434 : real(RK) , parameter :: tolerance = 1.e-12_RK
1435 : integer(IK) , parameter :: nd = 3_IK
1436 : integer(IK) , parameter :: np = 5_IK
1437 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1438 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1439 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1440 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1441 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1442 : ], shape = shape(Point) )
1443 :
1444 : real(RK) , parameter :: Mean(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1445 : integer(IK) :: i
1446 : integer(IK) , parameter :: Weight(np) = [(i,i=0,np-1)]
1447 : real(RK), parameter :: Variance_ref(nd) = [.1332603228384714_RK, .1036548486974186_RK, .1075836700131122_RK]
1448 : real(RK), allocatable :: Variance(:)
1449 : real(RK), allocatable :: Difference(:)
1450 3 : Variance = getVariance(nd,np,Mean,Point,Weight)
1451 12 : Difference = abs(Variance - Variance_ref) / Variance_ref
1452 12 : assertion = all(Difference <= tolerance)
1453 :
1454 : ! LCOV_EXCL_START
1455 : if (Test%isDebugMode .and. .not. assertion) then
1456 : write(Test%outputUnit,"(*(g0,:,', '))")
1457 : write(Test%outputUnit,"(*(g0,:,', '))") "Variance_ref ", Variance_ref
1458 : write(Test%outputUnit,"(*(g0,:,', '))") "Variance ", Variance
1459 : write(Test%outputUnit,"(*(g0,:,', '))") "Difference ", Difference
1460 : write(Test%outputUnit,"(*(g0,:,', '))")
1461 : end if
1462 : ! LCOV_EXCL_STOP
1463 :
1464 3 : end function test_getVariance_2D_2
1465 :
1466 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1467 :
1468 3 : function test_getSamCholFac_1() result(assertion)
1469 :
1470 3 : use Constants_mod, only: IK, RK
1471 : implicit none
1472 :
1473 : logical :: assertion
1474 : real(RK) , parameter :: tolerance = 1.e-12_RK
1475 : integer(IK) , parameter :: nd = 3_IK
1476 : integer(IK) , parameter :: np = 5_IK
1477 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1478 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1479 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1480 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1481 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1482 : ], shape = shape(Point) )
1483 :
1484 : real(RK) , parameter :: Mean(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1485 : !integer(IK) :: i
1486 : !integer(IK) , parameter :: Weight(np) = [(i,i=0,np-1)]
1487 : real(RK), parameter :: CholeskyLower_ref(nd,nd) = transpose(reshape([ +0.140203078481164_RK, +0.029035771029083_RK, -0.031754793421423_RK &
1488 : , +0.077545140669886_RK, +0.092838466390969_RK, -0.043469498536710_RK &
1489 : , -0.084806768876261_RK, -0.125205309782421_RK, +0.116582891117029_RK &
1490 : ] , shape = shape(CholeskyLower_ref) ))
1491 : real(RK), parameter :: Diagonal_ref(nd) = [ +0.374437015372631_RK, +0.294661191115248_RK, +0.306127969111099_RK ]
1492 : real(RK) :: CholeskyLower_diff(nd,nd)
1493 : real(RK) :: CholeskyLower(nd,nd)
1494 : real(RK) :: Diagonal_diff(nd)
1495 : real(RK) :: Diagonal(nd)
1496 :
1497 3 : call getSamCholFac(nd,np,Mean,Point,CholeskyLower,Diagonal)
1498 39 : CholeskyLower_diff = abs(CholeskyLower - CholeskyLower_ref) / abs(CholeskyLower_ref)
1499 12 : Diagonal_diff = abs(Diagonal - Diagonal_ref) / abs(Diagonal_ref)
1500 48 : assertion = all(CholeskyLower_diff <= tolerance) .and. all(Diagonal_diff <= tolerance)
1501 :
1502 : ! LCOV_EXCL_START
1503 : if (Test%isDebugMode .and. .not. assertion) then
1504 : write(Test%outputUnit,"(*(g0,:,', '))")
1505 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyLower_ref ", CholeskyLower_ref
1506 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyLower ", CholeskyLower
1507 : write(Test%outputUnit,"(*(g0,:,', '))") "CholeskyLower_diff ", CholeskyLower_diff
1508 : write(Test%outputUnit,"(*(g0,:,', '))")
1509 : write(Test%outputUnit,"(*(g0,:,', '))") "Diagonal_ref ", Diagonal_ref
1510 : write(Test%outputUnit,"(*(g0,:,', '))") "Diagonal ", Diagonal
1511 : write(Test%outputUnit,"(*(g0,:,', '))") "Diagonal_diff ", Diagonal_diff
1512 : write(Test%outputUnit,"(*(g0,:,', '))")
1513 : end if
1514 : ! LCOV_EXCL_STOP
1515 :
1516 3 : end function test_getSamCholFac_1
1517 :
1518 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1519 :
1520 3 : function test_getSamCovMean_1() result(assertion)
1521 :
1522 3 : use Constants_mod, only: IK, RK
1523 : implicit none
1524 :
1525 : logical :: assertion, assertionCurrent
1526 : real(RK) , parameter :: tolerance = 1.e-11_RK
1527 : integer(IK) , parameter :: nd = 3_IK
1528 : integer(IK) , parameter :: np = 5_IK
1529 : real(RK) , parameter :: Point(np,nd) = transpose(reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1530 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1531 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1532 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1533 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1534 : ], shape = [nd,np] ) )
1535 :
1536 : real(RK), parameter :: Mean_ref(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1537 : real(RK), parameter :: CovMat_ref(nd,nd) = reshape( [ +0.140203078481164_RK, +0.029035771029083_RK, -0.031754793421423_RK &
1538 : , +0.029035771029083_RK, +0.092838466390969_RK, -0.043469498536710_RK &
1539 : , -0.031754793421423_RK, -0.043469498536710_RK, +0.116582891117029_RK &
1540 : ] , shape = shape(CovMat_ref) )
1541 : real(RK), parameter :: InvCovMat_ref(nd,nd) = reshape( [ +7.831154269923546_RK, -1.757283899385404_RK, +1.477819211284223_RK &
1542 : , -1.757283899385404_RK, 13.444000278032743_RK, +4.534128105255763_RK &
1543 : , +1.477819211284223_RK, +4.534128105255763_RK, 10.670726269401085_RK &
1544 : ] , shape = shape(CovMat_ref) )
1545 : real(RK), parameter :: MahalSq_ref(np) = [3.178088257444105_RK, 1.653804994353691_RK, 2.669296951657121_RK, 1.8980344204538842_RK, 2.6007753760912014_RK]
1546 : real(RK), parameter :: sqrtDetInvCovMat_ref = 29.607059382128476_RK
1547 : real(RK) :: CovMat(nd,nd)
1548 : real(RK) :: CovMat_diff(nd,nd)
1549 3 : real(RK) :: Mean_diff(nd), MahalSq_diff(np), InvCovMat_diff(nd,nd), sqrtDetInvCovMat_diff
1550 3 : real(RK) :: Mean(nd), MahalSq(np), InvCovMat(nd,nd), sqrtDetInvCovMat
1551 :
1552 3 : assertion = .true.
1553 :
1554 3 : call getSamCovMean(np,nd,Point,CovMat,Mean,MahalSq,InvCovMat,sqrtDetInvCovMat)
1555 39 : CovMat_diff = abs(CovMat - CovMat_ref) / abs(CovMat_ref)
1556 39 : assertionCurrent = all(CovMat_diff <= tolerance)
1557 3 : assertion = assertion .and. assertionCurrent
1558 :
1559 : ! LCOV_EXCL_START
1560 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1561 : write(Test%outputUnit,"(*(g0,:,', '))")
1562 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMat_ref ", CovMat_ref
1563 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMat ", CovMat
1564 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMat_diff ", CovMat_diff
1565 : write(Test%outputUnit,"(*(g0,:,', '))")
1566 : end if
1567 : ! LCOV_EXCL_STOP
1568 :
1569 12 : Mean_diff = abs(Mean - Mean_ref) / abs(Mean_ref)
1570 12 : assertionCurrent = all(Mean_diff <= tolerance)
1571 3 : assertion = assertion .and. assertionCurrent
1572 :
1573 : ! LCOV_EXCL_START
1574 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1575 : write(Test%outputUnit,"(*(g0,:,', '))")
1576 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_ref ", Mean_ref
1577 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean ", Mean
1578 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_diff ", Mean_diff
1579 : write(Test%outputUnit,"(*(g0,:,', '))")
1580 : end if
1581 : ! LCOV_EXCL_STOP
1582 :
1583 18 : MahalSq_diff = abs(MahalSq - MahalSq_ref) / abs(MahalSq_ref)
1584 18 : assertionCurrent = all(MahalSq_diff <= tolerance)
1585 3 : assertion = assertion .and. assertionCurrent
1586 :
1587 : ! LCOV_EXCL_START
1588 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1589 : write(Test%outputUnit,"(*(g0,:,', '))")
1590 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq_ref ", MahalSq_ref
1591 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq ", MahalSq
1592 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq_diff ", MahalSq_diff
1593 : write(Test%outputUnit,"(*(g0,:,', '))")
1594 : end if
1595 : ! LCOV_EXCL_STOP
1596 :
1597 39 : InvCovMat_diff = abs(InvCovMat - InvCovMat_ref) / abs(InvCovMat_ref)
1598 39 : assertionCurrent = all(InvCovMat_diff <= tolerance)
1599 3 : assertion = assertion .and. assertionCurrent
1600 :
1601 : ! LCOV_EXCL_START
1602 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1603 : write(Test%outputUnit,"(*(g0,:,', '))")
1604 : write(Test%outputUnit,"(*(g0,:,', '))") "InvCovMat_ref ", InvCovMat_ref
1605 : write(Test%outputUnit,"(*(g0,:,', '))") "InvCovMat ", InvCovMat
1606 : write(Test%outputUnit,"(*(g0,:,', '))") "InvCovMat_diff ", InvCovMat_diff
1607 : write(Test%outputUnit,"(*(g0,:,', '))")
1608 : end if
1609 : ! LCOV_EXCL_STOP
1610 :
1611 3 : sqrtDetInvCovMat_diff = abs(sqrtDetInvCovMat - sqrtDetInvCovMat_ref) / abs(sqrtDetInvCovMat_ref)
1612 3 : assertionCurrent = sqrtDetInvCovMat_diff <= tolerance
1613 3 : assertion = assertion .and. assertionCurrent
1614 :
1615 : ! LCOV_EXCL_START
1616 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1617 : write(Test%outputUnit,"(*(g0,:,', '))")
1618 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvCovMat_ref ", sqrtDetInvCovMat_ref
1619 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvCovMat ", sqrtDetInvCovMat
1620 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvCovMat_diff ", sqrtDetInvCovMat_diff
1621 : write(Test%outputUnit,"(*(g0,:,', '))")
1622 : end if
1623 : ! LCOV_EXCL_STOP
1624 :
1625 3 : end function test_getSamCovMean_1
1626 :
1627 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1628 :
1629 3 : function test_getSamCovMeanTrans_1() result(assertion)
1630 :
1631 3 : use Constants_mod, only: IK, RK
1632 : implicit none
1633 :
1634 : logical :: assertion, assertionCurrent
1635 : real(RK) , parameter :: tolerance = 1.e-11_RK
1636 : integer(IK) , parameter :: nd = 3_IK
1637 : integer(IK) , parameter :: np = 5_IK
1638 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1639 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1640 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1641 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1642 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1643 : ], shape = shape(Point) )
1644 : real(RK), parameter :: Mean_ref(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1645 : real(RK), parameter :: CovMat_ref(nd,nd) = reshape( [ +0.140203078481164_RK, +0.029035771029083_RK, -0.031754793421423_RK &
1646 : , +0.029035771029083_RK, +0.092838466390969_RK, -0.043469498536710_RK &
1647 : , -0.031754793421423_RK, -0.043469498536710_RK, +0.116582891117029_RK &
1648 : ] , shape = shape(CovMat_ref) )
1649 : real(RK), parameter :: InvCovMat_ref(nd,nd) = reshape( [ +7.831154269923546_RK, -1.757283899385404_RK, +1.477819211284223_RK &
1650 : , -1.757283899385404_RK, 13.444000278032743_RK, +4.534128105255763_RK &
1651 : , +1.477819211284223_RK, +4.534128105255763_RK, 10.670726269401085_RK &
1652 : ] , shape = shape(CovMat_ref) )
1653 : real(RK), parameter :: MahalSq_ref(np) = [3.178088257444105_RK, 1.653804994353691_RK, 2.669296951657121_RK, 1.8980344204538842_RK, 2.6007753760912014_RK]
1654 : real(RK), parameter :: sqrtDetInvCovMat_ref = 29.607059382128476_RK
1655 : real(RK) :: CovMat(nd,nd)
1656 : real(RK) :: CovMat_diff(nd,nd)
1657 3 : real(RK) :: Mean_diff(nd), MahalSq_diff(np), InvCovMat_diff(nd,nd), sqrtDetInvCovMat_diff
1658 3 : real(RK) :: Mean(nd), MahalSq(np), InvCovMat(nd,nd), sqrtDetInvCovMat
1659 :
1660 3 : assertion = .true.
1661 :
1662 3 : call getSamCovMeanTrans(np,nd,Point,CovMat,Mean,MahalSq,InvCovMat,sqrtDetInvCovMat)
1663 39 : CovMat_diff = abs(CovMat - CovMat_ref) / abs(CovMat_ref)
1664 39 : assertionCurrent = all(CovMat_diff <= tolerance)
1665 3 : assertion = assertion .and. assertionCurrent
1666 :
1667 : ! LCOV_EXCL_START
1668 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1669 : write(Test%outputUnit,"(*(g0,:,', '))")
1670 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMat_ref ", CovMat_ref
1671 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMat ", CovMat
1672 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMat_diff ", CovMat_diff
1673 : write(Test%outputUnit,"(*(g0,:,', '))")
1674 : end if
1675 : ! LCOV_EXCL_STOP
1676 :
1677 12 : Mean_diff = abs(Mean - Mean_ref) / abs(Mean_ref)
1678 12 : assertionCurrent = all(Mean_diff <= tolerance)
1679 3 : assertion = assertion .and. assertionCurrent
1680 :
1681 : ! LCOV_EXCL_START
1682 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1683 : write(Test%outputUnit,"(*(g0,:,', '))")
1684 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_ref ", Mean_ref
1685 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean ", Mean
1686 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_diff ", Mean_diff
1687 : write(Test%outputUnit,"(*(g0,:,', '))")
1688 : end if
1689 : ! LCOV_EXCL_STOP
1690 :
1691 18 : MahalSq_diff = abs(MahalSq - MahalSq_ref) / abs(MahalSq_ref)
1692 18 : assertionCurrent = all(MahalSq_diff <= tolerance)
1693 3 : assertion = assertion .and. assertionCurrent
1694 :
1695 : ! LCOV_EXCL_START
1696 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1697 : write(Test%outputUnit,"(*(g0,:,', '))")
1698 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq_ref ", MahalSq_ref
1699 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq ", MahalSq
1700 : write(Test%outputUnit,"(*(g0,:,', '))") "MahalSq_diff ", MahalSq_diff
1701 : write(Test%outputUnit,"(*(g0,:,', '))")
1702 : end if
1703 : ! LCOV_EXCL_STOP
1704 :
1705 39 : InvCovMat_diff = abs(InvCovMat - InvCovMat_ref) / abs(InvCovMat_ref)
1706 39 : assertionCurrent = all(InvCovMat_diff <= tolerance)
1707 3 : assertion = assertion .and. assertionCurrent
1708 :
1709 : ! LCOV_EXCL_START
1710 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1711 : write(Test%outputUnit,"(*(g0,:,', '))")
1712 : write(Test%outputUnit,"(*(g0,:,', '))") "InvCovMat_ref ", InvCovMat_ref
1713 : write(Test%outputUnit,"(*(g0,:,', '))") "InvCovMat ", InvCovMat
1714 : write(Test%outputUnit,"(*(g0,:,', '))") "InvCovMat_diff ", InvCovMat_diff
1715 : write(Test%outputUnit,"(*(g0,:,', '))")
1716 : end if
1717 : ! LCOV_EXCL_STOP
1718 :
1719 3 : sqrtDetInvCovMat_diff = abs(sqrtDetInvCovMat - sqrtDetInvCovMat_ref) / abs(sqrtDetInvCovMat_ref)
1720 3 : assertionCurrent = sqrtDetInvCovMat_diff <= tolerance
1721 3 : assertion = assertion .and. assertionCurrent
1722 :
1723 : ! LCOV_EXCL_START
1724 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1725 : write(Test%outputUnit,"(*(g0,:,', '))")
1726 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvCovMat_ref ", sqrtDetInvCovMat_ref
1727 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvCovMat ", sqrtDetInvCovMat
1728 : write(Test%outputUnit,"(*(g0,:,', '))") "sqrtDetInvCovMat_diff ", sqrtDetInvCovMat_diff
1729 : write(Test%outputUnit,"(*(g0,:,', '))")
1730 : end if
1731 : ! LCOV_EXCL_STOP
1732 :
1733 3 : end function test_getSamCovMeanTrans_1
1734 :
1735 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1736 :
1737 3 : function test_getSamCovUpperMeanTrans_1() result(assertion)
1738 :
1739 3 : use Constants_mod, only: IK, RK
1740 : implicit none
1741 :
1742 : integer(IK) :: i, j
1743 : logical :: assertion, assertionCurrent
1744 : real(RK) , parameter :: tolerance = 1.e-12_RK
1745 : integer(IK) , parameter :: nd = 3_IK
1746 : integer(IK) , parameter :: np = 5_IK
1747 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1748 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1749 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1750 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1751 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1752 : ], shape = shape(Point) )
1753 : real(RK), parameter :: Mean_ref(nd) = [0.449401794055698_RK, 0.336001673693518_RK, 0.523806784754785_RK]
1754 : real(RK), parameter :: CovMatUpper_ref(nd,nd) = reshape([ +0.140203078481164_RK, +0.029035771029083_RK, -0.031754793421423_RK &
1755 : , +0.029035771029083_RK, +0.092838466390969_RK, -0.043469498536710_RK &
1756 : , -0.031754793421423_RK, -0.043469498536710_RK, +0.116582891117029_RK &
1757 : ] , shape = shape(CovMatUpper_ref) )
1758 : real(RK) :: CovMatUpper(nd,nd)
1759 : real(RK) :: CovMatUpper_diff(nd,nd)
1760 : real(RK) :: Mean_diff(nd)
1761 : real(RK) :: Mean(nd)
1762 :
1763 3 : assertion = .true.
1764 :
1765 3 : CovMatUpper = CovMatUpper_ref ! to avoid un-initialization error in the computation of `CovMatUpper_diff`.
1766 3 : call getSamCovUpperMeanTrans(np,nd,Point,CovMatUpper,Mean)
1767 39 : CovMatUpper_diff = abs(CovMatUpper - CovMatUpper_ref) / abs(CovMatUpper_ref)
1768 48 : assertionCurrent = all([ ((CovMatUpper_diff(i,j) <= tolerance, i=1,j), j=1,nd) ])
1769 3 : assertion = assertion .and. assertionCurrent
1770 :
1771 : ! LCOV_EXCL_START
1772 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1773 : write(Test%outputUnit,"(*(g0,:,', '))")
1774 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpper_ref ", ((CovMatUpper_ref(i,j), i=1,j), j=1,nd)
1775 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpper ", ((CovMatUpper(i,j), i=1,j), j=1,nd)
1776 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpper_diff ", ((CovMatUpper_diff(i,j), i=1,j), j=1,nd)
1777 : write(Test%outputUnit,"(*(g0,:,', '))")
1778 : end if
1779 : ! LCOV_EXCL_STOP
1780 :
1781 12 : Mean_diff = abs(Mean - Mean_ref) / abs(Mean_ref)
1782 12 : assertionCurrent = all(Mean_diff <= tolerance)
1783 3 : assertion = assertion .and. assertionCurrent
1784 :
1785 : ! LCOV_EXCL_START
1786 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1787 : write(Test%outputUnit,"(*(g0,:,', '))")
1788 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_ref ", Mean_ref
1789 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean ", Mean
1790 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_diff ", Mean_diff
1791 : write(Test%outputUnit,"(*(g0,:,', '))")
1792 : end if
1793 : ! LCOV_EXCL_STOP
1794 :
1795 3 : end function test_getSamCovUpperMeanTrans_1
1796 :
1797 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1798 :
1799 3 : function test_getWeiSamCovUppMeanTrans_1() result(assertion)
1800 :
1801 3 : use Constants_mod, only: IK, RK
1802 : implicit none
1803 :
1804 : integer(IK) :: i, j
1805 : logical :: assertion, assertionCurrent
1806 : real(RK) , parameter :: tolerance = 1.e-12_RK
1807 : integer(IK) , parameter :: nd = 3_IK
1808 : integer(IK) , parameter :: np = 5_IK
1809 : integer(IK) , parameter :: Weight(np) = [(i,i=1,np)]
1810 : real(RK) , parameter :: Point(nd,np) = reshape( [ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1811 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1812 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1813 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1814 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1815 : ], shape = shape(Point) )
1816 : real(RK), parameter :: Mean_ref(nd) = [.4565495333977344_RK, .4605581195561738_RK, .4823401092849646_RK]
1817 : real(RK), parameter :: CovMatUpper_ref(nd,nd) = reshape([ +.1256706333432408E+00_RK, +.4589762909514189E-01_RK, -.2021812246706341E-01_RK &
1818 : , +.4589762909514189E-01_RK, +.7653806291082801E-01_RK, -.6048750626649663E-01_RK &
1819 : , -.2021812246706341E-01_RK, -.6048750626649663E-01_RK, +.1006280226405955E+00_RK &
1820 : ] , shape = shape(CovMatUpper_ref) )
1821 : real(RK) :: CovMatUpper(nd,nd)
1822 : real(RK) :: CovMatUpper_diff(nd,nd)
1823 : real(RK) :: Mean_diff(nd)
1824 : real(RK) :: Mean(nd)
1825 :
1826 3 : assertion = .true.
1827 :
1828 3 : call getWeiSamCovUppMeanTrans(np,sum(Weight),nd,Point,Weight,CovMatUpper,Mean)
1829 3 : assertionCurrent = .true.
1830 12 : do j = 1, nd
1831 30 : do i = 1, j
1832 18 : CovMatUpper_diff(i,j) = abs(CovMatUpper(i,j) - CovMatUpper_ref(i,j)) / abs(CovMatUpper_ref(i,j))
1833 27 : assertionCurrent = assertionCurrent .and. CovMatUpper_diff(i,j) <= tolerance
1834 : end do
1835 : end do
1836 3 : assertion = assertion .and. assertionCurrent
1837 :
1838 : ! LCOV_EXCL_START
1839 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1840 : write(Test%outputUnit,"(*(g0,:,', '))")
1841 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpper_ref ", ((CovMatUpper_ref(i,j), i=1,j), j=1,nd)
1842 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpper ", ((CovMatUpper(i,j), i=1,j), j=1,nd)
1843 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpper_diff ", ((CovMatUpper_diff(i,j), i=1,j), j=1,nd)
1844 : write(Test%outputUnit,"(*(g0,:,', '))")
1845 : end if
1846 : ! LCOV_EXCL_STOP
1847 :
1848 12 : Mean_diff = abs(Mean - Mean_ref) / abs(Mean_ref)
1849 12 : assertionCurrent = all(Mean_diff <= tolerance)
1850 3 : assertion = assertion .and. assertionCurrent
1851 :
1852 : ! LCOV_EXCL_START
1853 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1854 : write(Test%outputUnit,"(*(g0,:,', '))")
1855 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_ref ", Mean_ref
1856 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean ", Mean
1857 : write(Test%outputUnit,"(*(g0,:,', '))") "Mean_diff ", Mean_diff
1858 : write(Test%outputUnit,"(*(g0,:,', '))")
1859 : end if
1860 : ! LCOV_EXCL_STOP
1861 :
1862 3 : end function test_getWeiSamCovUppMeanTrans_1
1863 :
1864 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1865 :
1866 3 : function test_mergeMeanCov_1() result(assertion)
1867 :
1868 3 : use Constants_mod, only: IK, RK
1869 : implicit none
1870 :
1871 : integer(IK) :: i, j
1872 : logical :: assertion, assertionCurrent
1873 : real(RK) , parameter :: tolerance = 1.e-12_RK
1874 : integer(IK) , parameter :: nd = 3_IK
1875 :
1876 : integer(IK) , parameter :: npA = 5_IK
1877 : real(RK) :: MeanA(nd), CovMatA(nd,nd)
1878 : real(RK) , parameter :: PointA(nd,npA) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1879 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1880 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1881 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1882 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1883 : ], shape = shape(PointA) )
1884 :
1885 : integer(IK) , parameter :: npB = 10_IK
1886 : real(RK) :: MeanB(nd), CovMatB(nd,nd)
1887 : real(RK) , parameter :: PointB(nd,npB) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1888 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1889 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1890 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1891 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1892 : , 1.706046088019609_RK, 2.031832846377421_RK, 3.276922984960890_RK &
1893 : , 1.046171390631154_RK, 2.097131781235848_RK, 3.823457828327293_RK &
1894 : , 1.694828622975817_RK, 2.317099480060861_RK, 3.950222048838355_RK &
1895 : , 1.034446080502909_RK, 2.438744359656398_RK, 3.381558457093008_RK &
1896 : , 1.765516788149002_RK, 2.795199901137063_RK, 3.186872604554379_RK &
1897 : ], shape = shape(PointB) )
1898 :
1899 :
1900 : integer(IK) , parameter :: npAB = npA + npB
1901 : real(RK) :: PointAB(nd,npAB)
1902 : real(RK) :: CovMatAB_diff(nd,nd)
1903 : real(RK) , parameter :: CovMatAB_refMATLAB(nd,nd) = reshape( [ 0.358269305364807_RK, 0.501078279929690_RK, 0.687067319924494_RK &
1904 : , 0.501078279929690_RK, 1.031956780716069_RK, 1.391311858397106_RK &
1905 : , 0.687067319924494_RK, 1.391311858397106_RK, 2.242785335243168_RK &
1906 : ], shape = shape(CovMatAB_refMATLAB) )
1907 : real(RK) :: CovMatAB_ref(nd,nd)
1908 : real(RK) :: CovMatAB(nd,nd)
1909 : real(RK) :: MeanAB_diff(nd)
1910 : real(RK) :: MeanAB_ref(nd)
1911 : real(RK) :: MeanAB(nd)
1912 :
1913 3 : assertion = .true.
1914 :
1915 3 : call getSamCovMeanTrans(npA,nd,PointA,CovMatA,MeanA)
1916 3 : call getSamCovMeanTrans(npB,nd,PointB,CovMatB,MeanB)
1917 :
1918 63 : PointAB(:,1:npA) = PointA
1919 123 : PointAB(:,npA+1:npAB) = PointB
1920 :
1921 3 : call getSamCovMeanTrans(npAB,nd,PointAB,CovMatAB_ref,MeanAB_ref)
1922 :
1923 : !do i = 1, 10000000
1924 3 : call mergeMeanCov(nd,npA,MeanA,CovMatA,npB,MeanB,CovMatB,MeanAB,CovMatAB)
1925 : !end do
1926 :
1927 39 : CovMatAB_diff = abs(CovMatAB - CovMatAB_ref) / abs(CovMatAB_ref)
1928 48 : assertionCurrent = all([ ((CovMatAB_diff(i,j) <= tolerance, i=1,j), j=1,nd) ])
1929 3 : assertion = assertion .and. assertionCurrent
1930 :
1931 : ! LCOV_EXCL_START
1932 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1933 : write(Test%outputUnit,"(*(g0,:,', '))")
1934 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatAB_refMATLAB ", ((CovMatAB_refMATLAB(i,j), i=1,nd), j=1,nd)
1935 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatAB_ref ", ((CovMatAB_ref(i,j), i=1,nd), j=1,nd)
1936 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatAB_diff ", ((CovMatAB_refMATLAB(i,j)-CovMatAB_ref(i,j), i=1,nd), j=1,nd)
1937 : write(Test%outputUnit,"(*(g0,:,', '))")
1938 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatAB_ref ", ((CovMatAB_ref(i,j), i=1,nd), j=1,nd)
1939 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatAB ", ((CovMatAB(i,j), i=1,nd), j=1,nd)
1940 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatAB_diff ", ((CovMatAB_diff(i,j), i=1,nd), j=1,nd)
1941 : write(Test%outputUnit,"(*(g0,:,', '))")
1942 : end if
1943 : ! LCOV_EXCL_STOP
1944 :
1945 12 : MeanAB_diff = abs(MeanAB - MeanAB_ref) / abs(MeanAB_ref)
1946 12 : assertionCurrent = all(MeanAB_diff <= tolerance)
1947 3 : assertion = assertion .and. assertionCurrent
1948 :
1949 : ! LCOV_EXCL_START
1950 : if (Test%isDebugMode .and. .not. assertionCurrent) then
1951 : write(Test%outputUnit,"(*(g0,:,', '))")
1952 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_ref ", MeanAB_ref
1953 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB ", MeanAB
1954 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_diff ", MeanAB_diff
1955 : write(Test%outputUnit,"(*(g0,:,', '))")
1956 : end if
1957 : ! LCOV_EXCL_STOP
1958 :
1959 3 : end function test_mergeMeanCov_1
1960 :
1961 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1962 :
1963 : ! function test_mergeMeanCovUpperSlow_1() result(assertion)
1964 : !
1965 : ! use Constants_mod, only: IK, RK
1966 : ! implicit none
1967 : !
1968 : ! integer(IK) :: i, j
1969 : ! logical :: assertion, assertionCurrent
1970 : ! real(RK) , parameter :: tolerance = 1.e-12_RK
1971 : ! integer(IK) , parameter :: nd = 3_IK
1972 : !
1973 : ! integer(IK) , parameter :: npA = 5_IK
1974 : ! real(RK) :: MeanA(nd), CovMatUpperA(nd,nd)
1975 : ! real(RK) , parameter :: PointA(nd,npA) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1976 : ! , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1977 : ! , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1978 : ! , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1979 : ! , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1980 : ! ], shape = shape(PointA) )
1981 : !
1982 : ! integer(IK) , parameter :: npB = 10_IK
1983 : ! real(RK) :: MeanB(nd), CovMatUpperB(nd,nd)
1984 : ! real(RK) , parameter :: PointB(nd,npB) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
1985 : ! , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
1986 : ! , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
1987 : ! , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
1988 : ! , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
1989 : ! , 1.706046088019609_RK, 2.031832846377421_RK, 3.276922984960890_RK &
1990 : ! , 1.046171390631154_RK, 2.097131781235848_RK, 3.823457828327293_RK &
1991 : ! , 1.694828622975817_RK, 2.317099480060861_RK, 3.950222048838355_RK &
1992 : ! , 1.034446080502909_RK, 2.438744359656398_RK, 3.381558457093008_RK &
1993 : ! , 1.765516788149002_RK, 2.795199901137063_RK, 3.186872604554379_RK &
1994 : ! ], shape = shape(PointB) )
1995 : !
1996 : !
1997 : ! integer(IK) , parameter :: npAB = npA + npB
1998 : ! real(RK) :: PointAB(nd,npAB)
1999 : ! real(RK) :: CovMatUpperAB_diff(nd,nd)
2000 : ! real(RK) , parameter :: CovMatUpperAB_refMATLAB(nd,nd) = reshape([ 0.358269305364807_RK, 0.501078279929690_RK, 0.687067319924494_RK &
2001 : ! , 0.501078279929690_RK, 1.031956780716069_RK, 1.391311858397106_RK &
2002 : ! , 0.687067319924494_RK, 1.391311858397106_RK, 2.242785335243168_RK &
2003 : ! ], shape = shape(CovMatUpperAB_refMATLAB) )
2004 : ! real(RK) :: CovMatUpperAB_ref(nd,nd)
2005 : ! real(RK) :: CovMatUpperAB(nd,nd)
2006 : ! real(RK) :: MeanAB_diff(nd)
2007 : ! real(RK) :: MeanAB_ref(nd)
2008 : ! real(RK) :: MeanAB(nd)
2009 : !
2010 : ! assertion = .true.
2011 : !
2012 : ! call getSamCovUpperMeanTrans(npA,nd,PointA,CovMatUpperA,MeanA)
2013 : ! call getSamCovUpperMeanTrans(npB,nd,PointB,CovMatUpperB,MeanB)
2014 : !
2015 : ! PointAB(:,1:npA) = PointA
2016 : ! PointAB(:,npA+1:npAB) = PointB
2017 : !
2018 : ! call getSamCovUpperMeanTrans(npAB,nd,PointAB,CovMatUpperAB_ref,MeanAB_ref)
2019 : !
2020 : ! !do i = 1, 10000000
2021 : ! call mergeMeanCovUpperSlow(nd,npA,MeanA,CovMatUpperA,npB,MeanB,CovMatUpperB,MeanAB,CovMatUpperAB)
2022 : ! !end do
2023 : ! assertionCurrent = .true.
2024 : ! do j = 1, nd
2025 : ! do i = 1, j
2026 : ! CovMatUpperAB_diff(i,j) = abs(CovMatUpperAB(i,j) - CovMatUpperAB_ref(i,j)) / abs(CovMatUpperAB_ref(i,j))
2027 : ! assertionCurrent = assertionCurrent .and. CovMatUpperAB_diff(i,j) <= tolerance
2028 : ! end do
2029 : ! end do
2030 : ! assertion = assertion .and. assertionCurrent
2031 : !
2032 : ! LCOV_EXCL_START
2033 : ! if (Test%isDebugMode .and. .not. assertionCurrent) then
2034 : ! write(Test%outputUnit,"(*(g0,:,', '))")
2035 : ! write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_refMATLAB ", ((CovMatUpperAB_refMATLAB(i,j), i=1,j), j=1,nd)
2036 : ! write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_ref ", ((CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2037 : ! write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_diff ", ((CovMatUpperAB_refMATLAB(i,j)-CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2038 : ! write(Test%outputUnit,"(*(g0,:,', '))")
2039 : ! write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_ref ", ((CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2040 : ! write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB ", ((CovMatUpperAB(i,j), i=1,j), j=1,nd)
2041 : ! write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_diff ", ((CovMatUpperAB_diff(i,j), i=1,j), j=1,nd)
2042 : ! write(Test%outputUnit,"(*(g0,:,', '))")
2043 : ! end if
2044 : ! LCOV_EXCL_STOP
2045 : !
2046 : ! MeanAB_diff = abs(MeanAB - MeanAB_ref) / abs(MeanAB_ref)
2047 : ! assertionCurrent = all(MeanAB_diff <= tolerance)
2048 : ! assertion = assertion .and. assertionCurrent
2049 : !
2050 : ! LCOV_EXCL_START
2051 : ! if (Test%isDebugMode .and. .not. assertionCurrent) then
2052 : ! write(Test%outputUnit,"(*(g0,:,', '))")
2053 : ! write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_ref ", MeanAB_ref
2054 : ! write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB ", MeanAB
2055 : ! write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_diff ", MeanAB_diff
2056 : ! write(Test%outputUnit,"(*(g0,:,', '))")
2057 : ! end if
2058 : ! LCOV_EXCL_STOP
2059 : !
2060 : ! end function test_mergeMeanCovUpperSlow_1
2061 :
2062 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2063 :
2064 3 : function test_mergeMeanCovUpper_1() result(assertion)
2065 :
2066 3 : use Constants_mod, only: IK, RK
2067 : implicit none
2068 :
2069 : integer(IK) :: i, j
2070 : logical :: assertion, assertionCurrent
2071 : real(RK) , parameter :: tolerance = 1.e-12_RK
2072 : integer(IK) , parameter :: nd = 3_IK
2073 :
2074 : integer(IK) , parameter :: npA = 5_IK
2075 : real(RK) :: MeanA(nd), CovMatUpperA(nd,nd)
2076 : real(RK) , parameter :: PointA(nd,npA) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
2077 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
2078 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
2079 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
2080 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
2081 : ], shape = shape(PointA) )
2082 :
2083 : integer(IK) , parameter :: npB = 10_IK
2084 : real(RK) :: MeanB(nd), CovMatUpperB(nd,nd)
2085 : real(RK) , parameter :: PointB(nd,npB) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
2086 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
2087 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
2088 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
2089 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
2090 : , 1.706046088019609_RK, 2.031832846377421_RK, 3.276922984960890_RK &
2091 : , 1.046171390631154_RK, 2.097131781235848_RK, 3.823457828327293_RK &
2092 : , 1.694828622975817_RK, 2.317099480060861_RK, 3.950222048838355_RK &
2093 : , 1.034446080502909_RK, 2.438744359656398_RK, 3.381558457093008_RK &
2094 : , 1.765516788149002_RK, 2.795199901137063_RK, 3.186872604554379_RK &
2095 : ], shape = shape(PointB) )
2096 :
2097 :
2098 : integer(IK) , parameter :: npAB = npA + npB
2099 : real(RK) :: PointAB(nd,npAB)
2100 : real(RK) :: CovMatUpperAB_diff(nd,nd)
2101 : real(RK) , parameter :: CovMatUpperAB_refMATLAB(nd,nd) = reshape([ 0.358269305364807_RK, 0.501078279929690_RK, 0.687067319924494_RK &
2102 : , 0.501078279929690_RK, 1.031956780716069_RK, 1.391311858397106_RK &
2103 : , 0.687067319924494_RK, 1.391311858397106_RK, 2.242785335243168_RK &
2104 : ], shape = shape(CovMatUpperAB_refMATLAB) )
2105 : real(RK) :: CovMatUpperAB_ref(nd,nd)
2106 : real(RK) :: CovMatUpperAB(nd,nd)
2107 : real(RK) :: MeanAB_diff(nd)
2108 : real(RK) :: MeanAB_ref(nd)
2109 : real(RK) :: MeanAB(nd)
2110 :
2111 3 : assertion = .true.
2112 :
2113 3 : call getSamCovUpperMeanTrans(npA,nd,PointA,CovMatUpperA,MeanA)
2114 3 : call getSamCovUpperMeanTrans(npB,nd,PointB,CovMatUpperB,MeanB)
2115 :
2116 63 : PointAB(:,1:npA) = PointA
2117 123 : PointAB(:,npA+1:npAB) = PointB
2118 :
2119 3 : call getSamCovUpperMeanTrans(npAB,nd,PointAB,CovMatUpperAB_ref,MeanAB_ref)
2120 :
2121 : !do i = 1, 10000000
2122 3 : call mergeMeanCovUpper(nd,npA,MeanA,CovMatUpperA,npB,MeanB,CovMatUpperB,MeanAB,CovMatUpperAB)
2123 : !end do
2124 :
2125 3 : assertionCurrent = .true.
2126 12 : do j = 1, nd
2127 30 : do i = 1, j
2128 18 : CovMatUpperAB_diff(i,j) = abs(CovMatUpperAB(i,j) - CovMatUpperAB_ref(i,j)) / abs(CovMatUpperAB_ref(i,j))
2129 27 : assertionCurrent = assertionCurrent .and. CovMatUpperAB_diff(i,j) <= tolerance
2130 : end do
2131 : end do
2132 3 : assertion = assertion .and. assertionCurrent
2133 :
2134 : ! LCOV_EXCL_START
2135 : if (Test%isDebugMode .and. .not. assertionCurrent) then
2136 : write(Test%outputUnit,"(*(g0,:,', '))")
2137 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_refMATLAB ", ((CovMatUpperAB_refMATLAB(i,j), i=1,nd), j=1,nd)
2138 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_ref ", ((CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2139 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_diff ", ((CovMatUpperAB_refMATLAB(i,j)-CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2140 : write(Test%outputUnit,"(*(g0,:,', '))")
2141 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_ref ", ((CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2142 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB ", ((CovMatUpperAB(i,j), i=1,j), j=1,nd)
2143 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_diff ", ((CovMatUpperAB_diff(i,j), i=1,j), j=1,nd)
2144 : write(Test%outputUnit,"(*(g0,:,', '))")
2145 : end if
2146 : ! LCOV_EXCL_STOP
2147 :
2148 12 : MeanAB_diff = abs(MeanAB - MeanAB_ref) / abs(MeanAB_ref)
2149 12 : assertionCurrent = all(MeanAB_diff <= tolerance)
2150 3 : assertion = assertion .and. assertionCurrent
2151 :
2152 : ! LCOV_EXCL_START
2153 : if (Test%isDebugMode .and. .not. assertionCurrent) then
2154 : write(Test%outputUnit,"(*(g0,:,', '))")
2155 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_ref ", MeanAB_ref
2156 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB ", MeanAB
2157 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_diff ", MeanAB_diff
2158 : write(Test%outputUnit,"(*(g0,:,', '))")
2159 : end if
2160 : ! LCOV_EXCL_STOP
2161 :
2162 3 : end function test_mergeMeanCovUpper_1
2163 :
2164 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2165 :
2166 3 : function test_mergeMeanCovUpperDense_1() result(assertion)
2167 :
2168 3 : use Constants_mod, only: IK, RK
2169 : implicit none
2170 :
2171 : integer(IK) :: i, j
2172 : logical :: assertion, assertionCurrent
2173 : real(RK) , parameter :: tolerance = 1.e-12_RK
2174 : integer(IK) , parameter :: nd = 3_IK
2175 :
2176 : integer(IK) , parameter :: npA = 5_IK
2177 : real(RK) :: MeanA(nd), CovMatUpperA(nd,nd)
2178 : real(RK) , parameter :: PointA(nd,npA) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
2179 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
2180 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
2181 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
2182 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
2183 : ], shape = shape(PointA) )
2184 :
2185 : integer(IK) , parameter :: npB = 10_IK
2186 : real(RK) :: MeanB(nd), CovMatUpperB(nd,nd)
2187 : real(RK) , parameter :: PointB(nd,npB) = reshape([ 0.706046088019609_RK, 0.031832846377421_RK, 0.276922984960890_RK &
2188 : , 0.046171390631154_RK, 0.097131781235848_RK, 0.823457828327293_RK &
2189 : , 0.694828622975817_RK, 0.317099480060861_RK, 0.950222048838355_RK &
2190 : , 0.034446080502909_RK, 0.438744359656398_RK, 0.381558457093008_RK &
2191 : , 0.765516788149002_RK, 0.795199901137063_RK, 0.186872604554379_RK &
2192 : , 1.706046088019609_RK, 2.031832846377421_RK, 3.276922984960890_RK &
2193 : , 1.046171390631154_RK, 2.097131781235848_RK, 3.823457828327293_RK &
2194 : , 1.694828622975817_RK, 2.317099480060861_RK, 3.950222048838355_RK &
2195 : , 1.034446080502909_RK, 2.438744359656398_RK, 3.381558457093008_RK &
2196 : , 1.765516788149002_RK, 2.795199901137063_RK, 3.186872604554379_RK &
2197 : ], shape = shape(PointB) )
2198 :
2199 :
2200 : integer(IK) , parameter :: npAB = npA + npB
2201 : real(RK) :: PointAB(nd,npAB)
2202 : real(RK) :: CovMatUpperAB_diff(nd,nd)
2203 : real(RK) :: MeanAB(nd), CovMatUpperAB(nd,nd)
2204 : real(RK) , parameter :: CovMatUpperAB_refMATLAB(nd,nd) = reshape([ 0.358269305364807_RK, 0.501078279929690_RK, 0.687067319924494_RK &
2205 : , 0.501078279929690_RK, 1.031956780716069_RK, 1.391311858397106_RK &
2206 : , 0.687067319924494_RK, 1.391311858397106_RK, 2.242785335243168_RK &
2207 : ], shape = shape(CovMatUpperAB_refMATLAB) )
2208 : real(RK) :: CovMatUpperAB_ref(nd,nd)
2209 : real(RK) :: MeanAB_diff(nd)
2210 : real(RK) :: MeanAB_ref(nd)
2211 :
2212 3 : assertion = .true.
2213 :
2214 3 : call getSamCovUpperMeanTrans(npA,nd,PointA,CovMatUpperA,MeanA)
2215 3 : call getSamCovUpperMeanTrans(npB,nd,PointB,CovMatUpperB,MeanB)
2216 :
2217 63 : PointAB(:,1:npA) = PointA
2218 123 : PointAB(:,npA+1:npAB) = PointB
2219 :
2220 3 : call getSamCovUpperMeanTrans(npAB,nd,PointAB,CovMatUpperAB_ref,MeanAB_ref)
2221 :
2222 3 : MeanAB = MeanB
2223 3 : CovMatUpperAB = CovMatUpperB
2224 : !do i = 1, 10000000-1
2225 3 : call mergeMeanCovUpperDense(nd,npA,MeanA,CovMatUpperA,npB,MeanAB,CovMatUpperAB)
2226 : !end do
2227 :
2228 3 : MeanAB = MeanB
2229 3 : CovMatUpperAB = CovMatUpperB
2230 3 : call mergeMeanCovUpperDense(nd,npA,MeanA,CovMatUpperA,npB,MeanAB,CovMatUpperAB)
2231 :
2232 3 : assertionCurrent = .true.
2233 12 : do j = 1, nd
2234 30 : do i = 1, j
2235 18 : CovMatUpperAB_diff(i,j) = abs(CovMatUpperAB(i,j) - CovMatUpperAB_ref(i,j)) / abs(CovMatUpperAB_ref(i,j))
2236 27 : assertionCurrent = assertionCurrent .and. CovMatUpperAB_diff(i,j) <= tolerance
2237 : end do
2238 : end do
2239 3 : assertion = assertion .and. assertionCurrent
2240 :
2241 : ! LCOV_EXCL_START
2242 : if (Test%isDebugMode .and. .not. assertionCurrent) then
2243 : write(Test%outputUnit,"(*(g0,:,', '))")
2244 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_refMATLAB ", ((CovMatUpperAB_refMATLAB(i,j), i=1,j), j=1,nd)
2245 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_ref ", ((CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2246 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_diff ", ((CovMatUpperAB_refMATLAB(i,j)-CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2247 : write(Test%outputUnit,"(*(g0,:,', '))")
2248 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_ref ", ((CovMatUpperAB_ref(i,j), i=1,j), j=1,nd)
2249 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB ", ((CovMatUpperAB(i,j), i=1,j), j=1,nd)
2250 : write(Test%outputUnit,"(*(g0,:,', '))") "CovMatUpperAB_diff ", ((CovMatUpperAB_diff(i,j), i=1,j), j=1,nd)
2251 : write(Test%outputUnit,"(*(g0,:,', '))")
2252 : end if
2253 : ! LCOV_EXCL_STOP
2254 :
2255 12 : MeanAB_diff = abs(MeanAB - MeanAB_ref) / abs(MeanAB_ref)
2256 12 : assertionCurrent = all(MeanAB_diff <= tolerance)
2257 3 : assertion = assertion .and. assertionCurrent
2258 :
2259 : ! LCOV_EXCL_START
2260 : if (Test%isDebugMode .and. .not. assertionCurrent) then
2261 : write(Test%outputUnit,"(*(g0,:,', '))")
2262 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_ref ", MeanAB_ref
2263 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB ", MeanAB
2264 : write(Test%outputUnit,"(*(g0,:,', '))") "MeanAB_diff ", MeanAB_diff
2265 : write(Test%outputUnit,"(*(g0,:,', '))")
2266 : end if
2267 : ! LCOV_EXCL_STOP
2268 :
2269 3 : end function test_mergeMeanCovUpperDense_1
2270 :
2271 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2272 :
2273 : ! \todo
2274 : ! What is the best method of testing for randomness?
2275 3 : function test_getRandGaus_1() result(assertion)
2276 3 : use Constants_mod, only: IK, RK
2277 : implicit none
2278 : logical :: assertion
2279 3 : real(RK) :: stdNormRnd
2280 3 : stdNormRnd = getRandGaus()
2281 3 : assertion = .true. !< There is really no easy testing of randomness. For now, we trust the giants.
2282 3 : end function test_getRandGaus_1
2283 :
2284 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2285 :
2286 : ! \todo
2287 : ! What is the best method of testing for randomness?
2288 3 : function test_getRandNorm_1() result(assertion)
2289 3 : use Constants_mod, only: IK, RK
2290 : implicit none
2291 : logical :: assertion
2292 3 : real(RK) :: normRnd
2293 3 : normRnd = getRandNorm(0._RK, 1._RK)
2294 3 : assertion = .true. !< There is really no easy testing of randomness. For now, we trust the giants.
2295 3 : end function test_getRandNorm_1
2296 :
2297 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2298 :
2299 : ! \todo
2300 : ! What is the best method of testing for randomness?
2301 3 : function test_getRandLogn_1() result(assertion)
2302 3 : use Constants_mod, only: IK, RK
2303 : implicit none
2304 : logical :: assertion
2305 3 : real(RK) :: lognRnd
2306 3 : lognRnd = getRandLogn(0._RK, 1._RK)
2307 3 : assertion = .true. !< There is really no easy testing of randomness. For now, we trust the giants.
2308 3 : end function test_getRandLogn_1
2309 :
2310 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2311 :
2312 : ! \todo
2313 : ! What is the best method of testing for randomness?
2314 3 : function test_getMVNDev_1() result(assertion)
2315 3 : use Constants_mod, only: IK, RK
2316 : implicit none
2317 : logical :: assertion
2318 : integer(IK) , parameter :: nd = 2_IK
2319 : real(RK) , parameter :: MeanVec(nd) = [ 1._RK, 2._RK ]
2320 : real(RK) , parameter :: CovMat(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2321 : real(RK) :: X(nd)
2322 3 : call getMVNDev(nd, MeanVec, CovMat, X)
2323 3 : assertion = .true. !< There is really no easy testing of randomness. For now, we trust the giants.
2324 3 : end function test_getMVNDev_1
2325 :
2326 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2327 :
2328 3 : function test_getMVUDev_1() result(assertion)
2329 3 : use Constants_mod, only: IK, RK
2330 : implicit none
2331 : logical :: assertion
2332 : integer(IK) , parameter :: nd = 2_IK
2333 : real(RK) , parameter :: MeanVec(nd) = [ 1._RK, 2._RK ]
2334 : real(RK) , parameter :: CovMat(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2335 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ +1.333333333333333_RK, -0.666666666666667_RK &
2336 : , -0.666666666666667_RK, +1.333333333333333_RK ] &
2337 : , shape = shape(InvCovMat) )
2338 : real(RK) :: X(nd), NormedPoint(nd)
2339 3 : call getMVUDev(nd, MeanVec, CovMat, X)
2340 9 : NormedPoint = X - MeanVec
2341 3 : assertion = isInsideEllipsoid(nd, NormedPoint, InvCovMat)
2342 3 : end function test_getMVUDev_1
2343 :
2344 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2345 :
2346 : ! \todo
2347 : ! What is the best method of testing for randomness?
2348 3 : function test_getRandMVN_1() result(assertion)
2349 3 : use Constants_mod, only: IK, RK
2350 : implicit none
2351 : logical :: assertion
2352 : integer(IK) , parameter :: nd = 2_IK
2353 : real(RK) , parameter :: MeanVec(nd) = [ 1._RK, 2._RK ]
2354 : real(RK) , parameter :: CovMat(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2355 : real(RK) , parameter :: CholeskyLower(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2356 : real(RK) , parameter :: Diagonal(nd) = [ 1._RK, 0.866025403784439_RK ]
2357 : real(RK) :: X(nd)
2358 3 : X = getRandMVN(nd, MeanVec, CholeskyLower, Diagonal)
2359 3 : assertion = .true. !< There is really no easy testing of randomness. For now, we trust the giants.
2360 3 : end function test_getRandMVN_1
2361 :
2362 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2363 :
2364 3 : function test_getRandMVU_1() result(assertion)
2365 3 : use Constants_mod, only: IK, RK
2366 : implicit none
2367 : logical :: assertion
2368 : integer(IK) , parameter :: nd = 2_IK
2369 : real(RK) , parameter :: MeanVec(nd) = [ 1._RK, 2._RK ]
2370 : real(RK) , parameter :: Diagonal(nd) = [ 1._RK, 0.866025403784439_RK ]
2371 : real(RK) , parameter :: CovMat(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2372 : real(RK) , parameter :: CholeskyLower(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2373 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ +1.333333333333333_RK, -0.666666666666667_RK &
2374 : , -0.666666666666667_RK, +1.333333333333333_RK ] &
2375 : , shape = shape(InvCovMat) )
2376 : real(RK) :: X(nd), NormedPoint(nd)
2377 3 : X = getRandMVU(nd, MeanVec, CholeskyLower, Diagonal)
2378 9 : NormedPoint = X - MeanVec
2379 3 : assertion = isInsideEllipsoid(nd, NormedPoint, InvCovMat)
2380 3 : end function test_getRandMVU_1
2381 :
2382 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2383 :
2384 3 : function test_isInsideEllipsoid_1() result(assertion)
2385 3 : use Constants_mod, only: IK, RK
2386 : implicit none
2387 : logical :: assertion
2388 : integer(IK) , parameter :: nd = 2_IK
2389 : real(RK) , parameter :: MeanVec(nd) = [ 1._RK, 2._RK ]
2390 : real(RK) , parameter :: Diagonal(nd) = [ 1._RK, 0.866025403784439_RK ]
2391 : real(RK) , parameter :: CovMat(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2392 : real(RK) , parameter :: CholeskyLower(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2393 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ +1.333333333333333_RK, -0.666666666666667_RK &
2394 : , -0.666666666666667_RK, +1.333333333333333_RK ] &
2395 : , shape = shape(InvCovMat) )
2396 : real(RK) :: X(nd), NormedPoint(nd)
2397 3 : X = getRandMVU(nd, MeanVec, CholeskyLower, Diagonal)
2398 9 : NormedPoint = X - MeanVec
2399 3 : assertion = isInsideEllipsoid(nd, NormedPoint, InvCovMat)
2400 3 : NormedPoint = [-1.e2_RK, 1.e2_RK]
2401 3 : assertion = assertion .and. .not. isInsideEllipsoid(nd, NormedPoint, InvCovMat)
2402 3 : end function test_isInsideEllipsoid_1
2403 :
2404 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2405 :
2406 3 : function test_getLogProbMVU_1() result(assertion)
2407 3 : use Constants_mod, only: IK, RK
2408 : implicit none
2409 : logical :: assertion
2410 : integer(IK) , parameter :: nd = 2_IK
2411 : real(RK) , parameter :: tolerance = 1.e-12_RK
2412 : real(RK) , parameter :: logSqrtDetCovMat = log(1._RK)
2413 : real(RK) , parameter :: logProbMVU_ref = -1.144729885849400_RK
2414 3 : real(RK) :: difference
2415 3 : real(RK) :: logProbMVU
2416 3 : logProbMVU = getLogProbMVU(nd, logSqrtDetCovMat)
2417 3 : difference = abs( (logProbMVU - logProbMVU_ref) / logProbMVU_ref)
2418 3 : assertion = difference < tolerance
2419 : ! LCOV_EXCL_START
2420 : if (Test%isDebugMode .and. .not. assertion) then
2421 : write(Test%outputUnit,"(*(g0,:,', '))")
2422 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMVU_ref ", logProbMVU_ref
2423 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbMVU ", logProbMVU
2424 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
2425 : write(Test%outputUnit,"(*(g0,:,', '))")
2426 : end if
2427 : ! LCOV_EXCL_STOP
2428 3 : end function test_getLogProbMVU_1
2429 :
2430 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2431 :
2432 3 : function test_getRandPointOnEllipsoid_1() result(assertion)
2433 3 : use Constants_mod, only: IK, RK
2434 : implicit none
2435 : logical :: assertion
2436 : integer(IK) , parameter :: nd = 2_IK
2437 : real(RK) , parameter :: tolerance = 1.e-12_RK
2438 : real(RK) , parameter :: MeanVec(nd) = [ 1._RK, 2._RK ]
2439 : real(RK) , parameter :: Diagonal(nd) = [ 1._RK, 0.866025403784439_RK ]
2440 : real(RK) , parameter :: CovMat(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2441 : real(RK) , parameter :: CholeskyLower(nd,nd) = reshape( [ 1._RK, 0.5_RK, 0.5_RK, 1._RK ], shape = shape(CovMat) )
2442 : real(RK) , parameter :: InvCovMat(nd,nd) = reshape( [ +1.333333333333333_RK, -0.666666666666667_RK &
2443 : , -0.666666666666667_RK, +1.333333333333333_RK ] &
2444 : , shape = shape(InvCovMat) )
2445 : real(RK) :: X(nd), NormedPoint(nd)
2446 3 : X = getRandPointOnEllipsoid(nd,MeanVec,CholeskyLower,Diagonal)
2447 9 : NormedPoint = X - MeanVec
2448 9 : assertion = dot_product(NormedPoint,matmul(InvCovMat,NormedPoint)) - 1._RK < tolerance
2449 3 : NormedPoint = [-1.e2_RK, 1.e2_RK]
2450 3 : assertion = assertion .and. .not. isInsideEllipsoid(nd, NormedPoint, InvCovMat)
2451 : ! LCOV_EXCL_START
2452 : if (Test%isDebugMode .and. .not. assertion) then
2453 : write(Test%outputUnit,"(*(g0,:,', '))")
2454 : write(Test%outputUnit,"(*(g0,:,', '))") "RandPointOnEllipsoid ", X
2455 : write(Test%outputUnit,"(*(g0,:,', '))") "distance from center ", dot_product(NormedPoint,matmul(InvCovMat,NormedPoint))
2456 : write(Test%outputUnit,"(*(g0,:,', '))") "expected distance ", 1._RK
2457 : write(Test%outputUnit,"(*(g0,:,', '))")
2458 : end if
2459 : ! LCOV_EXCL_STOP
2460 3 : end function test_getRandPointOnEllipsoid_1
2461 :
2462 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2463 :
2464 3 : function test_getLogProbLognSP_1() result(assertion)
2465 :
2466 3 : use Constants_mod, only: IK, RK
2467 : implicit none
2468 : logical :: assertion
2469 : integer(IK) , parameter :: nd = 2_IK
2470 : real(RK) , parameter :: tolerance = 1.e-12_RK
2471 : real(RK) , parameter :: logMean = 2._RK
2472 : real(RK) , parameter :: logPoint = log(5._RK)
2473 : real(RK) , parameter :: inverseVariance = 1.e-2_RK
2474 : real(RK) , parameter :: logSqrtInverseVariance = log(sqrt(inverseVariance))
2475 : real(RK) , parameter :: logProbLogn_ref = -4.831724232354038_RK
2476 3 : real(RK) :: difference
2477 3 : real(RK) :: logProbLogn
2478 :
2479 3 : logProbLogn = getLogProbLogn(logMean = logMean, inverseVariance = inverseVariance, logSqrtInverseVariance = logSqrtInverseVariance, logPoint = logPoint)
2480 3 : difference = abs( (logProbLogn - logProbLogn_ref) / logProbLogn_ref )
2481 3 : assertion = difference < tolerance
2482 :
2483 3 : logProbLogn = getLogProbLognorm(logMean = logMean, inverseVariance = inverseVariance, logSqrtInverseVariance = logSqrtInverseVariance, logPoint = logPoint)
2484 3 : difference = abs( (logProbLogn - logProbLogn_ref) / logProbLogn_ref )
2485 3 : assertion = assertion .and. difference < tolerance
2486 :
2487 : ! LCOV_EXCL_START
2488 : if (Test%isDebugMode .and. .not. assertion) then
2489 : write(Test%outputUnit,"(*(g0,:,', '))")
2490 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbLogn_ref ", logProbLogn_ref
2491 : write(Test%outputUnit,"(*(g0,:,', '))") "logProbLogn ", logProbLogn
2492 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
2493 : write(Test%outputUnit,"(*(g0,:,', '))")
2494 : end if
2495 : ! LCOV_EXCL_STOP
2496 :
2497 3 : end function test_getLogProbLognSP_1
2498 :
2499 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2500 :
2501 3 : function test_getLogProbLognMP_1() result(assertion)
2502 :
2503 3 : use Constants_mod, only: IK, RK
2504 : implicit none
2505 : logical :: assertion
2506 : integer(IK) , parameter :: nd = 2_IK, np = 3_IK
2507 : real(RK) , parameter :: tolerance = 1.e-12_RK
2508 : real(RK) , parameter :: logMean = 2._RK
2509 : real(RK) , parameter :: LogPoint(np) = log([5._RK, 6._RK, 7._RK])
2510 : real(RK) , parameter :: inverseVariance = 1.e-2_RK
2511 : real(RK) , parameter :: logSqrtInverseVariance = log(sqrt(inverseVariance))
2512 : real(RK) , parameter :: LogProbLogn_ref(np) = [ -4.831724232354038_RK, -5.013499916020054_RK, -5.167448403813908_RK ]
2513 : real(RK) :: LogProbLogn(np)
2514 : real(RK) :: difference(np)
2515 :
2516 3 : LogProbLogn = getLogProbLogn(np = np, logMean = logMean, inverseVariance = inverseVariance, logSqrtInverseVariance = logSqrtInverseVariance, logPoint = LogPoint)
2517 12 : difference = abs( (LogProbLogn - LogProbLogn_ref) / LogProbLogn_ref )
2518 12 : assertion = all(difference < tolerance)
2519 :
2520 3 : LogProbLogn = getLogProbLognorm(np = np, logMean = logMean, inverseVariance = inverseVariance, logSqrtInverseVariance = logSqrtInverseVariance, logPoint = LogPoint)
2521 12 : difference = abs( (LogProbLogn - LogProbLogn_ref) / LogProbLogn_ref )
2522 12 : assertion = assertion .and. all(difference < tolerance)
2523 :
2524 : ! LCOV_EXCL_START
2525 : if (Test%isDebugMode .and. .not. assertion) then
2526 : write(Test%outputUnit,"(*(g0,:,', '))")
2527 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbLogn_ref ", LogProbLogn_ref
2528 : write(Test%outputUnit,"(*(g0,:,', '))") "LogProbLogn ", LogProbLogn
2529 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
2530 : write(Test%outputUnit,"(*(g0,:,', '))")
2531 : end if
2532 : ! LCOV_EXCL_STOP
2533 :
2534 3 : end function test_getLogProbLognMP_1
2535 :
2536 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2537 :
2538 3 : function test_getRandRealLecuyer_1() result(assertion)
2539 3 : use Constants_mod, only: IK, RK
2540 : implicit none
2541 : logical :: assertion
2542 : integer(IK) , parameter :: np = 100
2543 : integer(IK) :: idum = 3333, i
2544 : real(RK) :: RandRealLecuyer(np)
2545 3 : assertion = .true.
2546 303 : do i = 1, np
2547 300 : RandRealLecuyer(i) = getRandRealLecuyer(idum)
2548 300 : assertion = assertion .and. RandRealLecuyer(i) <= 1._RK .and. RandRealLecuyer(i) >= 0._RK
2549 : ! LCOV_EXCL_START
2550 : if (Test%isDebugMode .and. .not. assertion) then
2551 : write(Test%outputUnit,"(*(g0,:,' '))")
2552 : write(Test%outputUnit,"(*(g0,:,' '))") "RandRealLecuyer(",i,") =", RandRealLecuyer(i)
2553 : write(Test%outputUnit,"(*(g0,:,' '))")
2554 : end if
2555 : ! LCOV_EXCL_STOP
2556 : end do
2557 3 : end function test_getRandRealLecuyer_1
2558 :
2559 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2560 :
2561 3 : function test_getRandIntLecuyer_1() result(assertion)
2562 3 : use Constants_mod, only: IK, RK
2563 : implicit none
2564 : logical :: assertion
2565 : integer(IK) , parameter :: lowerBound = -2, upperBound = 9
2566 : integer(IK) , parameter :: np = 100
2567 : integer(IK) :: idum = 3333, i
2568 : integer(IK) :: RandIntLecuyer(np)
2569 3 : assertion = .true.
2570 303 : do i = 1, np
2571 300 : RandIntLecuyer(i) = getRandIntLecuyer(lowerBound,upperBound,idum)
2572 300 : assertion = assertion .and. RandIntLecuyer(i) <= upperBound .and. RandIntLecuyer(i) >= lowerBound
2573 : ! LCOV_EXCL_START
2574 : if (Test%isDebugMode .and. .not. assertion) then
2575 : write(Test%outputUnit,"(*(g0,:,' '))")
2576 : write(Test%outputUnit,"(*(g0,:,' '))") "RandIntLecuyer(",i,") =", RandIntLecuyer(i)
2577 : write(Test%outputUnit,"(*(g0,:,' '))")
2578 : end if
2579 : ! LCOV_EXCL_STOP
2580 : end do
2581 3 : end function test_getRandIntLecuyer_1
2582 :
2583 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2584 :
2585 3 : function test_getRandUniform_1() result(assertion)
2586 3 : use Constants_mod, only: IK, RK
2587 : implicit none
2588 : logical :: assertion
2589 : integer(IK) , parameter :: np = 100
2590 : real(RK) , parameter :: lowerBound = -2._RK, upperBound = 9._RK
2591 : integer(IK) :: i
2592 : real(RK) :: RandUniform(np)
2593 3 : assertion = .true.
2594 303 : do i = 1, np
2595 300 : RandUniform(i) = getRandUniform(lowerBound,upperBound)
2596 300 : assertion = assertion .and. RandUniform(i) <= upperBound .and. RandUniform(i) >= lowerBound
2597 : ! LCOV_EXCL_START
2598 : if (Test%isDebugMode .and. .not. assertion) then
2599 : write(Test%outputUnit,"(*(g0,:,' '))")
2600 : write(Test%outputUnit,"(*(g0,:,' '))") "RandUniform(",i,") =", RandUniform(i)
2601 : write(Test%outputUnit,"(*(g0,:,' '))")
2602 : end if
2603 : ! LCOV_EXCL_STOP
2604 : end do
2605 3 : end function test_getRandUniform_1
2606 :
2607 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2608 :
2609 3 : function test_getRandGamma_1() result(assertion)
2610 3 : use Constants_mod, only: IK, RK, HUGE_RK
2611 : implicit none
2612 : logical :: assertion
2613 : integer(IK) , parameter :: np = 100
2614 : real(RK) , parameter :: alpha = 2._RK
2615 : real(RK) , parameter :: lowerBound = 0._RK, upperBound = HUGE_RK
2616 : integer(IK) :: i
2617 : real(RK) :: RandGamma(np)
2618 3 : assertion = .true.
2619 303 : do i = 1, np
2620 300 : RandGamma(i) = getRandGamma(alpha)
2621 300 : assertion = assertion .and. RandGamma(i) <= upperBound .and. RandGamma(i) >= lowerBound
2622 : ! LCOV_EXCL_START
2623 : if (Test%isDebugMode .and. .not. assertion) then
2624 : write(Test%outputUnit,"(*(g0,:,' '))")
2625 : write(Test%outputUnit,"(*(g0,:,' '))") "RandGamma(",i,") =", RandGamma(i)
2626 : write(Test%outputUnit,"(*(g0,:,' '))")
2627 : end if
2628 : ! LCOV_EXCL_STOP
2629 : end do
2630 3 : assertion = assertion .and. getRandGamma(alpha=-1._RK) < 0._RK
2631 3 : end function test_getRandGamma_1
2632 :
2633 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2634 :
2635 3 : function test_getRandGammaIntShape_1() result(assertion)
2636 3 : use Constants_mod, only: IK, RK, HUGE_RK
2637 : implicit none
2638 : logical :: assertion
2639 : integer(IK) , parameter :: np = 100
2640 : integer(IK) , parameter :: alpha = 2_IK
2641 : real(RK) , parameter :: lowerBound = 0._RK, upperBound = HUGE_RK
2642 : integer(IK) :: i
2643 : real(RK) :: RandGamma(np)
2644 3 : assertion = .true.
2645 303 : do i = 1, np
2646 300 : RandGamma(i) = getRandGammaIntShape(alpha)
2647 300 : assertion = assertion .and. RandGamma(i) <= upperBound .and. RandGamma(i) >= lowerBound
2648 : ! LCOV_EXCL_START
2649 : if (Test%isDebugMode .and. .not. assertion) then
2650 : write(Test%outputUnit,"(*(g0,:,' '))")
2651 : write(Test%outputUnit,"(*(g0,:,' '))") "RandGamma(",i,") =", RandGamma(i)
2652 : write(Test%outputUnit,"(*(g0,:,' '))")
2653 : end if
2654 : ! LCOV_EXCL_STOP
2655 : end do
2656 3 : assertion = assertion .and. getRandGammaIntShape(alpha=-1_IK) < 0._RK
2657 3 : end function test_getRandGammaIntShape_1
2658 :
2659 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2660 :
2661 3 : function test_getRandBeta_1() result(assertion)
2662 3 : use Constants_mod, only: IK, RK
2663 : implicit none
2664 : logical :: assertion
2665 : integer(IK) , parameter :: np = 100
2666 : real(RK) , parameter :: alpha = 2._RK, beta = 3._RK
2667 : real(RK) , parameter :: lowerBound = 0._RK, upperBound = 1._RK
2668 : integer(IK) :: i
2669 : real(RK) :: RandBeta(np)
2670 3 : assertion = .true.
2671 303 : do i = 1, np
2672 300 : RandBeta(i) = getRandBeta(alpha, beta)
2673 300 : assertion = assertion .and. RandBeta(i) <= upperBound .and. RandBeta(i) >= lowerBound
2674 : ! LCOV_EXCL_START
2675 : if (Test%isDebugMode .and. .not. assertion) then
2676 : write(Test%outputUnit,"(*(g0,:,' '))")
2677 : write(Test%outputUnit,"(*(g0,:,' '))") "RandBeta(",i,") =", RandBeta(i)
2678 : write(Test%outputUnit,"(*(g0,:,' '))")
2679 : end if
2680 : ! LCOV_EXCL_STOP
2681 : end do
2682 3 : assertion = assertion .and. getRandBeta(alpha=-1._RK, beta=+2._RK) < 0._RK
2683 3 : assertion = assertion .and. getRandBeta(alpha=+1._RK, beta=-2._RK) < 0._RK
2684 3 : assertion = assertion .and. getRandBeta(alpha=-1._RK, beta=-2._RK) < 0._RK
2685 3 : end function test_getRandBeta_1
2686 :
2687 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2688 :
2689 3 : function test_getRandExpWithInvMean_1() result(assertion)
2690 3 : use Constants_mod, only: IK, RK, HUGE_RK
2691 : implicit none
2692 : logical :: assertion
2693 : integer(IK) , parameter :: np = 100
2694 : real(RK) , parameter :: invMean = 0.1_RK
2695 : real(RK) , parameter :: lowerBound = 0._RK, upperBound = HUGE_RK
2696 : integer(IK) :: i
2697 : real(RK) :: RandExp(np)
2698 3 : assertion = .true.
2699 303 : do i = 1, np
2700 300 : RandExp(i) = getRandExpWithInvMean(invMean)
2701 300 : assertion = assertion .and. RandExp(i) <= upperBound .and. RandExp(i) >= lowerBound
2702 : ! LCOV_EXCL_START
2703 : if (Test%isDebugMode .and. .not. assertion) then
2704 : write(Test%outputUnit,"(*(g0,:,' '))")
2705 : write(Test%outputUnit,"(*(g0,:,' '))") "RandExp(",i,") =", RandExp(i)
2706 : write(Test%outputUnit,"(*(g0,:,' '))")
2707 : end if
2708 : ! LCOV_EXCL_STOP
2709 : end do
2710 3 : end function test_getRandExpWithInvMean_1
2711 :
2712 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2713 :
2714 3 : function test_getRandExp_1() result(assertion)
2715 3 : use Constants_mod, only: IK, RK, HUGE_RK
2716 : implicit none
2717 : logical :: assertion
2718 : integer(IK) , parameter :: np = 100
2719 : real(RK) , parameter :: lowerBound = 0._RK, upperBound = HUGE_RK
2720 : integer(IK) :: i
2721 : real(RK) :: RandExp(np)
2722 3 : assertion = .true.
2723 303 : do i = 1, np
2724 300 : RandExp(i) = getRandExp()
2725 300 : assertion = assertion .and. RandExp(i) <= upperBound .and. RandExp(i) >= lowerBound
2726 : ! LCOV_EXCL_START
2727 : if (Test%isDebugMode .and. .not. assertion) then
2728 : write(Test%outputUnit,"(*(g0,:,' '))")
2729 : write(Test%outputUnit,"(*(g0,:,' '))") "RandExp(",i,") =", RandExp(i)
2730 : write(Test%outputUnit,"(*(g0,:,' '))")
2731 : end if
2732 : ! LCOV_EXCL_STOP
2733 : end do
2734 3 : end function test_getRandExp_1
2735 :
2736 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2737 :
2738 3 : function test_getRandCorMat_1() result(assertion)
2739 3 : use Constants_mod, only: IK, RK, HUGE_RK
2740 : use Matrix_mod, only: isPosDef
2741 : implicit none
2742 : logical :: assertion
2743 : integer(IK) , parameter :: nd = 2_IK
2744 : integer(IK) , parameter :: np = 100_IK
2745 : real(RK) , parameter :: eta = 5._RK
2746 : real(RK) , parameter :: lowerBound = -1._RK, upperBound = 1._RK
2747 : integer(IK) :: i
2748 : real(RK) :: RandCorMat(nd,nd)
2749 3 : assertion = .true.
2750 303 : do i = 1, np
2751 300 : RandCorMat = getRandCorMat(nd,eta)
2752 3900 : assertion = assertion .and. all(RandCorMat <= upperBound) .and. all(RandCorMat >= lowerBound)
2753 300 : assertion = assertion .and. isPosDef(nd,RandCorMat)
2754 : ! LCOV_EXCL_START
2755 : if (Test%isDebugMode .and. .not. assertion) then
2756 : write(Test%outputUnit,"(*(g0,:,' '))")
2757 : write(Test%outputUnit,"(*(g0,:,' '))") "RandCorMat(:,:,) =", RandCorMat
2758 : write(Test%outputUnit,"(*(g0,:,' '))")
2759 : end if
2760 : ! LCOV_EXCL_STOP
2761 : end do
2762 : ! do not uncomment this as GNU Fortran 9.1 crashes on this in debug mode with error message:
2763 : ! Fortran runtime error: Dimension 1 of array 'randcormat' has extent 0 instead of 2...
2764 : !RandCorMat = getRandCorMat(nd = 0_IK, eta = 1._RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2765 : !RandCorMat = getRandCorMat(nd = -1_IK, eta = +2._RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2766 : !RandCorMat = getRandCorMat(nd = +1_IK, eta = +0._RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2767 3 : RandCorMat = getRandCorMat(nd = +2_IK, eta = +0._RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2768 3 : RandCorMat = getRandCorMat(nd = +2_IK, eta = -1._RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2769 3 : RandCorMat = getRandCorMat(nd = +2_IK, eta = -2._RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2770 3 : end function test_getRandCorMat_1
2771 :
2772 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2773 :
2774 3 : function test_getRandCorMatRejection_1() result(assertion)
2775 3 : use Constants_mod, only: IK, RK, HUGE_RK
2776 : use Matrix_mod, only: isPosDef
2777 : implicit none
2778 : logical :: assertion, assertionCurrent
2779 : integer(IK) , parameter :: nd = 2_IK
2780 : integer(IK) , parameter :: np = 100_IK
2781 : real(RK) , parameter :: minRho = -.3_RK, maxRho = .6_RK
2782 : real(RK) , parameter :: lowerBound = -1._RK, upperBound = 1._RK
2783 : integer(IK) :: i, j, k
2784 : real(RK) :: RandCorMat(nd,nd)
2785 3 : assertion = .true.
2786 303 : do i = 1, np
2787 300 : RandCorMat = getRandCorMatRejection(nd,minRho,maxRho)
2788 3900 : assertion = assertion .and. all(RandCorMat <= upperBound) .and. all(RandCorMat >= lowerBound)
2789 900 : do j = 1, nd
2790 2100 : do k = 1, nd
2791 1200 : if (j==k) then
2792 600 : assertionCurrent = RandCorMat(j,k) == upperBound
2793 : else
2794 600 : assertionCurrent = RandCorMat(j,k) <= maxRho .and. RandCorMat(j,k) >= minRho
2795 : end if
2796 1800 : assertion = assertion .and. assertionCurrent
2797 : end do
2798 : end do
2799 300 : assertionCurrent = isPosDef(nd,RandCorMat)
2800 : ! LCOV_EXCL_START
2801 : if (Test%isDebugMode .and. .not. assertion) then
2802 : write(Test%outputUnit,"(*(g0,:,' '))")
2803 : write(Test%outputUnit,"(*(g0,:,' '))") "RandCorMat(:,:,) =", RandCorMat
2804 : write(Test%outputUnit,"(*(g0,:,' '))")
2805 : end if
2806 : ! LCOV_EXCL_STOP
2807 303 : assertion = assertion .and. assertionCurrent
2808 : end do
2809 : ! do not uncomment this as GNU Fortran 9.1 crashes on this in debug mode with error message:
2810 : ! Fortran runtime error: Dimension 1 of array 'randcormat' has extent 1 instead of 2...
2811 : !RandCorMat = getRandCorMatRejection(nd = +1_IK, minRho = +.2_RK, maxRho = -.5_RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2812 3 : RandCorMat = getRandCorMatRejection(nd = +2_IK, minRho = -.2_RK, maxRho = -.5_RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2813 3 : RandCorMat = getRandCorMatRejection(nd = +2_IK, minRho = +.2_RK, maxRho = -.5_RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2814 : !RandCorMat = getRandCorMatRejection(nd = +0_IK, minRho = +.2_RK, maxRho = +.5_RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2815 : !RandCorMat = getRandCorMatRejection(nd = -3_IK, minRho = +.2_RK, maxRho = +.5_RK); assertion = assertion .and. RandCorMat(1,1) < 0._RK
2816 3 : end function test_getRandCorMatRejection_1
2817 :
2818 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2819 :
2820 3 : function test_getCorMatUpperFromCovMatUpper_1() result(assertion)
2821 3 : use Constants_mod, only: IK, RK
2822 : implicit none
2823 : logical :: assertion
2824 : integer(IK) , parameter :: nd = 2_IK
2825 : real(RK) , parameter :: tolerance = 1.e-12_RK
2826 : real(RK) , parameter :: CovMatUpper (nd,nd) = reshape( [ 2._RK, 0.00_RK, 0.50_RK, 2._RK ], shape = shape(CovMatUpper) )
2827 : real(RK) , parameter :: CorMatUpper_ref(nd,nd) = reshape( [ 1._RK, 0.00_RK, 0.25_RK, 1._RK ], shape = shape(CorMatUpper_ref) )
2828 : integer(IK) :: i, j
2829 : real(RK) :: CorMatUpper(nd,nd)
2830 : real(RK) :: Difference(nd,nd)
2831 :
2832 3 : CorMatUpper = getCorMatUpperFromCovMatUpper(nd, CovMatUpper)
2833 :
2834 3 : assertion = .true.
2835 9 : do j = 1, nd
2836 18 : do i = 1, j
2837 9 : Difference(i,j) = abs(CorMatUpper(i,j) - CorMatUpper_ref(i,j)) / abs(CorMatUpper_ref(i,j))
2838 15 : assertion = assertion .and. Difference(i,j) < tolerance
2839 : end do
2840 : end do
2841 :
2842 : ! LCOV_EXCL_START
2843 : if (Test%isDebugMode .and. .not. assertion) then
2844 : write(Test%outputUnit,"(*(g0,:,' '))")
2845 : write(Test%outputUnit,"(*(g0,:,' '))") "CorMatUpper_ref =", CorMatUpper_ref
2846 : write(Test%outputUnit,"(*(g0,:,' '))") "CorMatUpper =", CorMatUpper
2847 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
2848 : write(Test%outputUnit,"(*(g0,:,' '))")
2849 : end if
2850 : ! LCOV_EXCL_STOP
2851 :
2852 3 : end function test_getCorMatUpperFromCovMatUpper_1
2853 :
2854 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2855 :
2856 3 : function test_getCovMatUpperFromCorMatUpper_1() result(assertion)
2857 3 : use Constants_mod, only: IK, RK
2858 : implicit none
2859 : logical :: assertion
2860 : integer(IK) , parameter :: nd = 2_IK
2861 : real(RK) , parameter :: tolerance = 1.e-12_RK
2862 : real(RK) , parameter :: CovMatUpper_ref(nd,nd) = reshape( [ 2._RK, 0.00_RK, 0.50_RK, 2._RK ], shape = shape(CovMatUpper_ref) )
2863 : real(RK) , parameter :: CorMatUpper(nd,nd) = reshape( [ 1._RK, 0.00_RK, 0.25_RK, 1._RK ], shape = shape(CorMatUpper) )
2864 : real(RK) , parameter :: StdVec(nd) = [ 1.414213562373095_RK, 1.414213562373095_RK ]
2865 : integer(IK) :: i, j
2866 : real(RK) :: CovMatUpper(nd,nd)
2867 : real(RK) :: Difference(nd,nd)
2868 :
2869 3 : CovMatUpper = getCovMatUpperFromCorMatUpper(nd, StdVec, CorMatUpper)
2870 :
2871 3 : assertion = .true.
2872 9 : do j = 1, nd
2873 18 : do i = 1, j
2874 9 : Difference(i,j) = abs(CovMatUpper(i,j) - CovMatUpper_ref(i,j)) / abs(CovMatUpper_ref(i,j))
2875 15 : assertion = assertion .and. Difference(i,j) < tolerance
2876 : end do
2877 : end do
2878 :
2879 : ! LCOV_EXCL_START
2880 : if (Test%isDebugMode .and. .not. assertion) then
2881 : write(Test%outputUnit,"(*(g0,:,' '))")
2882 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatUpper_ref =", CovMatUpper_ref
2883 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatUpper =", CovMatUpper
2884 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
2885 : write(Test%outputUnit,"(*(g0,:,' '))")
2886 : end if
2887 : ! LCOV_EXCL_STOP
2888 :
2889 3 : end function test_getCovMatUpperFromCorMatUpper_1
2890 :
2891 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2892 :
2893 3 : function test_getCovMatUpperFromCorMatLower_1() result(assertion)
2894 3 : use Constants_mod, only: IK, RK
2895 : implicit none
2896 : logical :: assertion
2897 : integer(IK) , parameter :: nd = 2_IK
2898 : real(RK) , parameter :: tolerance = 1.e-12_RK
2899 : real(RK) , parameter :: CovMatUpper_ref(nd,nd) = reshape( [ 2._RK, 0.00_RK, 0.50_RK, 2._RK ], shape = shape(CovMatUpper_ref) )
2900 : real(RK) , parameter :: CorMatLower(nd,nd) = transpose(reshape( [ 1._RK, 0.00_RK, 0.25_RK, 1._RK ], shape = shape(CorMatLower) ))
2901 : real(RK) , parameter :: StdVec(nd) = [ 1.414213562373095_RK, 1.414213562373095_RK ]
2902 : integer(IK) :: i, j
2903 : real(RK) :: CovMatUpper(nd,nd)
2904 : real(RK) :: Difference(nd,nd)
2905 :
2906 3 : CovMatUpper = getCovMatUpperFromCorMatLower(nd, StdVec, CorMatLower)
2907 :
2908 3 : assertion = .true.
2909 9 : do j = 1, nd
2910 18 : do i = 1, j
2911 9 : Difference(i,j) = abs(CovMatUpper(i,j) - CovMatUpper_ref(i,j)) / abs(CovMatUpper_ref(i,j))
2912 15 : assertion = assertion .and. Difference(i,j) < tolerance
2913 : end do
2914 : end do
2915 :
2916 3 : if (Test%isDebugMode .and. .not. assertion) then
2917 : ! LCOV_EXCL_START
2918 : write(Test%outputUnit,"(*(g0,:,' '))")
2919 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatUpper_ref =", CovMatUpper_ref
2920 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatUpper =", CovMatUpper
2921 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
2922 : write(Test%outputUnit,"(*(g0,:,' '))")
2923 : end if
2924 : ! LCOV_EXCL_STOP
2925 :
2926 3 : end function test_getCovMatUpperFromCorMatLower_1
2927 :
2928 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2929 :
2930 3 : function test_getCovMatLowerFromCorMatUpper_1() result(assertion)
2931 3 : use Constants_mod, only: IK, RK
2932 : implicit none
2933 : logical :: assertion
2934 : integer(IK) , parameter :: nd = 2_IK
2935 : real(RK) , parameter :: tolerance = 1.e-12_RK
2936 : real(RK) , parameter :: CovMatLower_ref(nd,nd) = transpose( reshape( [ 2._RK, 0.00_RK, 0.50_RK, 2._RK ], shape = shape(CovMatLower_ref) ) )
2937 : real(RK) , parameter :: CorMatUpper(nd,nd) = reshape( [ 1._RK, 0.00_RK, 0.25_RK, 1._RK ], shape = shape(CorMatUpper) )
2938 : real(RK) , parameter :: StdVec(nd) = [ 1.414213562373095_RK, 1.414213562373095_RK ]
2939 : integer(IK) :: i, j
2940 : real(RK) :: CovMatLower(nd,nd)
2941 : real(RK) :: Difference(nd,nd)
2942 :
2943 3 : CovMatLower = getCovMatLowerFromCorMatUpper(nd, StdVec, CorMatUpper)
2944 :
2945 3 : assertion = .true.
2946 9 : do j = 1, nd
2947 18 : do i = j, nd
2948 9 : Difference(i,j) = abs(CovMatLower(i,j) - CovMatLower_ref(i,j)) / abs(CovMatLower_ref(i,j))
2949 15 : assertion = assertion .and. Difference(i,j) < tolerance
2950 : end do
2951 : end do
2952 :
2953 3 : if (Test%isDebugMode .and. .not. assertion) then
2954 : ! LCOV_EXCL_START
2955 : write(Test%outputUnit,"(*(g0,:,' '))")
2956 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatLower_ref =", CovMatLower_ref
2957 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatLower =", CovMatLower
2958 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
2959 : write(Test%outputUnit,"(*(g0,:,' '))")
2960 : end if
2961 : ! LCOV_EXCL_STOP
2962 :
2963 3 : end function test_getCovMatLowerFromCorMatUpper_1
2964 :
2965 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2966 :
2967 3 : function test_getCovMatLowerFromCorMatLower_1() result(assertion)
2968 3 : use Constants_mod, only: IK, RK
2969 : implicit none
2970 : logical :: assertion
2971 : integer(IK) , parameter :: nd = 2_IK
2972 : real(RK) , parameter :: tolerance = 1.e-12_RK
2973 : real(RK) , parameter :: CovMatLower_ref(nd,nd) = transpose( reshape( [ 2._RK, 0.00_RK, 0.50_RK, 2._RK ], shape = shape(CovMatLower_ref) ) )
2974 : real(RK) , parameter :: CorMatLower(nd,nd) = transpose( reshape( [ 1._RK, 0.00_RK, 0.25_RK, 1._RK ], shape = shape(CorMatLower) ) )
2975 : real(RK) , parameter :: StdVec(nd) = [ 1.414213562373095_RK, 1.414213562373095_RK ]
2976 : integer(IK) :: i, j
2977 : real(RK) :: CovMatLower(nd,nd)
2978 : real(RK) :: Difference(nd,nd)
2979 :
2980 3 : CovMatLower = getCovMatLowerFromCorMatLower(nd, StdVec, CorMatLower)
2981 :
2982 3 : assertion = .true.
2983 9 : do j = 1, nd
2984 18 : do i = j, nd
2985 9 : Difference(i,j) = abs(CovMatLower(i,j) - CovMatLower_ref(i,j)) / abs(CovMatLower_ref(i,j))
2986 15 : assertion = assertion .and. Difference(i,j) < tolerance
2987 : end do
2988 : end do
2989 :
2990 3 : if (Test%isDebugMode .and. .not. assertion) then
2991 : ! LCOV_EXCL_START
2992 : write(Test%outputUnit,"(*(g0,:,' '))")
2993 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatLower_ref =", CovMatLower_ref
2994 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMatLower =", CovMatLower
2995 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
2996 : write(Test%outputUnit,"(*(g0,:,' '))")
2997 : end if
2998 : ! LCOV_EXCL_STOP
2999 :
3000 3 : end function test_getCovMatLowerFromCorMatLower_1
3001 :
3002 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3003 :
3004 3 : function test_getCovMatFromCorMatUpper_1() result(assertion)
3005 3 : use Constants_mod, only: IK, RK
3006 : implicit none
3007 : logical :: assertion
3008 : integer(IK) , parameter :: nd = 2_IK
3009 : real(RK) , parameter :: tolerance = 1.e-12_RK
3010 : real(RK) , parameter :: CovMat_ref(nd,nd) = reshape( [ 2._RK, 0.50_RK, 0.50_RK, 2._RK ], shape = shape(CovMat_ref) )
3011 : real(RK) , parameter :: CorMatUpper(nd,nd) = reshape( [ 1._RK, 0.00_RK, 0.25_RK, 1._RK ], shape = shape(CorMatUpper) )
3012 : real(RK) , parameter :: StdVec(nd) = [ 1.414213562373095_RK, 1.414213562373095_RK ]
3013 : integer(IK) :: i, j
3014 : real(RK) :: CovMat(nd,nd)
3015 : real(RK) :: Difference(nd,nd)
3016 :
3017 3 : CovMat = getCovMatFromCorMatUpper(nd, StdVec, CorMatUpper)
3018 :
3019 3 : assertion = .true.
3020 9 : do j = 1, nd
3021 21 : do i = 1, nd
3022 12 : Difference(i,j) = abs(CovMat(i,j) - CovMat_ref(i,j)) / abs(CovMat_ref(i,j))
3023 18 : assertion = assertion .and. Difference(i,j) < tolerance
3024 : end do
3025 : end do
3026 :
3027 3 : if (Test%isDebugMode .and. .not. assertion) then
3028 : ! LCOV_EXCL_START
3029 : write(Test%outputUnit,"(*(g0,:,' '))")
3030 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMat_ref =", CovMat_ref
3031 : write(Test%outputUnit,"(*(g0,:,' '))") "CovMat =", CovMat
3032 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
3033 : write(Test%outputUnit,"(*(g0,:,' '))")
3034 : end if
3035 : ! LCOV_EXCL_STOP
3036 :
3037 3 : end function test_getCovMatFromCorMatUpper_1
3038 :
3039 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3040 :
3041 3 : function test_getLogProbGeo_1() result(assertion)
3042 3 : use Constants_mod, only: IK, RK
3043 : implicit none
3044 : logical :: assertion
3045 : integer(IK) :: i
3046 : integer(IK) , parameter :: numTrial = 10_IK
3047 : integer(IK) , parameter :: SuccessStep(numTrial) = [ (i, i = 1, numTrial) ]
3048 : real(RK) , parameter :: successProb = 0.7_RK
3049 : real(RK) , parameter :: tolerance = 1.e-12_RK
3050 : real(RK) , parameter :: LogProbGeo_ref(numTrial) = [ -.3566749439387324_RK &
3051 : , -1.560647748264668_RK &
3052 : , -2.764620552590604_RK &
3053 : , -3.968593356916540_RK &
3054 : , -5.172566161242476_RK &
3055 : , -6.376538965568412_RK &
3056 : , -7.580511769894348_RK &
3057 : , -8.784484574220283_RK &
3058 : , -9.988457378546219_RK &
3059 : , -11.19243018287215_RK ]
3060 : real(RK) :: LogProbGeo(numTrial)
3061 : real(RK) :: Difference(numTrial)
3062 :
3063 3 : LogProbGeo = getLogProbGeo(numTrial, SuccessStep, successProb)
3064 33 : Difference = abs(LogProbGeo - LogProbGeo_ref) / abs(LogProbGeo_ref)
3065 33 : assertion = all( Difference < tolerance )
3066 :
3067 3 : if (Test%isDebugMode .and. .not. assertion) then
3068 : ! LCOV_EXCL_START
3069 : write(Test%outputUnit,"(*(g0,:,' '))")
3070 : write(Test%outputUnit,"(*(g0,:,' '))") "LogProbGeo_ref =", LogProbGeo_ref
3071 : write(Test%outputUnit,"(*(g0,:,' '))") "LogProbGeo =", LogProbGeo
3072 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
3073 : write(Test%outputUnit,"(*(g0,:,' '))")
3074 : end if
3075 : ! LCOV_EXCL_STOP
3076 :
3077 3 : end function test_getLogProbGeo_1
3078 :
3079 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3080 :
3081 3 : function test_getLogProbGeoCyclic_1() result(assertion)
3082 3 : use Constants_mod, only: IK, RK
3083 : implicit none
3084 : logical :: assertion
3085 : integer(IK) :: i
3086 : integer(IK) , parameter :: numTrial = 10_IK
3087 : integer(IK) , parameter :: maxNumTrial = 3_IK
3088 : integer(IK) , parameter :: SuccessStep(numTrial) = [ (i, i = 1, numTrial) ]
3089 : real(RK) , parameter :: successProb = 0.7_RK
3090 : real(RK) , parameter :: tolerance = 1.e-12_RK
3091 : real(RK) , parameter :: LogProbGeoCyclic_ref(numTrial) = [ -0.32930374714260041_RK &
3092 : , -1.53327655146853630_RK &
3093 : , -2.73724935579447240_RK &
3094 : , -3.94122216012040830_RK &
3095 : , -5.14519496444634420_RK &
3096 : , -6.34916776877228010_RK &
3097 : , -7.55314057309821600_RK &
3098 : , -8.75711337742415100_RK &
3099 : , -9.96108618175008690_RK &
3100 : , -11.1650589860760230_RK ]
3101 : real(RK) :: LogProbGeoCyclic(numTrial)
3102 : real(RK) :: Difference(numTrial)
3103 :
3104 3 : LogProbGeoCyclic = getLogProbGeoCyclic(successProb, maxNumTrial, numTrial, SuccessStep)
3105 33 : Difference = abs(LogProbGeoCyclic - LogProbGeoCyclic_ref) / abs(LogProbGeoCyclic_ref)
3106 33 : assertion = all( Difference < tolerance )
3107 :
3108 3 : if (Test%isDebugMode .and. .not. assertion) then
3109 : ! LCOV_EXCL_START
3110 : write(Test%outputUnit,"(*(g0,:,' '))")
3111 : write(Test%outputUnit,"(*(g0,:,' '))") "LogProbGeoCyclic_ref =", LogProbGeoCyclic_ref
3112 : write(Test%outputUnit,"(*(g0,:,' '))") "LogProbGeoCyclic =", LogProbGeoCyclic
3113 : write(Test%outputUnit,"(*(g0,:,' '))") "Difference =", Difference
3114 : write(Test%outputUnit,"(*(g0,:,' '))")
3115 : end if
3116 : ! LCOV_EXCL_STOP
3117 :
3118 3 : end function test_getLogProbGeoCyclic_1
3119 :
3120 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3121 :
3122 3 : function test_getSNormPDF_1() result(assertion)
3123 3 : use Constants_mod, only: IK, RK
3124 : implicit none
3125 : logical :: assertion
3126 : integer(IK) , parameter :: nd = 2_IK
3127 : real(RK) , parameter :: tolerance = 1.e-12_RK
3128 : real(RK) , parameter :: snormPDF_ref = 0.004431848411938_RK
3129 3 : real(RK) :: difference
3130 3 : real(RK) :: snormPDF
3131 3 : snormPDF = getSNormPDF(3._RK)
3132 3 : difference = abs( (snormPDF - snormPDF_ref) / snormPDF_ref )
3133 3 : assertion = difference < tolerance
3134 3 : if (Test%isDebugMode .and. .not. assertion) then
3135 : ! LCOV_EXCL_START
3136 : write(Test%outputUnit,"(*(g0,:,', '))")
3137 : write(Test%outputUnit,"(*(g0,:,', '))") "snormPDF_ref ", snormPDF_ref
3138 : write(Test%outputUnit,"(*(g0,:,', '))") "snormPDF ", snormPDF
3139 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3140 : write(Test%outputUnit,"(*(g0,:,', '))")
3141 : end if
3142 : ! LCOV_EXCL_STOP
3143 3 : end function test_getSNormPDF_1
3144 :
3145 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3146 :
3147 3 : function test_getSNormCDF_SPR_1() result(assertion)
3148 : use iso_fortran_env, only: RK => real32
3149 3 : use Constants_mod, only: IK
3150 : implicit none
3151 : logical :: assertion
3152 : integer(IK) , parameter :: nd = 2_IK
3153 : real(RK) , parameter :: tolerance = 1.e-5_RK
3154 : real(RK) , parameter :: snormCDF_ref = 0.998650101968370_RK
3155 3 : real(RK) :: difference
3156 3 : real(RK) :: snormCDF
3157 3 : snormCDF = getSNormCDF(3._RK)
3158 3 : difference = abs( (snormCDF - snormCDF_ref) / snormCDF_ref )
3159 3 : assertion = difference < tolerance
3160 3 : if (Test%isDebugMode .and. .not. assertion) then
3161 : ! LCOV_EXCL_START
3162 : write(Test%outputUnit,"(*(g0,:,', '))")
3163 : write(Test%outputUnit,"(*(g0,:,', '))") "snormCDF_ref ", snormCDF_ref
3164 : write(Test%outputUnit,"(*(g0,:,', '))") "snormCDF ", snormCDF
3165 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3166 : write(Test%outputUnit,"(*(g0,:,', '))")
3167 : end if
3168 : ! LCOV_EXCL_STOP
3169 3 : end function test_getSNormCDF_SPR_1
3170 :
3171 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3172 :
3173 3 : function test_getSNormCDF_DPR_1() result(assertion)
3174 : use iso_fortran_env, only: RK => real64
3175 3 : use Constants_mod, only: IK
3176 : implicit none
3177 : logical :: assertion
3178 : integer(IK) , parameter :: nd = 2_IK
3179 : real(RK) , parameter :: tolerance = 1.e-12_RK
3180 : real(RK) , parameter :: snormCDF_ref = 0.998650101968370_RK
3181 3 : real(RK) :: difference
3182 3 : real(RK) :: snormCDF
3183 3 : snormCDF = getSNormCDF(3._RK)
3184 3 : difference = abs( (snormCDF - snormCDF_ref) / snormCDF_ref )
3185 3 : assertion = difference < tolerance
3186 3 : if (Test%isDebugMode .and. .not. assertion) then
3187 : ! LCOV_EXCL_START
3188 : write(Test%outputUnit,"(*(g0,:,', '))")
3189 : write(Test%outputUnit,"(*(g0,:,', '))") "snormCDF_ref ", snormCDF_ref
3190 : write(Test%outputUnit,"(*(g0,:,', '))") "snormCDF ", snormCDF
3191 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3192 : write(Test%outputUnit,"(*(g0,:,', '))")
3193 : end if
3194 : ! LCOV_EXCL_STOP
3195 3 : end function test_getSNormCDF_DPR_1
3196 :
3197 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3198 :
3199 3 : function test_getNormPDF_1() result(assertion)
3200 3 : use Constants_mod, only: IK, RK
3201 : implicit none
3202 : logical :: assertion
3203 : integer(IK) , parameter :: nd = 2_IK
3204 : real(RK) , parameter :: tolerance = 1.e-12_RK
3205 : real(RK) , parameter :: avg = 2._RK
3206 : real(RK) , parameter :: std = 5._RK
3207 : real(RK) , parameter :: val = 10._RK
3208 : real(RK) , parameter :: normPDF_ref = 0.022184166935891_RK
3209 3 : real(RK) :: difference
3210 3 : real(RK) :: normPDF
3211 3 : normPDF = getNormPDF(avg,std,std**2,val)
3212 3 : difference = abs( (normPDF - normPDF_ref) / normPDF_ref )
3213 3 : assertion = difference < tolerance
3214 3 : if (Test%isDebugMode .and. .not. assertion) then
3215 : ! LCOV_EXCL_START
3216 : write(Test%outputUnit,"(*(g0,:,', '))")
3217 : write(Test%outputUnit,"(*(g0,:,', '))") "normPDF_ref ", normPDF_ref
3218 : write(Test%outputUnit,"(*(g0,:,', '))") "normPDF ", normPDF
3219 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3220 : write(Test%outputUnit,"(*(g0,:,', '))")
3221 : end if
3222 : ! LCOV_EXCL_STOP
3223 3 : end function test_getNormPDF_1
3224 :
3225 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3226 :
3227 3 : function test_getNormCDF_1() result(assertion)
3228 3 : use Constants_mod, only: IK, RK
3229 : implicit none
3230 : logical :: assertion
3231 : integer(IK) , parameter :: nd = 2_IK
3232 : real(RK) , parameter :: tolerance = 1.e-12_RK
3233 : real(RK) , parameter :: avg = 2._RK
3234 : real(RK) , parameter :: std = 5._RK
3235 : real(RK) , parameter :: val = 10._RK
3236 : real(RK) , parameter :: normCDF_ref = 0.945200708300442_RK
3237 3 : real(RK) :: difference
3238 3 : real(RK) :: normCDF
3239 3 : normCDF = getNormCDF(avg,std,val)
3240 3 : difference = abs( (normCDF - normCDF_ref) / normCDF_ref )
3241 3 : assertion = difference < tolerance
3242 3 : if (Test%isDebugMode .and. .not. assertion) then
3243 : ! LCOV_EXCL_START
3244 : write(Test%outputUnit,"(*(g0,:,', '))")
3245 : write(Test%outputUnit,"(*(g0,:,', '))") "normCDF_ref ", normCDF_ref
3246 : write(Test%outputUnit,"(*(g0,:,', '))") "normCDF ", normCDF
3247 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3248 : write(Test%outputUnit,"(*(g0,:,', '))")
3249 : end if
3250 : ! LCOV_EXCL_STOP
3251 3 : end function test_getNormCDF_1
3252 :
3253 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3254 :
3255 3 : function test_getBetaCDF_SPR_1() result(assertion)
3256 : use iso_fortran_env, only: RK => real32
3257 3 : use Constants_mod, only: IK
3258 : implicit none
3259 : logical :: assertion
3260 : integer(IK) , parameter :: nd = 2_IK
3261 : real(RK) , parameter :: tolerance = 1.e-12_RK
3262 : real(RK) , parameter :: alpha = 2._RK
3263 : real(RK) , parameter :: beta = 5._RK
3264 : real(RK) , parameter :: val = 0.5_RK
3265 : real(RK) , parameter :: betaCDF_ref = 0.890625000000000_RK
3266 3 : real(RK) :: difference
3267 3 : real(RK) :: betaCDF
3268 6 : betaCDF = getBetaCDF(alpha,beta,val)
3269 3 : difference = abs( (betaCDF - betaCDF_ref) / betaCDF_ref )
3270 3 : assertion = difference < tolerance
3271 3 : if (Test%isDebugMode .and. .not. assertion) then
3272 : ! LCOV_EXCL_START
3273 : write(Test%outputUnit,"(*(g0,:,', '))")
3274 : write(Test%outputUnit,"(*(g0,:,', '))") "betaCDF_ref ", betaCDF_ref
3275 : write(Test%outputUnit,"(*(g0,:,', '))") "betaCDF ", betaCDF
3276 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3277 : write(Test%outputUnit,"(*(g0,:,', '))")
3278 : end if
3279 : ! LCOV_EXCL_STOP
3280 3 : assertion = assertion .and. getBetaCDF(alpha,beta,-.01_RK) < 0._RK .and. getBetaCDF(alpha,beta,+1.01_RK) < 0._RK
3281 3 : end function test_getBetaCDF_SPR_1
3282 :
3283 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3284 :
3285 3 : function test_getBetaCDF_DPR_1() result(assertion)
3286 : use iso_fortran_env, only: RK => real64
3287 3 : use Constants_mod, only: IK
3288 : implicit none
3289 : logical :: assertion
3290 : integer(IK) , parameter :: nd = 2_IK
3291 : real(RK) , parameter :: tolerance = 1.e-12_RK
3292 : real(RK) , parameter :: alpha = 2._RK
3293 : real(RK) , parameter :: beta = 5._RK
3294 : real(RK) , parameter :: val = 0.5_RK
3295 : real(RK) , parameter :: betaCDF_ref = 0.890625000000000_RK
3296 3 : real(RK) :: difference
3297 3 : real(RK) :: betaCDF
3298 6 : betaCDF = getBetaCDF(alpha,beta,val)
3299 3 : difference = abs( (betaCDF - betaCDF_ref) / betaCDF_ref )
3300 3 : assertion = difference < tolerance
3301 3 : if (Test%isDebugMode .and. .not. assertion) then
3302 : ! LCOV_EXCL_START
3303 : write(Test%outputUnit,"(*(g0,:,', '))")
3304 : write(Test%outputUnit,"(*(g0,:,', '))") "betaCDF_ref ", betaCDF_ref
3305 : write(Test%outputUnit,"(*(g0,:,', '))") "betaCDF ", betaCDF
3306 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3307 : write(Test%outputUnit,"(*(g0,:,', '))")
3308 : end if
3309 : ! LCOV_EXCL_STOP
3310 3 : assertion = assertion .and. getBetaCDF(alpha,beta,-.01_RK) < 0._RK .and. getBetaCDF(alpha,beta,+1.01_RK) < 0._RK
3311 3 : end function test_getBetaCDF_DPR_1
3312 :
3313 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3314 :
3315 3 : function test_getUniformCDF_1() result(assertion)
3316 3 : use Constants_mod, only: IK, RK
3317 : implicit none
3318 : logical :: assertion
3319 : integer(IK) , parameter :: nd = 2_IK
3320 : real(RK) , parameter :: tolerance = 1.e-12_RK
3321 : real(RK) , parameter :: avg = 2._RK
3322 : real(RK) , parameter :: std = 5._RK
3323 : real(RK) , parameter :: val = 10._RK
3324 : real(RK) , parameter :: uniformCDF_ref = val
3325 3 : real(RK) :: difference
3326 3 : real(RK) :: uniformCDF
3327 3 : uniformCDF = getUniformCDF(val)
3328 3 : difference = abs( (uniformCDF - uniformCDF_ref) / uniformCDF_ref )
3329 3 : assertion = difference < tolerance
3330 3 : if (Test%isDebugMode .and. .not. assertion) then
3331 : ! LCOV_EXCL_START
3332 : write(Test%outputUnit,"(*(g0,:,', '))")
3333 : write(Test%outputUnit,"(*(g0,:,', '))") "uniformCDF_ref ", uniformCDF_ref
3334 : write(Test%outputUnit,"(*(g0,:,', '))") "uniformCDF ", uniformCDF
3335 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3336 : write(Test%outputUnit,"(*(g0,:,', '))")
3337 : end if
3338 : ! LCOV_EXCL_STOP
3339 3 : end function test_getUniformCDF_1
3340 :
3341 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3342 :
3343 3 : function test_doKS1_1() result(assertion)
3344 :
3345 3 : use Constants_mod, only: RK, IK
3346 :
3347 : implicit none
3348 :
3349 : logical :: assertion
3350 : real(RK) , parameter :: tolerance = 1.e-7_RK
3351 :
3352 : real(RK), parameter :: probKS_ref = .3763758622852317E-01_RK
3353 : real(RK), parameter :: statKS_ref = .1955719390701096_RK
3354 3 : real(RK) :: statKS
3355 3 : real(RK) :: probKS
3356 3 : real(RK) :: difference
3357 3 : type(Err_type) :: Err
3358 :
3359 : call doKS1 ( np = lenRnd &
3360 : , Point = StdNormRnd1 &
3361 : , getCDF = getSNormCDF_DPR &
3362 : , statKS = statKS &
3363 : , probKS = probKS &
3364 : , Err = Err &
3365 3 : )
3366 3 : assertion = .not. Err%occurred
3367 3 : if (.not. assertion) return
3368 :
3369 3 : difference = abs( (probKS - probKS_ref) / probKS_ref )
3370 3 : assertion = difference < tolerance
3371 :
3372 3 : if (Test%isDebugMode .and. .not. assertion) then
3373 : ! LCOV_EXCL_START
3374 : write(Test%outputUnit,"(*(g0,:,', '))")
3375 : write(Test%outputUnit,"(*(g0,:,', '))") "probKS_ref :", probKS_ref
3376 : write(Test%outputUnit,"(*(g0,:,', '))") "probKS :", probKS
3377 : write(Test%outputUnit,"(*(g0,:,', '))") "difference :", difference
3378 : write(Test%outputUnit,"(*(g0,:,', '))")
3379 : end if
3380 : ! LCOV_EXCL_STOP
3381 :
3382 3 : difference = abs( (statKS - statKS_ref) / statKS_ref )
3383 3 : assertion = assertion .and. difference < tolerance
3384 :
3385 3 : if (Test%isDebugMode .and. .not. assertion) then
3386 : ! LCOV_EXCL_START
3387 : write(Test%outputUnit,"(*(g0,:,', '))")
3388 : write(Test%outputUnit,"(*(g0,:,', '))") "statKS_ref :", statKS_ref
3389 : write(Test%outputUnit,"(*(g0,:,', '))") "statKS :", statKS
3390 : write(Test%outputUnit,"(*(g0,:,', '))") "difference :", difference
3391 : write(Test%outputUnit,"(*(g0,:,', '))")
3392 : end if
3393 : ! LCOV_EXCL_STOP
3394 :
3395 3 : end function test_doKS1_1
3396 :
3397 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3398 :
3399 3 : function test_doUniformKS1_1() result(assertion)
3400 :
3401 3 : use Constants_mod, only: RK, IK
3402 :
3403 : implicit none
3404 :
3405 : logical :: assertion
3406 : real(RK) , parameter :: tolerance = 1.e-7_RK
3407 :
3408 : real(RK), parameter :: probKS_ref = .1982797523608350_RK
3409 : real(RK), parameter :: statKS_ref = .1491359075319200_RK
3410 3 : real(RK) :: statKS
3411 3 : real(RK) :: probKS
3412 3 : real(RK) :: difference
3413 3 : type(Err_type) :: Err
3414 :
3415 : call doUniformKS1 ( np = lenRnd &
3416 : , Point = UnifRnd &
3417 : , statKS = statKS &
3418 : , probKS = probKS &
3419 : , Err = Err &
3420 3 : )
3421 3 : assertion = .not. Err%occurred
3422 3 : if (.not. assertion) return
3423 :
3424 3 : difference = abs( (probKS - probKS_ref) / probKS_ref )
3425 3 : assertion = difference < tolerance
3426 :
3427 3 : if (Test%isDebugMode .and. .not. assertion) then
3428 : ! LCOV_EXCL_START
3429 : write(Test%outputUnit,"(*(g0,:,', '))")
3430 : write(Test%outputUnit,"(*(g0,:,', '))") "probKS_ref :", probKS_ref
3431 : write(Test%outputUnit,"(*(g0,:,', '))") "probKS :", probKS
3432 : write(Test%outputUnit,"(*(g0,:,', '))") "difference :", difference
3433 : write(Test%outputUnit,"(*(g0,:,', '))")
3434 : end if
3435 : ! LCOV_EXCL_STOP
3436 :
3437 3 : difference = abs( (statKS - statKS_ref) / statKS_ref )
3438 3 : assertion = assertion .and. difference < tolerance
3439 :
3440 3 : if (Test%isDebugMode .and. .not. assertion) then
3441 : ! LCOV_EXCL_START
3442 : write(Test%outputUnit,"(*(g0,:,', '))")
3443 : write(Test%outputUnit,"(*(g0,:,', '))") "statKS_ref :", statKS_ref
3444 : write(Test%outputUnit,"(*(g0,:,', '))") "statKS :", statKS
3445 : write(Test%outputUnit,"(*(g0,:,', '))") "difference :", difference
3446 : write(Test%outputUnit,"(*(g0,:,', '))")
3447 : end if
3448 : ! LCOV_EXCL_STOP
3449 :
3450 3 : end function test_doUniformKS1_1
3451 :
3452 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3453 :
3454 3 : function test_doSortedKS2_1() result(assertion)
3455 :
3456 3 : use Constants_mod, only: RK, IK
3457 : use Sort_mod, only: sortAscending
3458 :
3459 : implicit none
3460 :
3461 : logical :: assertion
3462 : integer(IK) , parameter :: lenRnd = 50_IK
3463 : real(RK) , parameter :: tolerance = 1.e-7_RK
3464 :
3465 : real(RK), parameter :: probKS_ref = 0.056045859714425_RK
3466 : real(RK), parameter :: statKS_ref = 0.260000000000000_RK
3467 3 : real(RK) :: difference, statKS, probKS
3468 3 : type(Err_type) :: Err
3469 :
3470 3 : call sortAscending( np = lenRnd, Point = StdNormRnd1, Err = Err )
3471 3 : assertion = .not. Err%occurred
3472 3 : if (.not. assertion) return
3473 :
3474 3 : call sortAscending( np = lenRnd, Point = StdNormRnd2, Err = Err )
3475 3 : assertion = .not. Err%occurred
3476 3 : if (.not. assertion) return
3477 :
3478 : call doSortedKS2( np1 = lenRnd &
3479 : , np2 = lenRnd &
3480 : , SortedPoint1 = StdNormRnd1 &
3481 : , SortedPoint2 = StdNormRnd2 &
3482 : , statKS = statKS &
3483 : , probKS = probKS &
3484 3 : )
3485 :
3486 3 : difference = 2 * abs(probKS - probKS_ref) / (probKS_ref + probKS)
3487 3 : assertion = assertion .and. difference < tolerance
3488 :
3489 3 : if (Test%isDebugMode .and. .not. assertion) then
3490 : ! LCOV_EXCL_START
3491 : write(Test%outputUnit,"(*(g0,:,', '))")
3492 : write(Test%outputUnit,"(*(g0,:,', '))") "probKS_ref :", probKS_ref
3493 : write(Test%outputUnit,"(*(g0,:,', '))") "probKS :", probKS
3494 : write(Test%outputUnit,"(*(g0,:,', '))") "difference :", difference
3495 : write(Test%outputUnit,"(*(g0,:,', '))")
3496 : end if
3497 : ! LCOV_EXCL_STOP
3498 :
3499 3 : difference = 2 * abs(statKS - statKS_ref) / (statKS_ref + statKS)
3500 3 : assertion = difference < tolerance
3501 :
3502 3 : if (Test%isDebugMode .and. .not. assertion) then
3503 : ! LCOV_EXCL_START
3504 : write(Test%outputUnit,"(*(g0,:,', '))")
3505 : write(Test%outputUnit,"(*(g0,:,', '))") "statKS_ref :", statKS_ref
3506 : write(Test%outputUnit,"(*(g0,:,', '))") "statKS :", statKS
3507 : write(Test%outputUnit,"(*(g0,:,', '))") "difference :", difference
3508 : write(Test%outputUnit,"(*(g0,:,', '))")
3509 : end if
3510 : ! LCOV_EXCL_STOP
3511 :
3512 3 : end function test_doSortedKS2_1
3513 :
3514 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3515 :
3516 3 : function test_getHist1D_1() result(assertion)
3517 :
3518 3 : use Constants_mod, only: RK, IK
3519 :
3520 : implicit none
3521 :
3522 : logical :: assertion
3523 : integer(IK) , parameter :: nxbin = 10_IK
3524 : real(RK) , parameter :: tolerance = 1.e-7_RK
3525 : real(RK) , parameter :: Density_ref(nxbin) = [ 1.000000000000000_RK &
3526 : , 1.000000000000000_RK &
3527 : , 4.000000000000000_RK &
3528 : , 6.000000000000000_RK &
3529 : , 11.00000000000000_RK &
3530 : , 15.00000000000000_RK &
3531 : , 8.000000000000000_RK &
3532 : , 1.000000000000000_RK &
3533 : , 2.000000000000000_RK &
3534 : , 1.000000000000000_RK ]
3535 : real(RK) , parameter :: Xbin_ref(nxbin) = [ -3.068150106908867_RK &
3536 : , -2.315881996736801_RK &
3537 : , -1.563613886564735_RK &
3538 : , -.8113457763926690_RK &
3539 : , -.5907766622060301e-01_RK &
3540 : , +.6931904439514631_RK &
3541 : , +1.445458554123529_RK &
3542 : , +2.197726664295595_RK &
3543 : , +2.949994774467661_RK &
3544 : , +3.702262884639727_RK ]
3545 : real(RK) :: Density(nxbin)
3546 : real(RK) :: Xbin(nxbin)
3547 : real(RK) :: Difference(nxbin)
3548 :
3549 : call getHist1D ( method = "count" &
3550 : , xmin = minval(StdNormRnd1) - 0.5_RK &
3551 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
3552 : , nxbin = nxbin &
3553 : , np = lenRnd &
3554 : , X = StdNormRnd1 &
3555 : , Xbin = Xbin &
3556 : , Density = Density &
3557 : , errorOccurred = assertion &
3558 303 : )
3559 3 : assertion = .not. assertion
3560 3 : if (.not. assertion) return
3561 :
3562 33 : Difference = abs( (Xbin - Xbin_ref) / Xbin_ref )
3563 33 : assertion = all(Difference < tolerance)
3564 :
3565 3 : if (Test%isDebugMode .and. .not. assertion) then
3566 : ! LCOV_EXCL_START
3567 : write(Test%outputUnit,"(*(g0,:,', '))")
3568 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_ref ", Xbin_ref
3569 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin ", Xbin
3570 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", Difference
3571 : write(Test%outputUnit,"(*(g0,:,', '))")
3572 : end if
3573 : ! LCOV_EXCL_STOP
3574 :
3575 33 : Difference = abs( (Density - Density_ref) / Density_ref )
3576 33 : assertion = assertion .and. all(Difference < tolerance)
3577 :
3578 3 : if (Test%isDebugMode .and. .not. assertion) then
3579 : ! LCOV_EXCL_START
3580 : write(Test%outputUnit,"(*(g0,:,', '))")
3581 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_ref", Density_ref
3582 : write(Test%outputUnit,"(*(g0,:,', '))") "Density ", Density
3583 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3584 : write(Test%outputUnit,"(*(g0,:,', '))")
3585 : end if
3586 : ! LCOV_EXCL_STOP
3587 :
3588 3 : end function test_getHist1D_1
3589 :
3590 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3591 :
3592 3 : function test_getHist1D_2() result(assertion)
3593 :
3594 3 : use Constants_mod, only: RK, IK
3595 :
3596 : implicit none
3597 :
3598 : logical :: assertion
3599 : integer(IK) , parameter :: nxbin = 10_IK
3600 : real(RK) , parameter :: tolerance = 1.e-7_RK
3601 : real(RK) , parameter :: Density_ref(nxbin) = [ .2000000000000000E-01_RK &
3602 : , .2000000000000000E-01_RK &
3603 : , .8000000000000000E-01_RK &
3604 : , .1200000000000000_RK &
3605 : , .2200000000000000_RK &
3606 : , .3000000000000000_RK &
3607 : , .1600000000000000_RK &
3608 : , .2000000000000000E-01_RK &
3609 : , .4000000000000000E-01_RK &
3610 : , .2000000000000000E-01_RK ]
3611 : real(RK) , parameter :: Xbin_ref(nxbin) = [ -3.068150106908867_RK &
3612 : , -2.315881996736801_RK &
3613 : , -1.563613886564735_RK &
3614 : , -.8113457763926690_RK &
3615 : , -.5907766622060301e-01_RK &
3616 : , +.6931904439514631_RK &
3617 : , +1.445458554123529_RK &
3618 : , +2.197726664295595_RK &
3619 : , +2.949994774467661_RK &
3620 : , +3.702262884639727_RK ]
3621 : real(RK) :: Density(nxbin)
3622 : real(RK) :: Xbin(nxbin)
3623 : real(RK) :: Difference(nxbin)
3624 :
3625 : call getHist1D ( method = "pdf" &
3626 : , xmin = minval(StdNormRnd1) - 0.5_RK &
3627 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
3628 : , nxbin = nxbin &
3629 : , np = lenRnd &
3630 : , X = StdNormRnd1 &
3631 : , Xbin = Xbin &
3632 : , Density = Density &
3633 : , errorOccurred = assertion &
3634 303 : )
3635 3 : assertion = .not. assertion
3636 3 : if (.not. assertion) return
3637 :
3638 33 : Difference = abs( (Xbin - Xbin_ref) / Xbin_ref )
3639 33 : assertion = all(Difference < tolerance)
3640 :
3641 3 : if (Test%isDebugMode .and. .not. assertion) then
3642 : ! LCOV_EXCL_START
3643 : write(Test%outputUnit,"(*(g0,:,', '))")
3644 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_ref ", Xbin_ref
3645 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin ", Xbin
3646 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", Difference
3647 : write(Test%outputUnit,"(*(g0,:,', '))")
3648 : end if
3649 : ! LCOV_EXCL_STOP
3650 :
3651 33 : Difference = abs( (Density - Density_ref) / Density_ref )
3652 33 : assertion = assertion .and. all(Difference < tolerance)
3653 :
3654 3 : if (Test%isDebugMode .and. .not. assertion) then
3655 : ! LCOV_EXCL_START
3656 : write(Test%outputUnit,"(*(g0,:,', '))")
3657 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_ref", Density_ref
3658 : write(Test%outputUnit,"(*(g0,:,', '))") "Density ", Density
3659 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", difference
3660 : write(Test%outputUnit,"(*(g0,:,', '))")
3661 : end if
3662 : ! LCOV_EXCL_STOP
3663 :
3664 3 : end function test_getHist1D_2
3665 :
3666 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3667 :
3668 3 : function test_getHist1D_3() result(assertion)
3669 :
3670 3 : use Constants_mod, only: RK, IK
3671 :
3672 : implicit none
3673 :
3674 : logical :: assertion
3675 : integer(IK) , parameter :: nxbin = 10_IK
3676 : real(RK) :: Density(nxbin)
3677 : real(RK) :: Xbin(nxbin)
3678 :
3679 : call getHist1D ( method = "nonsense" &
3680 : , xmin = minval(StdNormRnd1) - 0.5_RK &
3681 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
3682 : , nxbin = nxbin &
3683 : , np = lenRnd &
3684 : , X = StdNormRnd1 &
3685 : , Xbin = Xbin &
3686 : , Density = Density &
3687 : , errorOccurred = assertion &
3688 303 : )
3689 :
3690 3 : end function test_getHist1D_3
3691 :
3692 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3693 :
3694 3 : function test_getHist2D_1() result(assertion)
3695 :
3696 3 : use Constants_mod, only: RK, IK
3697 :
3698 : implicit none
3699 :
3700 : logical :: assertion
3701 : integer(IK) , parameter :: nxbin = 8_IK
3702 : integer(IK) , parameter :: nybin = 7_IK
3703 : real(RK) , parameter :: tolerance = 1.e-7_RK
3704 : real(RK) , parameter :: Density_ref(nybin,nxbin) = reshape( [ .000000000000000_RK &
3705 : , 1.000000000000000_RK &
3706 : , .000000000000000_RK &
3707 : , .000000000000000_RK &
3708 : , .000000000000000_RK &
3709 : , .000000000000000_RK &
3710 : , .000000000000000_RK &
3711 : , .000000000000000_RK &
3712 : , 1.000000000000000_RK &
3713 : , .000000000000000_RK &
3714 : , 1.000000000000000_RK &
3715 : , .000000000000000_RK &
3716 : , .000000000000000_RK &
3717 : , .000000000000000_RK &
3718 : , .000000000000000_RK &
3719 : , 2.000000000000000_RK &
3720 : , 2.000000000000000_RK &
3721 : , 3.000000000000000_RK &
3722 : , 1.000000000000000_RK &
3723 : , 1.000000000000000_RK &
3724 : , .000000000000000_RK &
3725 : , .000000000000000_RK &
3726 : , 2.000000000000000_RK &
3727 : , 3.000000000000000_RK &
3728 : , 2.000000000000000_RK &
3729 : , 2.000000000000000_RK &
3730 : , 1.000000000000000_RK &
3731 : , 1.000000000000000_RK &
3732 : , .000000000000000_RK &
3733 : , 3.000000000000000_RK &
3734 : , 5.000000000000000_RK &
3735 : , 5.000000000000000_RK &
3736 : , 3.000000000000000_RK &
3737 : , 1.000000000000000_RK &
3738 : , .000000000000000_RK &
3739 : , .000000000000000_RK &
3740 : , 1.000000000000000_RK &
3741 : , 2.000000000000000_RK &
3742 : , 2.000000000000000_RK &
3743 : , 1.000000000000000_RK &
3744 : , .000000000000000_RK &
3745 : , 1.000000000000000_RK &
3746 : , 1.000000000000000_RK &
3747 : , .000000000000000_RK &
3748 : , 1.000000000000000_RK &
3749 : , .000000000000000_RK &
3750 : , .000000000000000_RK &
3751 : , .000000000000000_RK &
3752 : , .000000000000000_RK &
3753 : , 1.000000000000000_RK &
3754 : , .000000000000000_RK &
3755 : , .000000000000000_RK &
3756 : , .000000000000000_RK &
3757 : , .000000000000000_RK &
3758 : , .000000000000000_RK &
3759 : , .000000000000000_RK &
3760 : ], shape = shape(Density_ref) )
3761 : real(RK) , parameter :: Xbin_ref(nxbin) = [ -2.974116593137359_RK &
3762 : , -2.033781455422276_RK &
3763 : , -1.093446317707194_RK &
3764 : , -.1531111799921113_RK &
3765 : , .7872239577229714_RK &
3766 : , 1.727559095438054_RK &
3767 : , 2.667894233153136_RK &
3768 : , 3.608229370868219_RK ]
3769 : real(RK) , parameter :: Ybin_ref(nybin) = [ -2.038843334246188_RK &
3770 : , -1.250484167036584_RK &
3771 : , -.4621249998269794_RK &
3772 : , .3262341673826251_RK &
3773 : , 1.114593334592229_RK &
3774 : , 1.902952501801833_RK &
3775 : , 2.691311669011438_RK ]
3776 : real(RK) :: Density(nybin,nxbin)
3777 : real(RK) :: Xbin(nxbin),Xbin_diff(nxbin)
3778 : real(RK) :: Ybin(nybin),Ybin_diff(nybin)
3779 : real(RK) :: Density_diff(nybin,nxbin)
3780 :
3781 : call getHist2D ( histType = "count" &
3782 : , xmin = minval(StdNormRnd1) - 0.5_RK &
3783 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
3784 : , ymin = minval(StdNormRnd2) - 0.5_RK &
3785 : , ymax = maxval(StdNormRnd2) + 0.5_RK &
3786 : , nxbin = nxbin &
3787 : , nybin = nybin &
3788 : , np = lenRnd &
3789 : , X = StdNormRnd1 &
3790 : , Y = StdNormRnd2 &
3791 : , Xbin = Xbin &
3792 : , Ybin = Ybin &
3793 : , Density = Density &
3794 : , errorOccurred = assertion &
3795 603 : )
3796 3 : assertion = .not. assertion
3797 3 : if (.not. assertion) return
3798 :
3799 27 : Xbin_diff = abs( (Xbin - Xbin_ref) / Xbin_ref )
3800 27 : assertion = all(Xbin_diff < tolerance)
3801 :
3802 3 : if (Test%isDebugMode .and. .not. assertion) then
3803 : ! LCOV_EXCL_START
3804 : write(Test%outputUnit,"(*(g0,:,', '))")
3805 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_ref ", Xbin_ref
3806 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin ", Xbin
3807 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_diff ", Xbin_diff
3808 : write(Test%outputUnit,"(*(g0,:,', '))")
3809 : end if
3810 : ! LCOV_EXCL_STOP
3811 :
3812 24 : Ybin_diff = abs( (Ybin - Ybin_ref) / Ybin_ref )
3813 24 : assertion = all(Ybin_diff < tolerance)
3814 :
3815 3 : if (Test%isDebugMode .and. .not. assertion) then
3816 : ! LCOV_EXCL_START
3817 : write(Test%outputUnit,"(*(g0,:,', '))")
3818 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_ref ", Ybin_ref
3819 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin ", Ybin
3820 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_diff ", Ybin_diff
3821 : write(Test%outputUnit,"(*(g0,:,', '))")
3822 : end if
3823 : ! LCOV_EXCL_STOP
3824 :
3825 195 : Density_diff = abs(Density - Density_ref)
3826 195 : assertion = assertion .and. all(Density_diff < tolerance)
3827 :
3828 3 : if (Test%isDebugMode .and. .not. assertion) then
3829 : ! LCOV_EXCL_START
3830 : write(Test%outputUnit,"(*(g0,:,', '))")
3831 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_ref ", Density_ref
3832 : write(Test%outputUnit,"(*(g0,:,', '))") "Density ", Density
3833 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_diff ", Density_diff
3834 : write(Test%outputUnit,"(*(g0,:,', '))")
3835 : end if
3836 : ! LCOV_EXCL_STOP
3837 :
3838 3 : end function test_getHist2D_1
3839 :
3840 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3841 :
3842 3 : function test_getHist2D_2() result(assertion)
3843 :
3844 3 : use Constants_mod, only: RK, IK
3845 :
3846 : implicit none
3847 :
3848 : logical :: assertion
3849 : integer(IK) , parameter :: nxbin = 8_IK
3850 : integer(IK) , parameter :: nybin = 7_IK
3851 : real(RK) , parameter :: tolerance = 1.e-7_RK
3852 : real(RK) , parameter :: Density_ref(nybin,nxbin) = reshape( [ .000000000000000_RK &
3853 : , .2000000000000000E-01_RK &
3854 : , .000000000000000_RK &
3855 : , .000000000000000_RK &
3856 : , .000000000000000_RK &
3857 : , .000000000000000_RK &
3858 : , .000000000000000_RK &
3859 : , .000000000000000_RK &
3860 : , .2000000000000000E-01_RK &
3861 : , .000000000000000_RK &
3862 : , .2000000000000000E-01_RK &
3863 : , .000000000000000_RK &
3864 : , .000000000000000_RK &
3865 : , .000000000000000_RK &
3866 : , .000000000000000_RK &
3867 : , .4000000000000000E-01_RK &
3868 : , .4000000000000000E-01_RK &
3869 : , .6000000000000000E-01_RK &
3870 : , .2000000000000000E-01_RK &
3871 : , .2000000000000000E-01_RK &
3872 : , .000000000000000_RK &
3873 : , .000000000000000_RK &
3874 : , .4000000000000000E-01_RK &
3875 : , .6000000000000000E-01_RK &
3876 : , .4000000000000000E-01_RK &
3877 : , .4000000000000000E-01_RK &
3878 : , .2000000000000000E-01_RK &
3879 : , .2000000000000000E-01_RK &
3880 : , .000000000000000_RK &
3881 : , .6000000000000000E-01_RK &
3882 : , .1000000000000000_RK &
3883 : , .1000000000000000_RK &
3884 : , .6000000000000000E-01_RK &
3885 : , .2000000000000000E-01_RK &
3886 : , .000000000000000_RK &
3887 : , .000000000000000_RK &
3888 : , .2000000000000000E-01_RK &
3889 : , .4000000000000000E-01_RK &
3890 : , .4000000000000000E-01_RK &
3891 : , .2000000000000000E-01_RK &
3892 : , .000000000000000_RK &
3893 : , .2000000000000000E-01_RK &
3894 : , .2000000000000000E-01_RK &
3895 : , .000000000000000_RK &
3896 : , .2000000000000000E-01_RK &
3897 : , .000000000000000_RK &
3898 : , .000000000000000_RK &
3899 : , .000000000000000_RK &
3900 : , .000000000000000_RK &
3901 : , .2000000000000000E-01_RK &
3902 : , .000000000000000_RK &
3903 : , .000000000000000_RK &
3904 : , .000000000000000_RK &
3905 : , .000000000000000_RK &
3906 : , .000000000000000_RK &
3907 : , .000000000000000_RK &
3908 : ], shape = shape(Density_ref) )
3909 : real(RK) , parameter :: Xbin_ref(nxbin) = [ -2.974116593137359_RK &
3910 : , -2.033781455422276_RK &
3911 : , -1.093446317707194_RK &
3912 : , -.1531111799921113_RK &
3913 : , .7872239577229714_RK &
3914 : , 1.727559095438054_RK &
3915 : , 2.667894233153136_RK &
3916 : , 3.608229370868219_RK ]
3917 : real(RK) , parameter :: Ybin_ref(nybin) = [ -2.038843334246188_RK &
3918 : , -1.250484167036584_RK &
3919 : , -.4621249998269794_RK &
3920 : , .3262341673826251_RK &
3921 : , 1.114593334592229_RK &
3922 : , 1.902952501801833_RK &
3923 : , 2.691311669011438_RK ]
3924 : real(RK) :: Density(nybin,nxbin)
3925 : real(RK) :: Xbin(nxbin),Xbin_diff(nxbin)
3926 : real(RK) :: Ybin(nybin),Ybin_diff(nybin)
3927 : real(RK) :: Density_diff(nybin,nxbin)
3928 :
3929 : call getHist2D ( histType = "pdf" &
3930 : , xmin = minval(StdNormRnd1) - 0.5_RK &
3931 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
3932 : , ymin = minval(StdNormRnd2) - 0.5_RK &
3933 : , ymax = maxval(StdNormRnd2) + 0.5_RK &
3934 : , nxbin = nxbin &
3935 : , nybin = nybin &
3936 : , np = lenRnd &
3937 : , X = StdNormRnd1 &
3938 : , Y = StdNormRnd2 &
3939 : , Xbin = Xbin &
3940 : , Ybin = Ybin &
3941 : , Density = Density &
3942 : , errorOccurred = assertion &
3943 603 : )
3944 3 : assertion = .not. assertion
3945 3 : if (.not. assertion) return
3946 :
3947 27 : Xbin_diff = abs( (Xbin - Xbin_ref) / Xbin_ref )
3948 27 : assertion = all(Xbin_diff < tolerance)
3949 :
3950 3 : if (Test%isDebugMode .and. .not. assertion) then
3951 : ! LCOV_EXCL_START
3952 : write(Test%outputUnit,"(*(g0,:,', '))")
3953 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_ref ", Xbin_ref
3954 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin ", Xbin
3955 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_diff ", Xbin_diff
3956 : write(Test%outputUnit,"(*(g0,:,', '))")
3957 : end if
3958 : ! LCOV_EXCL_STOP
3959 :
3960 24 : Ybin_diff = abs( (Ybin - Ybin_ref) / Ybin_ref )
3961 24 : assertion = all(Ybin_diff < tolerance)
3962 :
3963 3 : if (Test%isDebugMode .and. .not. assertion) then
3964 : ! LCOV_EXCL_START
3965 : write(Test%outputUnit,"(*(g0,:,', '))")
3966 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_ref ", Ybin_ref
3967 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin ", Ybin
3968 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_diff ", Ybin_diff
3969 : write(Test%outputUnit,"(*(g0,:,', '))")
3970 : end if
3971 : ! LCOV_EXCL_STOP
3972 :
3973 195 : Density_diff = abs(Density - Density_ref)
3974 195 : assertion = assertion .and. all(Density_diff < tolerance)
3975 :
3976 3 : if (Test%isDebugMode .and. .not. assertion) then
3977 : ! LCOV_EXCL_START
3978 : write(Test%outputUnit,"(*(g0,:,', '))")
3979 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_ref ", Density_ref
3980 : write(Test%outputUnit,"(*(g0,:,', '))") "Density ", Density
3981 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_diff ", Density_diff
3982 : write(Test%outputUnit,"(*(g0,:,', '))")
3983 : end if
3984 : ! LCOV_EXCL_STOP
3985 :
3986 3 : end function test_getHist2D_2
3987 :
3988 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3989 :
3990 3 : function test_getHist2D_3() result(assertion)
3991 :
3992 3 : use Constants_mod, only: RK, IK
3993 :
3994 : implicit none
3995 :
3996 : logical :: assertion
3997 : integer(IK) , parameter :: nxbin = 8_IK
3998 : integer(IK) , parameter :: nybin = 7_IK
3999 : real(RK) , parameter :: tolerance = 1.e-7_RK
4000 : real(RK) , parameter :: Density_ref(nybin,nxbin) = reshape( [ .000000000000000_RK &
4001 : , 1.000000000000000_RK &
4002 : , .000000000000000_RK &
4003 : , .000000000000000_RK &
4004 : , .000000000000000_RK &
4005 : , .000000000000000_RK &
4006 : , .000000000000000_RK &
4007 : , .000000000000000_RK &
4008 : , .5000000000000000_RK &
4009 : , .000000000000000_RK &
4010 : , .5000000000000000_RK &
4011 : , .000000000000000_RK &
4012 : , .000000000000000_RK &
4013 : , .000000000000000_RK &
4014 : , .000000000000000_RK &
4015 : , .2222222222222222_RK &
4016 : , .2222222222222222_RK &
4017 : , .3333333333333333_RK &
4018 : , .1111111111111111_RK &
4019 : , .1111111111111111_RK &
4020 : , .000000000000000_RK &
4021 : , .000000000000000_RK &
4022 : , .1818181818181818_RK &
4023 : , .2727272727272727_RK &
4024 : , .1818181818181818_RK &
4025 : , .1818181818181818_RK &
4026 : , .9090909090909091E-01_RK &
4027 : , .9090909090909091E-01_RK &
4028 : , .000000000000000_RK &
4029 : , .1764705882352941_RK &
4030 : , .2941176470588235_RK &
4031 : , .2941176470588235_RK &
4032 : , .1764705882352941_RK &
4033 : , .5882352941176471E-01_RK &
4034 : , .000000000000000_RK &
4035 : , .000000000000000_RK &
4036 : , .1428571428571428_RK &
4037 : , .2857142857142857_RK &
4038 : , .2857142857142857_RK &
4039 : , .1428571428571428_RK &
4040 : , .000000000000000_RK &
4041 : , .1428571428571428_RK &
4042 : , .5000000000000000_RK &
4043 : , .000000000000000_RK &
4044 : , .5000000000000000_RK &
4045 : , .000000000000000_RK &
4046 : , .000000000000000_RK &
4047 : , .000000000000000_RK &
4048 : , .000000000000000_RK &
4049 : , 1.000000000000000_RK &
4050 : , .000000000000000_RK &
4051 : , .000000000000000_RK &
4052 : , .000000000000000_RK &
4053 : , .000000000000000_RK &
4054 : , .000000000000000_RK &
4055 : , .000000000000000_RK &
4056 : ], shape = shape(Density_ref) )
4057 : real(RK) , parameter :: Xbin_ref(nxbin) = [ -2.974116593137359_RK &
4058 : , -2.033781455422276_RK &
4059 : , -1.093446317707194_RK &
4060 : , -.1531111799921113_RK &
4061 : , .7872239577229714_RK &
4062 : , 1.727559095438054_RK &
4063 : , 2.667894233153136_RK &
4064 : , 3.608229370868219_RK ]
4065 : real(RK) , parameter :: Ybin_ref(nybin) = [ -2.038843334246188_RK &
4066 : , -1.250484167036584_RK &
4067 : , -.4621249998269794_RK &
4068 : , .3262341673826251_RK &
4069 : , 1.114593334592229_RK &
4070 : , 1.902952501801833_RK &
4071 : , 2.691311669011438_RK ]
4072 : real(RK) :: Density(nybin,nxbin)
4073 : real(RK) :: Xbin(nxbin),Xbin_diff(nxbin)
4074 : real(RK) :: Ybin(nybin),Ybin_diff(nybin)
4075 : real(RK) :: Density_diff(nybin,nxbin)
4076 :
4077 : call getHist2D ( histType = "pdf(y|x)" &
4078 : , xmin = minval(StdNormRnd1) - 0.5_RK &
4079 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
4080 : , ymin = minval(StdNormRnd2) - 0.5_RK &
4081 : , ymax = maxval(StdNormRnd2) + 0.5_RK &
4082 : , nxbin = nxbin &
4083 : , nybin = nybin &
4084 : , np = lenRnd &
4085 : , X = StdNormRnd1 &
4086 : , Y = StdNormRnd2 &
4087 : , Xbin = Xbin &
4088 : , Ybin = Ybin &
4089 : , Density = Density &
4090 : , errorOccurred = assertion &
4091 603 : )
4092 3 : assertion = .not. assertion
4093 3 : if (.not. assertion) return
4094 :
4095 27 : Xbin_diff = abs( (Xbin - Xbin_ref) / Xbin_ref )
4096 27 : assertion = all(Xbin_diff < tolerance)
4097 :
4098 3 : if (Test%isDebugMode .and. .not. assertion) then
4099 : ! LCOV_EXCL_START
4100 : write(Test%outputUnit,"(*(g0,:,', '))")
4101 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_ref ", Xbin_ref
4102 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin ", Xbin
4103 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_diff ", Xbin_diff
4104 : write(Test%outputUnit,"(*(g0,:,', '))")
4105 : end if
4106 : ! LCOV_EXCL_STOP
4107 :
4108 24 : Ybin_diff = abs( (Ybin - Ybin_ref) / Ybin_ref )
4109 24 : assertion = all(Ybin_diff < tolerance)
4110 :
4111 3 : if (Test%isDebugMode .and. .not. assertion) then
4112 : ! LCOV_EXCL_START
4113 : write(Test%outputUnit,"(*(g0,:,', '))")
4114 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_ref ", Ybin_ref
4115 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin ", Ybin
4116 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_diff ", Ybin_diff
4117 : write(Test%outputUnit,"(*(g0,:,', '))")
4118 : end if
4119 : ! LCOV_EXCL_STOP
4120 :
4121 195 : Density_diff = abs(Density - Density_ref)
4122 195 : assertion = assertion .and. all(Density_diff < tolerance)
4123 :
4124 3 : if (Test%isDebugMode .and. .not. assertion) then
4125 : ! LCOV_EXCL_START
4126 : write(Test%outputUnit,"(*(g0,:,', '))")
4127 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_ref ", Density_ref
4128 : write(Test%outputUnit,"(*(g0,:,', '))") "Density ", Density
4129 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_diff ", Density_diff
4130 : write(Test%outputUnit,"(*(g0,:,', '))")
4131 : end if
4132 : ! LCOV_EXCL_STOP
4133 :
4134 3 : end function test_getHist2D_3
4135 :
4136 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4137 :
4138 3 : function test_getHist2D_4() result(assertion)
4139 :
4140 3 : use Constants_mod, only: RK, IK
4141 :
4142 : implicit none
4143 :
4144 : logical :: assertion
4145 : integer(IK) , parameter :: nxbin = 8_IK
4146 : integer(IK) , parameter :: nybin = 7_IK
4147 : real(RK) , parameter :: tolerance = 1.e-7_RK
4148 : real(RK) , parameter :: Density_ref(nybin,nxbin) = reshape( [ .000000000000000_RK &
4149 : , .1000000000000000_RK &
4150 : , .000000000000000_RK &
4151 : , .000000000000000_RK &
4152 : , .000000000000000_RK &
4153 : , .000000000000000_RK &
4154 : , .000000000000000_RK &
4155 : , .000000000000000_RK &
4156 : , .1000000000000000_RK &
4157 : , .000000000000000_RK &
4158 : , .7692307692307693E-01_RK &
4159 : , .000000000000000_RK &
4160 : , .000000000000000_RK &
4161 : , .000000000000000_RK &
4162 : , .000000000000000_RK &
4163 : , .2000000000000000_RK &
4164 : , .1538461538461539_RK &
4165 : , .2307692307692308_RK &
4166 : , .1428571428571428_RK &
4167 : , .3333333333333333_RK &
4168 : , .000000000000000_RK &
4169 : , .000000000000000_RK &
4170 : , .2000000000000000_RK &
4171 : , .2307692307692308_RK &
4172 : , .1538461538461539_RK &
4173 : , .2857142857142857_RK &
4174 : , .3333333333333333_RK &
4175 : , .5000000000000000_RK &
4176 : , .000000000000000_RK &
4177 : , .3000000000000000_RK &
4178 : , .3846153846153846_RK &
4179 : , .3846153846153846_RK &
4180 : , .4285714285714285_RK &
4181 : , .3333333333333333_RK &
4182 : , .000000000000000_RK &
4183 : , .000000000000000_RK &
4184 : , .1000000000000000_RK &
4185 : , .1538461538461539_RK &
4186 : , .1538461538461539_RK &
4187 : , .1428571428571428_RK &
4188 : , .000000000000000_RK &
4189 : , .5000000000000000_RK &
4190 : , .5000000000000000_RK &
4191 : , .000000000000000_RK &
4192 : , .7692307692307693E-01_RK &
4193 : , .000000000000000_RK &
4194 : , .000000000000000_RK &
4195 : , .000000000000000_RK &
4196 : , .000000000000000_RK &
4197 : , .5000000000000000_RK &
4198 : , .000000000000000_RK &
4199 : , .000000000000000_RK &
4200 : , .000000000000000_RK &
4201 : , .000000000000000_RK &
4202 : , .000000000000000_RK &
4203 : , .000000000000000_RK &
4204 : ], shape = shape(Density_ref) )
4205 : real(RK) , parameter :: Xbin_ref(nxbin) = [ -2.974116593137359_RK &
4206 : , -2.033781455422276_RK &
4207 : , -1.093446317707194_RK &
4208 : , -.1531111799921113_RK &
4209 : , .7872239577229714_RK &
4210 : , 1.727559095438054_RK &
4211 : , 2.667894233153136_RK &
4212 : , 3.608229370868219_RK ]
4213 : real(RK) , parameter :: Ybin_ref(nybin) = [ -2.038843334246188_RK &
4214 : , -1.250484167036584_RK &
4215 : , -.4621249998269794_RK &
4216 : , .3262341673826251_RK &
4217 : , 1.114593334592229_RK &
4218 : , 1.902952501801833_RK &
4219 : , 2.691311669011438_RK ]
4220 : real(RK) :: Density(nybin,nxbin)
4221 : real(RK) :: Xbin(nxbin),Xbin_diff(nxbin)
4222 : real(RK) :: Ybin(nybin),Ybin_diff(nybin)
4223 : real(RK) :: Density_diff(nybin,nxbin)
4224 :
4225 : call getHist2D ( histType = "pdf(x|y)" &
4226 : , xmin = minval(StdNormRnd1) - 0.5_RK &
4227 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
4228 : , ymin = minval(StdNormRnd2) - 0.5_RK &
4229 : , ymax = maxval(StdNormRnd2) + 0.5_RK &
4230 : , nxbin = nxbin &
4231 : , nybin = nybin &
4232 : , np = lenRnd &
4233 : , X = StdNormRnd1 &
4234 : , Y = StdNormRnd2 &
4235 : , Xbin = Xbin &
4236 : , Ybin = Ybin &
4237 : , Density = Density &
4238 : , errorOccurred = assertion &
4239 603 : )
4240 3 : assertion = .not. assertion
4241 3 : if (.not. assertion) return
4242 :
4243 27 : Xbin_diff = abs( (Xbin - Xbin_ref) / Xbin_ref )
4244 27 : assertion = all(Xbin_diff < tolerance)
4245 :
4246 3 : if (Test%isDebugMode .and. .not. assertion) then
4247 : ! LCOV_EXCL_START
4248 : write(Test%outputUnit,"(*(g0,:,', '))")
4249 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_ref ", Xbin_ref
4250 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin ", Xbin
4251 : write(Test%outputUnit,"(*(g0,:,', '))") "Xbin_diff ", Xbin_diff
4252 : write(Test%outputUnit,"(*(g0,:,', '))")
4253 : end if
4254 : ! LCOV_EXCL_STOP
4255 :
4256 24 : Ybin_diff = abs( (Ybin - Ybin_ref) / Ybin_ref )
4257 24 : assertion = all(Ybin_diff < tolerance)
4258 :
4259 3 : if (Test%isDebugMode .and. .not. assertion) then
4260 : ! LCOV_EXCL_START
4261 : write(Test%outputUnit,"(*(g0,:,', '))")
4262 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_ref ", Ybin_ref
4263 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin ", Ybin
4264 : write(Test%outputUnit,"(*(g0,:,', '))") "Ybin_diff ", Ybin_diff
4265 : write(Test%outputUnit,"(*(g0,:,', '))")
4266 : end if
4267 : ! LCOV_EXCL_STOP
4268 :
4269 195 : Density_diff = abs(Density - Density_ref)
4270 195 : assertion = assertion .and. all(Density_diff < tolerance)
4271 :
4272 3 : if (Test%isDebugMode .and. .not. assertion) then
4273 : ! LCOV_EXCL_START
4274 : write(Test%outputUnit,"(*(g0,:,', '))")
4275 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_ref ", Density_ref
4276 : write(Test%outputUnit,"(*(g0,:,', '))") "Density ", Density
4277 : write(Test%outputUnit,"(*(g0,:,', '))") "Density_diff ", Density_diff
4278 : write(Test%outputUnit,"(*(g0,:,', '))")
4279 : end if
4280 : ! LCOV_EXCL_STOP
4281 :
4282 3 : end function test_getHist2D_4
4283 :
4284 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4285 :
4286 3 : function test_getHist2D_5() result(assertion)
4287 :
4288 3 : use Constants_mod, only: RK, IK
4289 :
4290 : implicit none
4291 :
4292 : logical :: assertion
4293 : integer(IK) , parameter :: nxbin = 8_IK
4294 : integer(IK) , parameter :: nybin = 7_IK
4295 : real(RK) :: Density(nybin,nxbin)
4296 : real(RK) :: Xbin(nxbin)
4297 : real(RK) :: Ybin(nybin)
4298 :
4299 : call getHist2D ( histType = "nonsense" &
4300 : , xmin = minval(StdNormRnd1) - 0.5_RK &
4301 : , xmax = maxval(StdNormRnd1) + 0.5_RK &
4302 : , ymin = minval(StdNormRnd2) - 0.5_RK &
4303 : , ymax = maxval(StdNormRnd2) + 0.5_RK &
4304 : , nxbin = nxbin &
4305 : , nybin = nybin &
4306 : , np = lenRnd &
4307 : , X = StdNormRnd1 &
4308 : , Y = StdNormRnd2 &
4309 : , Xbin = Xbin &
4310 : , Ybin = Ybin &
4311 : , Density = Density &
4312 : , errorOccurred = assertion &
4313 603 : )
4314 :
4315 3 : end function test_getHist2D_5
4316 :
4317 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4318 :
4319 3 : function test_getQuantile_1() result(assertion)
4320 :
4321 3 : use Constants_mod, only: RK, IK
4322 :
4323 : implicit none
4324 :
4325 : logical :: assertion
4326 : integer(IK) , parameter :: nq = 9_IK
4327 : integer(IK) , parameter :: Weight(lenRnd) = [ 2_IK, 2_IK, 1_IK, 0_IK, 0_IK, 1_IK, 2_IK, 2_IK, 1_IK, 0_IK, 0_IK, 0_IK, 1_IK, 2_IK &
4328 : , 2_IK, 1_IK, 0_IK, 0_IK, 1_IK, 2_IK, 2_IK, 1_IK, 0_IK, 0_IK, 1_IK, 2_IK, 2_IK, 1_IK &
4329 : , 0_IK, 0_IK, 1_IK, 2_IK, 2_IK, 2_IK, 1_IK, 0_IK, 0_IK, 1_IK, 2_IK, 2_IK, 1_IK, 0_IK &
4330 : , 0_IK, 1_IK, 2_IK, 2_IK, 1_IK, 0_IK, 0_IK, 1_IK ]
4331 : real(RK) , parameter :: SortedQuantileProbability(nq) = [0._RK, 0.05_RK, 0.1_RK, 0.25_RK, 0.5_RK, 0.75_RK, 0.9_RK, 0.95_RK, 1._RK]
4332 : real(RK) , parameter :: Quantile_ref(nq) = [ -2.944284161994900_RK &
4333 : , -2.258846861003650_RK &
4334 : , -1.711516418853700_RK &
4335 : , -.3034409247860160_RK &
4336 : , .3251905394561980_RK &
4337 : , 1.034693009917860_RK &
4338 : , 1.489697607785470_RK &
4339 : , 1.630235289164730_RK &
4340 : , 3.578396939725760_RK ]
4341 : real(RK) :: Quantile(nq)
4342 : real(RK) :: Difference(nq)
4343 :
4344 : Quantile = getQuantile ( np = lenRnd &
4345 : , nq = nq &
4346 : , SortedQuantileProbability = SortedQuantileProbability &
4347 : , Point = StdNormRnd1 &
4348 : , Weight = Weight &
4349 : , sumWeight = sum(Weight) &
4350 3 : )
4351 :
4352 30 : Difference = (Quantile - Quantile_ref)
4353 30 : assertion = all(Difference==0._RK)
4354 :
4355 3 : if (Test%isDebugMode .and. .not. assertion) then
4356 : ! LCOV_EXCL_START
4357 : write(Test%outputUnit,"(*(g0,:,', '))")
4358 : write(Test%outputUnit,"(*(g0,:,', '))") "Quantile_ref ", Quantile_ref
4359 : write(Test%outputUnit,"(*(g0,:,', '))") "Quantile ", Quantile
4360 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", Difference
4361 : write(Test%outputUnit,"(*(g0,:,', '))")
4362 : end if
4363 : ! LCOV_EXCL_STOP
4364 :
4365 3 : end function test_getQuantile_1
4366 :
4367 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4368 :
4369 3 : function test_getQuantile_2() result(assertion)
4370 :
4371 3 : use Constants_mod, only: RK, IK
4372 :
4373 : implicit none
4374 :
4375 : logical :: assertion
4376 : integer(IK) , parameter :: nq = 9_IK
4377 : real(RK) , parameter :: SortedQuantileProbability(nq) = [0._RK, 0.05_RK, 0.1_RK, 0.25_RK, 0.5_RK, 0.75_RK, 0.9_RK, 0.95_RK, 1._RK]
4378 : real(RK) , parameter :: Quantile_ref(nq) = [ -2.944284161994900_RK &
4379 : , -1.711516418853700_RK &
4380 : , -1.307688296305270_RK &
4381 : , -.4335920223056840_RK &
4382 : , .3192067391655020_RK &
4383 : , 1.034693009917860_RK &
4384 : , 1.489697607785470_RK &
4385 : , 2.769437029884880_RK &
4386 : , 3.578396939725760_RK ]
4387 : real(RK) :: Quantile(nq)
4388 : real(RK) :: Difference(nq)
4389 :
4390 : Quantile = getQuantile ( np = lenRnd &
4391 : , nq = nq &
4392 : , SortedQuantileProbability = SortedQuantileProbability &
4393 : , Point = StdNormRnd1 &
4394 3 : )
4395 :
4396 30 : Difference = (Quantile - Quantile_ref)
4397 30 : assertion = all(Difference==0._RK)
4398 :
4399 3 : if (Test%isDebugMode .and. .not. assertion) then
4400 : ! LCOV_EXCL_START
4401 : write(Test%outputUnit,"(*(g0,:,', '))")
4402 : write(Test%outputUnit,"(*(g0,:,', '))") "Quantile_ref ", Quantile_ref
4403 : write(Test%outputUnit,"(*(g0,:,', '))") "Quantile ", Quantile
4404 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", Difference
4405 : write(Test%outputUnit,"(*(g0,:,', '))")
4406 : end if
4407 : ! LCOV_EXCL_STOP
4408 :
4409 3 : end function test_getQuantile_2
4410 :
4411 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4412 :
4413 : !> \brief
4414 : !> Test whether [flatten_2D()](@ref statistics_mod::flatten_2d) can successfully flatten an input weighted 2D array.
4415 3 : function test_flatten_1() result(assertion)
4416 :
4417 3 : use Constants_mod, only: RK, IK
4418 :
4419 : implicit none
4420 :
4421 : integer(IK) :: i
4422 : logical :: assertion
4423 : integer(IK) , parameter :: nd = 2_IK
4424 : integer(IK) , parameter :: np = 4_IK
4425 : integer(IK) , parameter :: Weight(np) = [-1_IK, 2_IK, 0_IK, 1_IK]
4426 : integer(IK) , parameter :: sumWeight = 3_IK
4427 : real(RK) , parameter :: tolerance = 1.e-14_RK
4428 : real(RK) , parameter :: Point(nd,np) = reshape( [( real(i,RK),i=1,nd*np )], shape = shape(Point) )
4429 : real(RK) , parameter :: FlattenedPoint_ref(nd,3) = reshape([ 3._RK, 4._RK, 3._RK, 4._RK, 7._RK, 8._RK, 7._RK, 8._RK ], shape = shape(FlattenedPoint_ref) )
4430 : real(RK), allocatable :: FlattenedPoint(:,:)
4431 : real(RK), allocatable :: Difference(:,:)
4432 :
4433 30 : FlattenedPoint = flatten_2D(nd, np, Point, Weight)
4434 :
4435 30 : Difference = abs(FlattenedPoint - FlattenedPoint_ref)
4436 9 : assertion = all(shape(FlattenedPoint)==[nd,sumWeight])
4437 30 : assertion = assertion .and. all(Difference < tolerance)
4438 :
4439 3 : if (Test%isDebugMode .and. .not. assertion) then
4440 : ! LCOV_EXCL_START
4441 : write(Test%outputUnit,"(*(g0,:,', '))")
4442 : write(Test%outputUnit,"(*(g0,:,', '))") "FlattenedPoint_ref ", FlattenedPoint_ref
4443 : write(Test%outputUnit,"(*(g0,:,', '))") "FlattenedPoint ", FlattenedPoint
4444 : write(Test%outputUnit,"(*(g0,:,', '))") "difference ", Difference
4445 : write(Test%outputUnit,"(*(g0,:,', '))")
4446 : end if
4447 : ! LCOV_EXCL_STOP
4448 :
4449 3 : end function test_flatten_1
4450 :
4451 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4452 :
4453 : end module Test_Statistics_mod
|
