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
44 : !> This file implements the body of the `Setup_smod` submodules of the `ParaDRAM_mod` and `ParaDISE_mod` modules.
45 : !>
46 : !> \remark
47 : !> This module requires preprocessing, prior to compilation.
48 : !>
49 : !> \author Amir Shahmoradi
50 :
51 : #if defined PARADRAM
52 :
53 : #define ParaDXXX ParaDRAM
54 : #define ParaDXXX_type ParaDRAM_type
55 : #define ParaDXXXProposalAbstract_mod ParaDRAMProposalAbstract_mod
56 :
57 : #elif defined PARADISE
58 :
59 : #define ParaDXXX ParaDISE
60 : #define ParaDXXX_type ParaDISE_type
61 : #define ParaDXXXProposalAbstract_mod ParaDISEProposalAbstract_mod
62 :
63 : #else
64 : #error "Unrecognized sampler in ParaDXXX_mod@Setup_mod.inc.f90"
65 : #endif
66 :
67 : !submodule (ParaDRAM_mod) Setup_smod
68 :
69 : !use Constants_mod, only: IK, RK ! gfortran 9.3 compile crashes with this line
70 : implicit none
71 :
72 : character(*), parameter :: SUBMODULE_NAME = MODULE_NAME // "@Setup_smod"
73 :
74 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
76 :
77 : contains
78 :
79 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81 :
82 : !> \brief
83 : !> This procedure is a method of [ParaDRAM_type](@ref paradram_type) and [ParaDISE_type](@ref paradise_type) classes.
84 : !> Setup the sampler and run the simulation.
85 : !>
86 : !> \remark
87 : !> There are only two mandatory pieces of information that must be provided when calling this subroutine.
88 : !> The rest of the parameters can be also set from within an external input file, whose path can be given
89 : !> as an input argument `inputFile` to this subroutine.
90 : !>
91 : !> @param[in] ndim : The number of dimensions of the domain of the objective function to be sampled.
92 : !> @param[in] getLogFunc : The first point.
93 : !> @param[in] inputFile : The path to the external input file (**optional**).
94 : !>
95 : !> \remark
96 : !> For detailed descriptions of the rest of the input parameters of this subroutine,
97 : !> see: https://www.cdslab.org/paramonte/notes/usage/paradram/specifications/
98 : !>
99 : !> \remark
100 : !> This procedure requires preprocessing.
101 1047 : module subroutine runSampler( self &
102 : , ndim &
103 : , getLogFunc &
104 : , inputFile &
105 : ! ParaMonte variables
106 : , sampleSize &
107 : , randomSeed &
108 : , description &
109 : , outputFileName &
110 : , outputDelimiter &
111 : , chainFileFormat &
112 0 : , variableNameList &
113 : , restartFileFormat &
114 : , outputColumnWidth &
115 : , overwriteRequested &
116 : , outputRealPrecision &
117 : , silentModeRequested &
118 78 : , domainLowerLimitVec &
119 78 : , domainUpperLimitVec &
120 : , parallelizationModel &
121 : , progressReportPeriod &
122 : , targetAcceptanceRate &
123 : , mpiFinalizeRequested &
124 : , maxNumDomainCheckToWarn &
125 : , maxNumDomainCheckToStop &
126 : ! ParaMCMC variables
127 : , chainSize &
128 12 : , startPointVec &
129 : , sampleRefinementCount &
130 : , sampleRefinementMethod &
131 : , randomStartPointRequested &
132 24 : , randomStartPointDomainLowerLimitVec &
133 24 : , randomStartPointDomainUpperLimitVec &
134 : ! ParaDRAM variables
135 : , scaleFactor &
136 : , proposalModel &
137 18 : , proposalStartStdVec &
138 18 : , proposalStartCorMat &
139 12 : , proposalStartCovMat &
140 : , adaptiveUpdateCount &
141 : , adaptiveUpdatePeriod &
142 : , greedyAdaptationCount &
143 : , delayedRejectionCount &
144 : , burninAdaptationMeasure &
145 1047 : , delayedRejectionScaleFactorVec &
146 : ) ! result(self)
147 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
148 : !!DEC$ ATTRIBUTES DLLEXPORT :: ParaDXXX_type
149 : !DEC$ ATTRIBUTES DLLEXPORT :: runSampler
150 : #endif
151 : use, intrinsic :: iso_fortran_env, only: output_unit, stat_stopped_image
152 : #if defined PARADRAM
153 : use ParaDRAMProposalUniform_mod, only: ProposalUniform_type => Proposal_type
154 : use ParaDRAMProposalNormal_mod, only: ProposalNormal_type => Proposal_type
155 : #elif defined PARADISE
156 : use ParaDISEProposalUniform_mod, only: ProposalUniform_type => Proposal_type
157 : use ParaDISEProposalNormal_mod, only: ProposalNormal_type => Proposal_type
158 : #endif
159 : use ParaMonteLogFunc_mod, only: getLogFunc_proc
160 : use Decoration_mod, only: GENERIC_OUTPUT_FORMAT
161 : use Decoration_mod, only: GENERIC_TABBED_FORMAT
162 : use Decoration_mod, only: getGenericFormat
163 : use Decoration_mod, only: INDENT
164 : use Statistics_mod, only: getCorMatUpperFromCovMatUpper
165 : use Statistics_mod, only: getWeiSamCovUppMeanTrans
166 : use Statistics_mod, only: getQuantile
167 : use Statistics_mod, only: getMean
168 : use ParaMonte_mod, only: QPROB
169 : use Constants_mod, only: IK, RK, NLC, PMSM, UNDEFINED, POSINF_RK
170 : use DateTime_mod, only: getNiceDateTime
171 : use String_mod, only: num2str
172 : use Matrix_mod, only: getEye
173 :
174 : implicit none
175 :
176 : ! self
177 : class(ParaDXXX_type), intent(inout) :: self
178 :
179 : ! mandatory variables
180 : integer(IK) , intent(in) :: ndim
181 : procedure(getLogFunc_proc) :: getLogFunc
182 :
183 : character(*), intent(in), optional :: inputFile
184 :
185 : ! ParaMonte variables
186 : integer(IK) , intent(in), optional :: sampleSize
187 : integer(IK) , intent(in), optional :: randomSeed
188 : character(*), intent(in), optional :: description
189 : character(*), intent(in), optional :: outputFileName
190 : character(*), intent(in), optional :: outputDelimiter
191 : character(*), intent(in), optional :: chainFileFormat
192 : character(*), intent(in), optional :: variableNameList(ndim)
193 : character(*), intent(in), optional :: restartFileFormat
194 : integer(IK) , intent(in), optional :: outputColumnWidth
195 : logical , intent(in), optional :: overwriteRequested
196 : integer(IK) , intent(in), optional :: outputRealPrecision
197 : real(RK) , intent(in), optional :: domainLowerLimitVec(ndim)
198 : real(RK) , intent(in), optional :: domainUpperLimitVec(ndim)
199 : logical , intent(in), optional :: silentModeRequested
200 : character(*), intent(in), optional :: parallelizationModel
201 : integer(IK) , intent(in), optional :: progressReportPeriod
202 : real(RK) , intent(in), optional :: targetAcceptanceRate(2)
203 : logical , intent(in), optional :: mpiFinalizeRequested
204 : integer(IK) , intent(in), optional :: maxNumDomainCheckToWarn
205 : integer(IK) , intent(in), optional :: maxNumDomainCheckToStop
206 :
207 : ! ParaMCMC variables
208 : integer(IK) , intent(in), optional :: chainSize
209 : real(RK) , intent(in), optional :: startPointVec(ndim)
210 : integer(IK) , intent(in), optional :: sampleRefinementCount
211 : character(*), intent(in), optional :: sampleRefinementMethod
212 : logical , intent(in), optional :: randomStartPointRequested
213 : real(RK) , intent(in), optional :: randomStartPointDomainLowerLimitVec(ndim)
214 : real(RK) , intent(in), optional :: randomStartPointDomainUpperLimitVec(ndim)
215 :
216 : ! ParaDRAM variables
217 : character(*), intent(in), optional :: scaleFactor
218 : character(*), intent(in), optional :: proposalModel
219 : real(RK) , intent(in), optional :: proposalStartStdVec(ndim)
220 : real(RK) , intent(in), optional :: proposalStartCorMat(ndim,ndim)
221 : real(RK) , intent(in), optional :: proposalStartCovMat(ndim,ndim)
222 : integer(IK) , intent(in), optional :: adaptiveUpdateCount
223 : integer(IK) , intent(in), optional :: adaptiveUpdatePeriod
224 : integer(IK) , intent(in), optional :: greedyAdaptationCount
225 : integer(IK) , intent(in), optional :: delayedRejectionCount
226 : real(RK) , intent(in), optional :: burninAdaptationMeasure
227 : real(RK) , intent(in), optional :: delayedRejectionScaleFactorVec(:)
228 :
229 : character(*), parameter :: PROCEDURE_NAME = SUBMODULE_NAME // "@runSampler()"
230 :
231 : !type(ProposalSymmetric_type), target:: ProposalSymmetric
232 : integer(IK) :: i, iq, effectiveSampleSize
233 1047 : character(:), allocatable :: msg, formatStr, formatStrInt, formatStrReal, formatAllReal
234 : real(RK) , allocatable :: ContiguousChain(:,:) ! used to avoid temporary array creation and the compiler warning message in debug mode
235 1047 : real(RK) :: mcmcSamplingEfficiency
236 :
237 :
238 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
239 : ! Initialize SpecBase variables, then check existence of inputFile and open it and return the unit file, if it exists
240 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
241 :
242 : call self%setupParaMonte( nd = ndim &
243 : , name = PMSM%ParaDXXX &
244 : , inputFile = inputFile &
245 1047 : )
246 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
247 1047 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
248 : #endif
249 1047 : if (self%Err%occurred) then
250 : ! LCOV_EXCL_START
251 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
252 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
253 : return
254 : ! LCOV_EXCL_STOP
255 : end if
256 :
257 :
258 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
259 : ! Initialize ParaMCMC variables
260 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
261 :
262 1047 : call self%setupParaMCMC()
263 :
264 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
265 : ! Initialize ParaDXXX variables
266 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267 :
268 1047 : self%SpecDRAM = SpecDRAM_type( self%nd%val, self%name ) ! , self%SpecMCMC%ChainSize%def )
269 :
270 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
271 : ! read variables from input file if it exists
272 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
273 :
274 1047 : call self%getSpecFromInputFile(ndim)
275 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
276 1047 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
277 : #endif
278 1047 : if (self%Err%occurred) then
279 : ! LCOV_EXCL_START
280 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
281 : call self%abort( Err = self%Err, prefix = self%brand, newline = "\n", outputUnit = self%LogFile%unit )
282 : return
283 : ! LCOV_EXCL_STOP
284 : end if
285 :
286 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
287 : ! read variables from argument list if needed
288 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
289 :
290 1047 : call self%setWarnAboutProcArgHasPriority()
291 1047 : if (self%procArgNeeded) then
292 : call self%SpecBase%setFromInputArgs ( Err = self%Err & ! LCOV_EXCL_LINE
293 : , sampleSize = sampleSize & ! LCOV_EXCL_LINE
294 : , randomSeed = randomSeed & ! LCOV_EXCL_LINE
295 : , description = description & ! LCOV_EXCL_LINE
296 : , outputFileName = outputFileName & ! LCOV_EXCL_LINE
297 : , outputDelimiter = outputDelimiter & ! LCOV_EXCL_LINE
298 : , chainFileFormat = chainFileFormat & ! LCOV_EXCL_LINE
299 : , variableNameList = variableNameList & ! LCOV_EXCL_LINE
300 : , restartFileFormat = restartFileFormat & ! LCOV_EXCL_LINE
301 : , outputColumnWidth = outputColumnWidth & ! LCOV_EXCL_LINE
302 : , overwriteRequested = overwriteRequested & ! LCOV_EXCL_LINE
303 : , outputRealPrecision = outputRealPrecision & ! LCOV_EXCL_LINE
304 : , silentModeRequested = silentModeRequested & ! LCOV_EXCL_LINE
305 : , domainLowerLimitVec = domainLowerLimitVec & ! LCOV_EXCL_LINE
306 : , domainUpperLimitVec = domainUpperLimitVec & ! LCOV_EXCL_LINE
307 : , parallelizationModel = parallelizationModel & ! LCOV_EXCL_LINE
308 : , progressReportPeriod = progressReportPeriod & ! LCOV_EXCL_LINE
309 : , targetAcceptanceRate = targetAcceptanceRate & ! LCOV_EXCL_LINE
310 : , mpiFinalizeRequested = mpiFinalizeRequested & ! LCOV_EXCL_LINE
311 : , maxNumDomainCheckToWarn = maxNumDomainCheckToWarn & ! LCOV_EXCL_LINE
312 : , maxNumDomainCheckToStop = maxNumDomainCheckToStop & ! LCOV_EXCL_LINE
313 1041 : )
314 : call self%SpecMCMC%setFromInputArgs ( chainSize = chainSize & ! LCOV_EXCL_LINE
315 : , scaleFactor = scaleFactor & ! LCOV_EXCL_LINE
316 : , startPointVec = startPointVec & ! LCOV_EXCL_LINE
317 : , proposalModel = proposalModel & ! LCOV_EXCL_LINE
318 : , proposalStartStdVec = proposalStartStdVec & ! LCOV_EXCL_LINE
319 : , proposalStartCorMat = proposalStartCorMat & ! LCOV_EXCL_LINE
320 : , proposalStartCovMat = proposalStartCovMat & ! LCOV_EXCL_LINE
321 : , sampleRefinementCount = sampleRefinementCount & ! LCOV_EXCL_LINE
322 : , sampleRefinementMethod = sampleRefinementMethod & ! LCOV_EXCL_LINE
323 : , randomStartPointRequested = randomStartPointRequested & ! LCOV_EXCL_LINE
324 : , randomStartPointDomainLowerLimitVec = randomStartPointDomainLowerLimitVec & ! LCOV_EXCL_LINE
325 : , randomStartPointDomainUpperLimitVec = randomStartPointDomainUpperLimitVec & ! LCOV_EXCL_LINE
326 1041 : )
327 : call self%SpecDRAM%setFromInputArgs ( adaptiveUpdateCount = adaptiveUpdateCount & ! LCOV_EXCL_LINE
328 : , adaptiveUpdatePeriod = adaptiveUpdatePeriod & ! LCOV_EXCL_LINE
329 : , greedyAdaptationCount = greedyAdaptationCount & ! LCOV_EXCL_LINE
330 : , delayedRejectionCount = delayedRejectionCount & ! LCOV_EXCL_LINE
331 : , burninAdaptationMeasure = burninAdaptationMeasure & ! LCOV_EXCL_LINE
332 : , delayedRejectionScaleFactorVec = delayedRejectionScaleFactorVec & ! LCOV_EXCL_LINE
333 1041 : )
334 : end if
335 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
336 1047 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
337 : #endif
338 1047 : if (self%Err%occurred) then
339 : ! LCOV_EXCL_START
340 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
341 : call self%abort( Err = self%Err, prefix = self%brand, newline = "\n", outputUnit = self%LogFile%unit )
342 : return
343 : ! LCOV_EXCL_STOP
344 : end if
345 :
346 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
347 : ! Now depending on the requested parallelism type, determine the process leader/rooter type and open output files
348 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
349 :
350 1047 : self%Image%isLeader = self%Image%isFirst .or. self%SpecBase%ParallelizationModel%isMultiChain
351 1047 : self%Image%isRooter = .not. self%Image%isLeader
352 :
353 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354 : ! setup log file and open output files
355 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
356 :
357 1047 : call self%setupOutputFiles()
358 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
359 1047 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
360 : #endif
361 1047 : if (self%Err%occurred) then
362 : ! LCOV_EXCL_START
363 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
364 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
365 : return
366 : ! LCOV_EXCL_STOP
367 : end if
368 :
369 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
370 : ! report variable values to the log file
371 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
372 :
373 1035 : if (self%isFreshRun) then
374 969 : call self%SpecBase%reportValues(prefix=self%brand,outputUnit=self%LogFile%unit,isLeaderImage=self%Image%isLeader)
375 969 : call self%SpecMCMC%reportValues(prefix=self%brand,outputUnit=self%LogFile%unit,isLeaderImage=self%Image%isLeader,methodName=self%name,splashModeRequested=self%SpecBase%SilentModeRequested%isFalse)
376 969 : call self%SpecDRAM%reportValues(prefix=self%brand,outputUnit=self%LogFile%unit,isLeaderImage=self%Image%isLeader,splashModeRequested=self%SpecBase%SilentModeRequested%isFalse)
377 : end if
378 :
379 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
380 : ! perform variable value sanity checks
381 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
382 :
383 1035 : self%Err%msg = ""
384 1035 : self%Err%occurred = .false.
385 :
386 1035 : call self%SpecBase%checkForSanity(Err = self%Err, methodName = self%name)
387 1035 : call self%SpecMCMC%checkForSanity(SpecBase = self%SpecBase, methodName = self%name, Err = self%Err, nd = ndim)
388 1035 : call self%SpecDRAM%checkForSanity(Err = self%Err, methodName = self%name, nd = ndim)
389 :
390 : !if (self%Image%isLeader) then
391 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
392 1035 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
393 : #endif
394 1035 : if (self%Err%occurred) then
395 348 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
396 348 : call self%abort( Err = self%Err, prefix = self%brand, newline = "\n", outputUnit = self%LogFile%unit )
397 348 : return
398 : end if
399 : !end if
400 :
401 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
402 : ! setup output files
403 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
404 :
405 : ! setup chain file
406 :
407 : block
408 : use ParaMonteChainFileContents_mod, only: ChainFileContents_type
409 1101 : self%Chain = ChainFileContents_type( ndim = ndim, variableNameList = self%SpecBase%VariableNameList%Val )
410 : end block
411 :
412 687 : if (self%isFreshRun) then
413 627 : if (self%Image%isLeader) then
414 : call self%Chain%writeHeader ( ndim = ndim &
415 : , chainFileUnit = self%ChainFile%unit &
416 : , isBinary = self%SpecBase%ChainFileFormat%isBinary &
417 : , chainFileFormat = self%ChainFile%format &
418 237 : )
419 237 : flush(self%ChainFile%unit)
420 : end if
421 : else
422 60 : call self%Chain%getLenHeader(ndim=ndim,isBinary=self%SpecBase%ChainFileFormat%isBinary,chainFileFormat=self%ChainFile%format)
423 : end if
424 :
425 : ! setup progress file
426 :
427 687 : if (self%Image%isLeader) then
428 269 : if (.not.self%isFreshRun) then
429 32 : read(self%TimeFile%unit, *, iostat = self%Err%stat) ! read the header line of the time file, only by leader images
430 : end if
431 269 : if (self%isFreshRun .or. is_iostat_end(self%Err%stat)) then
432 237 : write( self%TimeFile%unit, self%TimeFile%format ) "NumFuncCallTotal" &
433 237 : , "NumFuncCallAccepted" &
434 237 : , "MeanAcceptanceRateSinceStart" &
435 237 : , "MeanAcceptanceRateSinceLastReport" &
436 237 : , "TimeElapsedSinceLastReportInSeconds" &
437 237 : , "TimeElapsedSinceStartInSeconds" &
438 474 : , "TimeRemainedToFinishInSeconds"
439 237 : flush(self%TimeFile%unit)
440 : end if
441 : end if
442 :
443 : ! The format setup of setupOutputFiles() uses the generic g0 edit descriptor. Here the format is revised to be more specific.
444 : ! g0 edit descriptor format is slightly more arbitrary and compiler-dependent.
445 :
446 687 : if (self%SpecBase%OutputColumnWidth%val>0) then
447 : self%TimeFile%format = "(" // &
448 : "2I" // self%SpecBase%OutputColumnWidth%str // ",'" // self%SpecBase%OutputDelimiter%val // &
449 : "',5(E" // self%SpecBase%OutputColumnWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3" // &
450 6 : ",:,'" // self%SpecBase%OutputDelimiter%val // "'))"
451 : self%ChainFile%format = "(" // &
452 : "2(I" // self%SpecBase%OutputColumnWidth%str // &
453 : ",'" // self%SpecBase%OutputDelimiter%val // "')," // &
454 : "2(E" // self%SpecBase%OutputColumnWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3" // &
455 : ",'" // self%SpecBase%OutputDelimiter%val // "')," // &
456 : "2(I" // self%SpecBase%OutputColumnWidth%str // &
457 : ",'" // self%SpecBase%OutputDelimiter%val // "')" // &
458 : "," // num2str(3_IK+self%nd%val) // &
459 : "(E" // self%SpecBase%OutputColumnWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3" // &
460 : ",:,'" // self%SpecBase%OutputDelimiter%val // "')" // &
461 6 : ")"
462 : end if
463 :
464 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
465 : ! setup the proposal distribution
466 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
467 :
468 687 : if (allocated(self%Proposal)) deallocate(self%Proposal)
469 : if (self%SpecMCMC%ProposalModel%isNormal) then ! LCOV_EXCL_LINE
470 : allocate( self%Proposal, source = ProposalNormal_type ( ndim = ndim & ! LCOV_EXCL_LINE
471 : , SpecBase = self%SpecBase & ! LCOV_EXCL_LINE
472 : , SpecMCMC = self%SpecMCMC & ! LCOV_EXCL_LINE
473 : , SpecDRAM = self%SpecDRAM & ! LCOV_EXCL_LINE
474 : , Image = self%Image & ! LCOV_EXCL_LINE
475 : , name = self%name & ! LCOV_EXCL_LINE
476 : , brand = self%brand & ! LCOV_EXCL_LINE
477 : , LogFile = self%LogFile & ! LCOV_EXCL_LINE
478 : , RestartFile = self%RestartFile & ! LCOV_EXCL_LINE
479 : , isFreshRun = self%isFreshRun & ! LCOV_EXCL_LINE
480 567 : ) )
481 : elseif (self%SpecMCMC%ProposalModel%isUniform) then ! LCOV_EXCL_LINE
482 : allocate( self%Proposal, source = ProposalUniform_type ( ndim = ndim & ! LCOV_EXCL_LINE
483 : , SpecBase = self%SpecBase & ! LCOV_EXCL_LINE
484 : , SpecMCMC = self%SpecMCMC & ! LCOV_EXCL_LINE
485 : , SpecDRAM = self%SpecDRAM & ! LCOV_EXCL_LINE
486 : , Image = self%Image & ! LCOV_EXCL_LINE
487 : , name = self%name & ! LCOV_EXCL_LINE
488 : , brand = self%brand & ! LCOV_EXCL_LINE
489 : , LogFile = self%LogFile & ! LCOV_EXCL_LINE
490 : , RestartFile = self%RestartFile & ! LCOV_EXCL_LINE
491 : , isFreshRun = self%isFreshRun & ! LCOV_EXCL_LINE
492 120 : ) )
493 : else
494 : ! LCOV_EXCL_START
495 : self%Err%occurred = .true.
496 : self%Err%msg = PROCEDURE_NAME // ": Internal error occurred. Unsupported proposal distribution for " // self%name // "."
497 : call self%abort( Err = self%Err, prefix = self%brand, newline = "\n", outputUnit = self%LogFile%unit )
498 : error stop
499 : return
500 : ! LCOV_EXCL_STOP
501 : end if
502 :
503 : #if (defined MATLAB_ENABLED || defined PYTHON_ENABLED || defined R_ENABLED) && !defined CAF_ENABLED && !defined MPI_ENABLED
504 : block
505 : #if defined PARADRAM
506 : use ParaDXXXProposalAbstract_mod, only: ProposalErr
507 : #elif defined PARADISE
508 : use ParaDXXXProposalAbstract_mod, only: ProposalErr
509 : #endif
510 :
511 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
512 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
513 : #endif
514 : if (ProposalErr%occurred) then
515 : self%Err%occurred = .true.
516 : self%Err%msg = PROCEDURE_NAME // ProposalErr%msg
517 : call self%abort( Err = self%Err, prefix = self%brand, newline = "\n", outputUnit = self%LogFile%unit )
518 : return
519 : end if
520 : end block
521 : #endif
522 :
523 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
524 : ! run ParaDXXX kernel
525 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
526 :
527 687 : if (self%isFreshRun .and. self%Image%isLeader) then
528 : call self%Decor%writeDecoratedText ( text = "\nStarting the " // self%name // " sampling - " // getNiceDateTime() // "\n" &
529 : , marginTop = 1_IK &
530 : , marginBot = 1_IK &
531 : , newline = "\n" &
532 237 : , outputUnit = self%LogFile%unit )
533 : end if
534 :
535 687 : call self%runKernel( getLogFunc = getLogFunc )
536 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
537 687 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
538 : #endif
539 : ! relevant only for MATLAB / Python
540 687 : if (self%Err%occurred) then
541 6 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
542 6 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
543 6 : return
544 : end if
545 :
546 681 : if (self%isFreshRun .and. self%Image%isLeader) then
547 : call self%Decor%writeDecoratedText ( text = "\nExiting the " // self%name // " sampling - " // getNiceDateTime() // "\n" &
548 : , marginTop = 1_IK &
549 : , marginBot = 1_IK &
550 : , newline = "\n" &
551 263 : , outputUnit = self%LogFile%unit )
552 : end if
553 :
554 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
555 : ! start ParaDXXX post-processing
556 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
557 :
558 : ! WARNING: For any error that happens within the leader if block below, the control must be passed to after the block
559 : ! WARNING: This is essential for error handling during testing in parallel mode.
560 :
561 681 : blockLeaderPostProcessing: if (self%Image%isLeader) then
562 :
563 263 : formatStrInt = "('" // INDENT // "',1A" // self%LogFile%maxColWidth%str // ",*(I" // self%LogFile%maxColWidth%str // "))"
564 263 : formatStrReal = "('" // INDENT // "',1A" // self%LogFile%maxColWidth%str // ",*(E" // self%LogFile%maxColWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3))"
565 :
566 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
567 :
568 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
569 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.numFuncCall.accepted"
570 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
571 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%NumFunCall%accepted
572 263 : msg = "This is the total number of accepted function calls (unique samples)."
573 263 : call self%reportDesc(msg)
574 :
575 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
576 :
577 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
578 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.numFuncCall.acceptedRejected"
579 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
580 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%NumFunCall%acceptedRejected
581 263 : msg = "This is the total number of accepted or rejected function calls."
582 263 : call self%reportDesc(msg)
583 :
584 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
585 :
586 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
587 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.numFuncCall.acceptedRejectedDelayed"
588 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
589 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%NumFunCall%acceptedRejectedDelayed
590 263 : msg = "This is the total number of accepted or rejected or delayed-rejection (if any requested) function calls."
591 263 : call self%reportDesc(msg)
592 :
593 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
594 :
595 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
596 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.numFuncCall.acceptedRejectedDelayedUnused"
597 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
598 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%NumFunCall%acceptedRejectedDelayedUnused
599 263 : msg = "This is the total number of accepted or rejected or unused function calls (by all processes, including delayed rejections, if any requested)."
600 263 : call self%reportDesc(msg)
601 :
602 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
603 :
604 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
605 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.efficiency.meanAcceptanceRate"
606 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
607 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted)
608 263 : msg = "This is the average MCMC acceptance rate."
609 263 : call self%reportDesc(msg)
610 :
611 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
612 :
613 263 : mcmcSamplingEfficiency = real(self%Stats%NumFunCall%accepted,kind=RK) / real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
614 :
615 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
616 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.efficiency.acceptedOverAcceptedRejected"
617 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
618 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) mcmcSamplingEfficiency
619 : msg = "This is the MCMC sampling efficiency given the accepted and rejected function calls, that is, &
620 263 : &the number of accepted function calls divided by the number of (accepted + rejected) function calls."
621 263 : call self%reportDesc(msg)
622 :
623 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
624 :
625 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
626 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.efficiency.acceptedOverAcceptedRejectedDelayed"
627 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
628 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) real(self%Stats%NumFunCall%accepted,kind=RK) / real(self%Stats%NumFunCall%acceptedRejectedDelayed,kind=RK)
629 : msg = "This is the MCMC sampling efficiency given the accepted, rejected, and delayed-rejection (if any requested) function calls, that is, &
630 263 : &the number of accepted function calls divided by the number of (accepted + rejected + delayed-rejection) function calls."
631 263 : call self%reportDesc(msg)
632 :
633 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
634 :
635 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
636 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.efficiency.acceptedOverAcceptedRejectedDelayedUnused"
637 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
638 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) real(self%Stats%NumFunCall%accepted,kind=RK) / real(self%Stats%NumFunCall%acceptedRejectedDelayedUnused,kind=RK)
639 : msg = "This is the MCMC sampling efficiency given the accepted, rejected, delayed-rejection (if any requested), and unused function calls, that is, &
640 263 : &the number of accepted function calls divided by the number of (accepted + rejected + delayed-rejection + unused) function calls."
641 263 : call self%reportDesc(msg)
642 :
643 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
644 :
645 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
646 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.time.total"
647 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
648 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Timer%Time%total
649 263 : msg = "This is the total runtime in seconds."
650 263 : call self%reportDesc(msg)
651 :
652 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
653 :
654 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
655 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.time.perFuncCallAccepted"
656 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
657 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Timer%Time%total / real(self%Stats%NumFunCall%accepted,kind=RK)
658 263 : msg = "This is the average effective time cost of each accepted function call, in seconds."
659 263 : call self%reportDesc(msg)
660 :
661 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
662 :
663 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
664 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.time.perFuncCallAcceptedRejected"
665 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
666 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Timer%Time%total / real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
667 263 : msg = "This is the average effective time cost of each accepted or rejected function call, in seconds."
668 263 : call self%reportDesc(msg)
669 :
670 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
671 :
672 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
673 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.time.perFuncCallAcceptedRejectedDelayed"
674 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
675 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Timer%Time%total / real(self%Stats%NumFunCall%acceptedRejectedDelayed,kind=RK)
676 263 : msg = "This is the average effective time cost of each accepted or rejected function call (including delayed-rejections, if any requested), in seconds."
677 263 : call self%reportDesc(msg)
678 :
679 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
680 :
681 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
682 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.time.perFuncCallAcceptedRejectedDelayedUnused"
683 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
684 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Timer%Time%total / real(self%Stats%NumFunCall%acceptedRejectedDelayedUnused,kind=RK)
685 263 : msg = "This is the average effective time cost of each accepted or rejected or unused function call (including delayed-rejections, if any requested), in seconds."
686 263 : call self%reportDesc(msg)
687 :
688 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
689 :
690 263 : if (self%Image%count==1_IK) then
691 : msg = UNDEFINED ! LCOV_EXCL_LINE
692 : else ! LCOV_EXCL_LINE
693 : msg = num2str( self%Stats%avgCommTimePerFunCall ) ! LCOV_EXCL_LINE
694 : end if
695 :
696 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
697 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.time.perInterProcessCommunication"
698 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
699 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) msg
700 263 : msg = "This is the average time cost of inter-process communications per used (accepted or rejected or delayed-rejection) function call, in seconds."
701 263 : call self%reportDesc(msg)
702 :
703 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
704 :
705 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
706 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.time.perFuncCall"
707 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
708 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%avgTimePerFunCalInSec
709 263 : msg = "This is the average pure time cost of each function call, in seconds."
710 263 : call self%reportDesc(msg)
711 :
712 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
713 :
714 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
715 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.current.numProcess"
716 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
717 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Image%count
718 : #if defined CAF_ENABLED || defined MPI_ENABLED
719 263 : if (self%Image%count==1_IK) then
720 : msg = self%name // " " // &
721 : "is being used in parallel mode but with only one processor. This is computationally inefficient. &
722 0 : &Consider using the serial version of the code or provide more processes at runtime if it is beneficial."
723 : call self%note ( prefix = self%brand &
724 : , outputUnit = output_unit &
725 : , newline = NLC &
726 : , marginTop = 3_IK &
727 : , marginBot = 0_IK &
728 0 : , msg = msg )
729 : end if
730 : #endif
731 263 : msg = "This is the number of processes (images) used in this simulation."
732 263 : call self%reportDesc(msg)
733 :
734 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
735 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin speedup compute %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
736 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
737 :
738 : blockSpeedup: block
739 :
740 : use Parallelism_mod, only: ForkJoin_type
741 : integer(IK) :: imageCount
742 263 : character(:), allocatable :: formatIn, formatScaling
743 : logical :: isForkJoinParallelism
744 263 : type(ForkJoin_type) :: ForkJoin
745 263 : character(:), allocatable :: undefinedInfinity
746 :
747 263 : if (self%SpecBase%ParallelizationModel%isMultiChain) then
748 54 : undefinedInfinity = "+INFINITY"
749 : else
750 209 : undefinedInfinity = UNDEFINED
751 : end if
752 :
753 263 : isForkJoinParallelism = self%Image%count > 1_IK .and. self%SpecBase%ParallelizationModel%isSinglChain
754 :
755 263 : formatScaling = "('" // INDENT // "',10(E" // self%LogFile%maxColWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3))" ! ,:,','
756 :
757 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
758 : ! compute the effective MCMC efficiency from the processor contributions and the current strong scaling
759 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
760 :
761 263 : imageCount = 1_IK
762 263 : if (isForkJoinParallelism) imageCount = self%Image%count
763 : ForkJoin = ForkJoin_type( processCount = imageCount & ! LCOV_EXCL_LINE
764 : , lenProcessID = self%Stats%NumFunCall%accepted & ! LCOV_EXCL_LINE
765 : , ProcessID = self%Chain%ProcessID & ! LCOV_EXCL_LINE
766 : , successProb = mcmcSamplingEfficiency & ! LCOV_EXCL_LINE
767 : , seqSecTime = epsilon(1._RK) & ! LCOV_EXCL_LINE ! time cost of the sequential section of the code, which is negligible here
768 : , parSecTime = self%Stats%avgTimePerFunCalInSec & ! LCOV_EXCL_LINE
769 : , comSecTime = self%Stats%avgCommTimePerFunCall & ! LCOV_EXCL_LINE
770 263 : )
771 263 : if (ForkJoin%Err%occurred) then
772 : ! LCOV_EXCL_START
773 : self%Err = ForkJoin%Err
774 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
775 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
776 : exit blockLeaderPostProcessing
777 : return
778 : ! LCOV_EXCL_STOP
779 : end if
780 :
781 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
782 :
783 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
784 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.processContribution"
785 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
786 263 : block
787 : integer, parameter :: NCOL = 40
788 263 : character(:), allocatable :: formatInteger
789 263 : formatInteger = "('"//INDENT//"',"//num2str(NCOL)//"(I0,:,' '))"
790 1052 : do imageCount = 1, ForkJoin%Contribution%count, NCOL
791 526 : write(self%LogFile%unit,formatInteger) ForkJoin%Contribution%Frequency( imageCount : min( imageCount+NCOL-1, ForkJoin%Contribution%count ) )
792 : end do
793 : end block
794 : msg = "These are contributions of individual processes to the construction of the MCMC chain. &
795 : &Essentially, they represent the total number of accepted states by the corresponding processor, &
796 : &starting from the first processor to the last. This information is mostly informative in parallel &
797 263 : &Fork-Join (singleChain) simulations."
798 263 : call self%reportDesc(msg)
799 :
800 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
801 :
802 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
803 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.processContribution.geometricFit.successProbNormFac"
804 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
805 263 : if (isForkJoinParallelism) then
806 209 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) ForkJoin%SuccessProb%PowellMinimum%xmin(1), exp(ForkJoin%SuccessProb%PowellMinimum%xmin(2))
807 : else
808 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) UNDEFINED
809 : end if
810 : msg = "These are the parameters of the Geometric fit to the distribution of the processor contributions &
811 : &to the construction of the MCMC chain (the processor contributions are reported in the first column &
812 : &of the output chain file. The fit has the following form: "//NLC//NLC// &
813 : " ProcessConstribution(i) = successProbNormFac(1) * successProbNormFac(2) * (1-successProbNormFac(1))^(i-1)"//NLC// &
814 : " / (1 - (1 - successProbNormFac(1))^numProcess)"//NLC//NLC// &
815 : "where i is the ID of the processor (starting from index 1), numProcess is the total number of processes &
816 : &used in the simulation, and successProbNormFac(1) is equivalent to an effective MCMC sampling efficiency &
817 : &computed from contributions of individual processes to the MCMC chain and successProbNormFac(2) is a &
818 263 : &normalization constant."
819 263 : call self%reportDesc(msg)
820 :
821 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
822 :
823 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
824 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.current.speedup"
825 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
826 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) ForkJoin%Speedup%current
827 263 : msg = "This is the estimated maximum speedup gained via "//self%SpecBase%ParallelizationModel%val//" parallelization model compared to serial mode."
828 263 : call self%reportDesc(msg)
829 :
830 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
831 :
832 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
833 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.optimal.current.numProcess"
834 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
835 263 : if (isForkJoinParallelism) then
836 209 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) ForkJoin%Speedup%Maximum%nproc
837 : else
838 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) undefinedInfinity
839 : end if
840 263 : msg = "This is the predicted optimal number of physical computing processes for "//self%SpecBase%ParallelizationModel%val//" parallelization model, given the current MCMC sampling efficiency."
841 263 : call self%reportDesc(msg)
842 :
843 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
844 :
845 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
846 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.optimal.current.speedup"
847 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
848 263 : if (isForkJoinParallelism) then
849 209 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) ForkJoin%Speedup%Maximum%value
850 209 : msg = ""
851 209 : if (ForkJoin%Speedup%current<1._RK) then
852 209 : formatIn = "(g0.6)"
853 : msg = "The time cost of calling the user-provided objective function must be at least " // &
854 : num2str(1._RK/ForkJoin%Speedup%current,formatIn) // " times more (that is, ~" // &
855 : num2str(10**6*self%Stats%avgTimePerFunCalInSec/ForkJoin%Speedup%current,formatIn) // &
856 : " microseconds) to see any performance benefits from " // &
857 209 : self%SpecBase%ParallelizationModel%val // " parallelization model for this simulation. "
858 209 : if (ForkJoin%Speedup%Maximum%nproc==1_IK) then
859 : msg = msg// "Consider simulating in the serial mode in the future (if used within &
860 209 : &the same computing environment and with the same configuration as used here)."
861 : else
862 : msg = msg// "Consider simulating on " // num2str(ForkJoin%Speedup%Maximum%nproc) // " processors in the future &
863 0 : &(if used within the same computing environment and with the same configuration as used here)."
864 : end if
865 : call self%note ( prefix = self%brand &
866 : , outputUnit = output_unit &
867 : , newline = NLC &
868 : , marginTop = 3_IK &
869 : , marginBot = 0_IK &
870 209 : , msg = msg )
871 209 : msg = NLC // msg
872 : end if
873 : else
874 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) undefinedInfinity
875 : end if
876 263 : msg = "This is the predicted optimal maximum speedup gained via "//self%SpecBase%ParallelizationModel%val//" parallelization model, given the current MCMC sampling efficiency." // msg
877 263 : call self%reportDesc(msg)
878 :
879 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
880 :
881 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
882 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.optimal.current.scaling.strong.speedup"
883 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
884 263 : if (isForkJoinParallelism) then
885 418 : do imageCount = 1, ForkJoin%Speedup%count, 10
886 418 : write(self%LogFile%unit,formatScaling) ForkJoin%Speedup%Scaling(imageCount:min(imageCount+9_IK,ForkJoin%Speedup%count))
887 : end do
888 : else
889 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) UNDEFINED
890 : end if
891 : msg = "This is the predicted strong-scaling speedup behavior of the "//self%SpecBase%ParallelizationModel%val//" parallelization model, &
892 263 : &given the current MCMC sampling efficiency, for increasing numbers of processes, starting from a single process."
893 263 : call self%reportDesc(msg)
894 :
895 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
896 : ! compute the absolute optimal parallelism efficiency under any MCMC sampling efficiency
897 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
898 :
899 263 : if (isForkJoinParallelism) then
900 : ForkJoin = ForkJoin_type( processCount = self%Image%count & ! LCOV_EXCL_LINE
901 : , lenProcessID = self%Stats%NumFunCall%accepted & ! LCOV_EXCL_LINE
902 : , ProcessID = self%Chain%ProcessID & ! LCOV_EXCL_LINE
903 : , successProb = 0._RK & ! LCOV_EXCL_LINE
904 : , seqSecTime = epsilon(1._RK) & ! LCOV_EXCL_LINE ! time cost of the sequential section of the code, which is negligible here
905 : , parSecTime = self%Stats%avgTimePerFunCalInSec & ! LCOV_EXCL_LINE
906 : , comSecTime = self%Stats%avgCommTimePerFunCall & ! LCOV_EXCL_LINE
907 209 : )
908 209 : if (ForkJoin%Err%occurred) then
909 : ! LCOV_EXCL_START
910 : self%Err = ForkJoin%Err
911 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
912 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
913 : exit blockLeaderPostProcessing
914 : return
915 : ! LCOV_EXCL_STOP
916 : end if
917 : end if ! otherwise, use the previously-generated ForkJoin
918 :
919 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
920 :
921 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
922 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.optimal.absolute.numProcess"
923 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) ! LCOV_EXCL_LINE
924 263 : if (isForkJoinParallelism) then
925 209 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) ForkJoin%Speedup%Maximum%nproc
926 : else
927 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) undefinedInfinity
928 : end if
929 263 : msg = "This is the predicted absolute optimal number of physical computing processes for "//self%SpecBase%ParallelizationModel%val//" parallelization model, under any MCMC sampling efficiency."
930 263 : call self%reportDesc(msg)
931 :
932 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
933 :
934 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
935 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.optimal.absolute.speedup"
936 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) ! LCOV_EXCL_LINE
937 263 : if (isForkJoinParallelism) then
938 209 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) ForkJoin%Speedup%Maximum%value
939 : else
940 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) undefinedInfinity
941 : end if
942 : msg = "This is the predicted absolute optimal maximum speedup gained via "//self%SpecBase%ParallelizationModel%val//" parallelization model, under any MCMC sampling efficiency. &
943 : &This simulation will likely NOT benefit from any additional computing processors beyond the predicted absolute optimal number, " // num2str(ForkJoin%Speedup%Maximum%nproc) // ", &
944 : &in the above. This is true for any value of MCMC sampling efficiency. Keep in mind that the predicted absolute optimal number of processors is just an estimate &
945 : &whose accuracy depends on many runtime factors, including the topology of the communication network being used, the number of processors per node, &
946 263 : &and the number of tasks to each processor or node."
947 263 : call self%reportDesc(msg)
948 :
949 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
950 :
951 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
952 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.parallelism.optimal.absolute.scaling.strong.speedup"
953 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) ! LCOV_EXCL_LINE
954 263 : if (isForkJoinParallelism) then
955 418 : do imageCount = 1, ForkJoin%Speedup%count, 10
956 418 : write(self%LogFile%unit,formatScaling) ForkJoin%Speedup%Scaling(imageCount:min(imageCount+9_IK,ForkJoin%Speedup%count))
957 : end do
958 : else
959 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) UNDEFINED
960 : end if
961 : msg = "This is the predicted absolute strong-scaling speedup behavior of the "//self%SpecBase%ParallelizationModel%val//" parallelization model, &
962 263 : &under any MCMC sampling efficiency, for increasing numbers of processes, starting from a single process."
963 526 : call self%reportDesc(msg)
964 :
965 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
966 :
967 : end block blockSpeedup
968 :
969 :
970 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
971 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end speedup compute %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
972 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
973 :
974 :
975 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
976 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.compact.burnin.location.likelihoodBased"
977 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
978 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%BurninLoc%compact
979 263 : msg = "This is the burnin location in the compact chain, based on the occurrence likelihood."
980 263 : call self%reportDesc(msg)
981 :
982 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
983 :
984 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
985 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.compact.burnin.location.adaptationBased"
986 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
987 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%AdaptationBurninLoc%compact
988 263 : msg = "This is the burnin location in the compact chain, based on the value of burninAdaptationMeasure simulation specification."
989 263 : call self%reportDesc(msg)
990 :
991 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
992 :
993 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
994 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.burnin.location.likelihoodBased"
995 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
996 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%BurninLoc%verbose
997 263 : msg = "This is the burnin location in the verbose (Markov) chain, based on the occurrence likelihood."
998 263 : call self%reportDesc(msg)
999 :
1000 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1001 :
1002 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1003 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.burnin.location.adaptationBased"
1004 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1005 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%AdaptationBurninLoc%verbose
1006 263 : msg = "This is the burnin location in the verbose (Markov) chain, based on the value of burninAdaptationMeasure simulation specification."
1007 263 : call self%reportDesc(msg)
1008 :
1009 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1010 :
1011 : ! reset BurninLoc to the maximum value
1012 :
1013 263 : if (self%Stats%AdaptationBurninLoc%compact>self%Stats%BurninLoc%compact) then
1014 1 : self%Stats%BurninLoc%compact = self%Stats%AdaptationBurninLoc%compact
1015 1 : self%Stats%BurninLoc%verbose = self%Stats%AdaptationBurninLoc%verbose
1016 : end if
1017 :
1018 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1019 :
1020 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1021 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.logFunc.max"
1022 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1023 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%LogFuncMode%val
1024 263 : msg = "This is the maximum logFunc value (the maximum of the user-specified objective function)."
1025 263 : call self%reportDesc(msg)
1026 :
1027 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1028 :
1029 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1030 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.compact.logFunc.max.location"
1031 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1032 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%LogFuncMode%Loc%compact
1033 263 : msg = "This is the location of the first occurrence of the maximum logFunc in the compact chain."
1034 263 : call self%reportDesc(msg)
1035 :
1036 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1037 :
1038 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1039 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.logFunc.max.location"
1040 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1041 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%LogFuncMode%Loc%verbose
1042 263 : msg = "This is the location of the first occurrence of the maximum logFunc in the verbose (Markov) chain."
1043 263 : call self%reportDesc(msg)
1044 :
1045 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1046 :
1047 263 : formatAllReal = "('" // INDENT // "',*(E" // self%LogFile%maxColWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3))"
1048 :
1049 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1050 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.logFunc.max.state"
1051 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1052 683 : write(self%LogFile%unit,self%LogFile%format) (trim(adjustl(self%SpecBase%VariableNameList%Val(i))), i=1, ndim)
1053 263 : write(self%LogFile%unit,formatAllReal) self%Stats%LogFuncMode%Crd
1054 263 : msg = "This is the coordinates, within the domain of the user-specified objective function, where the maximum logFunc occurs."
1055 263 : call self%reportDesc(msg)
1056 :
1057 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1058 : ! Compute the statistical properties of the MCMC chain
1059 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1060 :
1061 263 : if (self%Image%isFirst) then
1062 : call self%note ( prefix = self%brand &
1063 : , outputUnit = output_unit &
1064 : , newline = NLC &
1065 : , marginTop = 3_IK &
1066 227 : , msg = "Computing the statistical properties of the Markov chain..." )
1067 : end if
1068 :
1069 : call self%Decor%writeDecoratedText ( text = "\nThe statistical properties of the Markov chain\n" &
1070 : , marginTop = 1_IK &
1071 : , marginBot = 1_IK &
1072 : , newline = "\n" &
1073 263 : , outputUnit = self%LogFile%unit )
1074 :
1075 94289 : self%Stats%Chain%count = sum(self%Chain%Weight(self%Stats%BurninLoc%compact:self%Chain%count%compact))
1076 :
1077 : ! compute the covariance and correlation upper-triangle matrices
1078 :
1079 263 : if (allocated(self%Stats%Chain%Mean)) deallocate(self%Stats%Chain%Mean); allocate(self%Stats%Chain%Mean(ndim))
1080 263 : if (allocated(self%Stats%Chain%CovMat)) deallocate(self%Stats%Chain%CovMat); allocate(self%Stats%Chain%CovMat(ndim,ndim))
1081 : call getWeiSamCovUppMeanTrans ( np = self%Chain%count%compact - self%Stats%BurninLoc%compact + 1_IK & ! LCOV_EXCL_LINE
1082 : , sumWeight = self%Stats%Chain%count & ! LCOV_EXCL_LINE
1083 : , nd = ndim & ! LCOV_EXCL_LINE
1084 : , Point = self%Chain%State(1:ndim,self%Stats%BurninLoc%compact:self%Chain%count%compact) & ! LCOV_EXCL_LINE
1085 : , Weight = self%Chain%Weight(self%Stats%BurninLoc%compact:self%Chain%count%compact) & ! LCOV_EXCL_LINE
1086 : , CovMatUpper = self%Stats%Chain%CovMat & ! LCOV_EXCL_LINE
1087 : , Mean = self%Stats%Chain%Mean & ! LCOV_EXCL_LINE
1088 263 : )
1089 :
1090 2585 : self%Stats%Chain%CorMat = getCorMatUpperFromCovMatUpper(nd=ndim,CovMatUpper=self%Stats%Chain%CovMat)
1091 :
1092 : ! transpose the covariance and correlation matrices
1093 :
1094 683 : do i = 1, ndim
1095 420 : self%Stats%Chain%CorMat(i,i) = 1._RK
1096 577 : self%Stats%Chain%CorMat(i+1:ndim,i) = self%Stats%Chain%CorMat(i,i+1:ndim)
1097 840 : self%Stats%Chain%CovMat(i+1:ndim,i) = self%Stats%Chain%CovMat(i,i+1:ndim)
1098 : end do
1099 :
1100 : ! compute the quantiles
1101 :
1102 263 : if (allocated(self%Stats%Chain%Quantile)) deallocate(self%Stats%Chain%Quantile)
1103 263 : allocate(self%Stats%Chain%Quantile(QPROB%count,ndim))
1104 158342 : ContiguousChain = transpose(self%Chain%State(1:ndim,self%Stats%BurninLoc%compact:self%Chain%count%compact)) ! avoid temporary array creation and the warning message in debug mode
1105 683 : do i = 1, ndim
1106 2520 : self%Stats%Chain%Quantile(1:QPROB%count,i) = getQuantile( np = self%Chain%count%compact - self%Stats%BurninLoc%compact + 1_IK &
1107 : , nq = QPROB%count & ! LCOV_EXCL_LINE
1108 : , SortedQuantileProbability = QPROB%Value & ! LCOV_EXCL_LINE
1109 : , Point = ContiguousChain(:,i) & ! LCOV_EXCL_LINE
1110 1680 : , Weight = self%Chain%Weight(self%Stats%BurninLoc%compact:self%Chain%count%compact) &
1111 : , sumWeight = self%Stats%Chain%count & ! LCOV_EXCL_LINE
1112 4463 : )
1113 : end do
1114 :
1115 : ! report the MCMC chain statistics
1116 :
1117 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1118 :
1119 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1120 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.length.burninExcluded"
1121 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1122 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%Chain%count
1123 263 : msg = "This is the length of the verbose (Markov) Chain excluding burnin."
1124 263 : call self%reportDesc(msg)
1125 :
1126 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1127 :
1128 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1129 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.avgStd"
1130 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1131 263 : write(self%LogFile%unit,self%LogFile%format) "variableName", "Mean", "Standard Deviation"
1132 683 : do i = 1, ndim
1133 683 : write(self%LogFile%unit,formatStrReal) trim(adjustl(self%SpecBase%VariableNameList%Val(i))), self%Stats%Chain%Mean(i), sqrt(self%Stats%Chain%CovMat(i,i))
1134 : end do
1135 263 : msg = "This is the mean and standard deviation of the verbose (Markov) chain variables."
1136 263 : call self%reportDesc(msg)
1137 :
1138 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1139 :
1140 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1141 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.covmat"
1142 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1143 683 : write(self%LogFile%unit,self%LogFile%format) "", (trim(adjustl(self%SpecBase%VariableNameList%Val(i))),i=1,ndim)
1144 683 : do i = 1, ndim
1145 683 : write(self%LogFile%unit,formatStrReal) trim(adjustl(self%SpecBase%VariableNameList%Val(i))), self%Stats%Chain%CovMat(1:ndim,i)
1146 : end do
1147 263 : msg = "This is the covariance matrix of the verbose (Markov) chain."
1148 263 : call self%reportDesc(msg)
1149 :
1150 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1151 :
1152 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1153 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.cormat"
1154 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1155 683 : write(self%LogFile%unit,self%LogFile%format) "", (trim(adjustl(self%SpecBase%VariableNameList%Val(i))),i=1,ndim)
1156 683 : do i = 1, ndim
1157 683 : write(self%LogFile%unit,formatStrReal) trim(adjustl(self%SpecBase%VariableNameList%Val(i))), self%Stats%Chain%CorMat(1:ndim,i)
1158 : end do
1159 263 : msg = "This is the correlation matrix of the verbose (Markov) chain."
1160 263 : call self%reportDesc(msg)
1161 :
1162 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1163 :
1164 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1165 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.verbose.quantile"
1166 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1167 683 : write(self%LogFile%unit,self%LogFile%format) "Quantile", (trim(adjustl(self%SpecBase%VariableNameList%Val(i))),i=1,ndim)
1168 2630 : do iq = 1, QPROB%count
1169 6410 : write(self%LogFile%unit,formatStrReal) trim(adjustl(QPROB%Name(iq))), (self%Stats%Chain%Quantile(iq,i),i=1,ndim)
1170 : end do
1171 263 : msg = "This is the quantiles table of the variables of the verbose (Markov) chain."
1172 263 : call self%reportDesc(msg)
1173 :
1174 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1175 : ! Generate the i.i.d. sample statistics and output file (if requested)
1176 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1177 :
1178 : ! report refined sample statistics, and generate output refined sample if requested.
1179 :
1180 263 : if (self%Image%isFirst) then
1181 : call self%note ( prefix = self%brand &
1182 : , outputUnit = output_unit &
1183 : , newline = NLC &
1184 227 : , msg = "Computing the final decorrelated sample size..." )
1185 : end if
1186 :
1187 : call self%RefinedChain%get ( CFC = self%Chain &
1188 : , Err = self%Err &
1189 : , burninLoc = self%Stats%BurninLoc%compact &
1190 : , sampleRefinementCount = self%SpecMCMC%SampleRefinementCount%val &
1191 : , sampleRefinementMethod = self%SpecMCMC%SampleRefinementMethod%val &
1192 263 : )
1193 :
1194 263 : if (self%Err%occurred) then
1195 : ! LCOV_EXCL_START
1196 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
1197 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
1198 : exit blockLeaderPostProcessing
1199 : return
1200 : ! LCOV_EXCL_STOP
1201 : end if
1202 :
1203 : ! compute the maximum integrated autocorrelation times for each variable
1204 :
1205 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1206 :
1207 263 : formatStr = "('" // INDENT // "',2I" // self%LogFile%maxColWidth%str // ",*(E" // self%LogFile%maxColWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3))"
1208 :
1209 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1210 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.iac"
1211 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1212 683 : write(self%LogFile%unit,self%LogFile%format) "RefinementStage","SampleSize","IAC_SampleLogFunc",("IAC_"//trim(adjustl(self%SpecBase%VariableNameList%Val(i))),i=1,ndim)
1213 980 : do i = 0, self%RefinedChain%numRefinement
1214 980 : write(self%LogFile%unit,formatStr) i, self%RefinedChain%Count(i)%Verbose, self%RefinedChain%IAC(0:ndim,i)
1215 : end do
1216 263 : msg = "This is the table of the Integrated Autocorrelation (IAC) of individual variables in the verbose (Markov) chain, at increasing stages of chain refinements."
1217 263 : if (self%RefinedChain%numRefinement==0_IK) then
1218 : msg = msg // NLC // "The user-specified sampleRefinementCount ("// num2str(self%SpecMCMC%SampleRefinementCount%val) // ") &
1219 : &is too small to ensure an accurate computation of the decorrelated i.i.d. effective sample size. &
1220 12 : &No refinement of the Markov chain was performed."
1221 : end if
1222 263 : call self%reportDesc(msg)
1223 :
1224 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1225 :
1226 : ! report the final Effective Sample Size (ESS) based on IAC
1227 :
1228 : !blockEffectiveSampleSize: associate( effectiveSampleSize => sum(self%RefinedChain%Weight(1:self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact)) )
1229 263 : effectiveSampleSize = self%RefinedChain%Count(self%RefinedChain%numRefinement)%verbose
1230 :
1231 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1232 :
1233 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1234 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.ess"
1235 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1236 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) effectiveSampleSize
1237 263 : msg = "This is the estimated Effective (decorrelated) Sample Size (ESS) of the final refined chain."
1238 263 : call self%reportDesc(msg)
1239 :
1240 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1241 :
1242 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1243 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.efficiency.essOverAccepted"
1244 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1245 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) real(effectiveSampleSize,kind=RK) / real(self%Stats%NumFunCall%accepted,kind=RK)
1246 : msg = "This is the effective MCMC sampling efficiency given the accepted function calls, that is, &
1247 263 : &the final refined effective sample size (ESS) divided by the number of accepted function calls."
1248 263 : call self%reportDesc(msg)
1249 :
1250 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1251 :
1252 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1253 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.efficiency.essOverAcceptedRejected"
1254 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1255 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) real(effectiveSampleSize,kind=RK) / real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
1256 : msg = "This is the effective MCMC sampling efficiency given the accepted and rejected function calls, that is, &
1257 263 : &the final refined effective sample size (ESS) divided by the number of (accepted + rejected) function calls."
1258 263 : call self%reportDesc(msg)
1259 :
1260 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1261 :
1262 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1263 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.efficiency.essOverAcceptedRejectedDelayed"
1264 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1265 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) real(effectiveSampleSize,kind=RK) / real(self%Stats%NumFunCall%acceptedRejectedDelayed,kind=RK)
1266 : msg = "This is the effective MCMC sampling efficiency given the accepted, rejected, and delayed-rejection (if any requested) function calls, that is, &
1267 263 : &the final refined effective sample size (ESS) divided by the number of (accepted + rejected + delayed-rejection) function calls."
1268 263 : call self%reportDesc(msg)
1269 :
1270 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1271 :
1272 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1273 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.efficiency.essOverAcceptedRejectedDelayedUnused"
1274 263 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1275 263 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) real(effectiveSampleSize,kind=RK) / real(self%Stats%NumFunCall%acceptedRejectedDelayedUnused,kind=RK)
1276 : msg = "This is the effective MCMC sampling efficiency given the accepted, rejected, delayed-rejection (if any requested), and unused function calls, that is, &
1277 263 : &the final refined effective sample size (ESS) divided by the number of (accepted + rejected + delayed-rejection + unused) function calls."
1278 263 : call self%reportDesc(msg)
1279 :
1280 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1281 :
1282 : !end associate blockEffectiveSampleSize
1283 :
1284 : ! generate output refined sample if requested
1285 :
1286 263 : blockSampleFileGeneration: if (self%SpecBase%SampleSize%val==0_IK) then
1287 :
1288 : call self%note ( prefix = self%brand & ! LCOV_EXCL_LINE
1289 : , outputUnit = self%LogFile%unit & ! LCOV_EXCL_LINE
1290 : , newline = NLC & ! LCOV_EXCL_LINE
1291 2 : , msg = "Skipping the generation of the decorrelated sample and output file, as requested by the user..." )
1292 :
1293 : else blockSampleFileGeneration
1294 :
1295 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1296 : ! sample file generation report
1297 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1298 :
1299 : ! report to the report-file(s)
1300 :
1301 : call self%note ( prefix = self%brand & ! LCOV_EXCL_LINE
1302 : , outputUnit = self%LogFile%unit & ! LCOV_EXCL_LINE
1303 : , newline = NLC & ! LCOV_EXCL_LINE
1304 261 : , msg = "Generating the output " // self%SampleFile%suffix // " file:"//NLC // self%SampleFile%Path%original )
1305 :
1306 :
1307 261 : if (self%Image%isFirst) then
1308 :
1309 : ! print the message for the generating the output sample file on the first image
1310 :
1311 : call self%note ( prefix = self%brand & ! LCOV_EXCL_LINE
1312 : , outputUnit = output_unit & ! LCOV_EXCL_LINE
1313 : , newline = NLC & ! LCOV_EXCL_LINE
1314 : , marginBot = 0_IK & ! LCOV_EXCL_LINE
1315 225 : , msg = "Generating the output " // self%SampleFile%suffix // " file:" )
1316 :
1317 : call self%note ( prefix = self%brand &
1318 : , outputUnit = output_unit &
1319 : , newline = NLC &
1320 : , marginTop = 0_IK &
1321 : , marginBot = 0_IK &
1322 225 : , msg = self%SampleFile%Path%original )
1323 :
1324 : ! print the message for the generating the output sample file on the rest of the images in order
1325 :
1326 : #if defined CAF_ENABLED || defined MPI_ENABLED
1327 225 : if (self%SpecBase%ParallelizationModel%isMultiChain) then
1328 : block
1329 : use String_mod, only: replaceStr, num2str
1330 : integer(IK) :: imageID
1331 54 : do imageID = 2, self%Image%count
1332 : call self%note ( prefix = self%brand &
1333 : , outputUnit = output_unit &
1334 : , newline = NLC &
1335 : , marginTop = 0_IK &
1336 : , marginBot = 0_IK &
1337 54 : , msg = replaceStr( string = self%SampleFile%Path%original, search = "process_1", substitute = "process_"//num2str(imageID) ) )
1338 : end do
1339 : end block
1340 : end if
1341 : #endif
1342 225 : call self%Decor%write(output_unit,0,1)
1343 :
1344 : end if
1345 :
1346 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1347 : ! begin sample file generation
1348 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1349 :
1350 261 : if (self%SpecBase%SampleSize%val/=-1_IK) then
1351 :
1352 85 : if (self%SpecBase%SampleSize%val<0_IK) then
1353 26 : self%SpecBase%SampleSize%abs = abs(self%SpecBase%SampleSize%abs) * self%RefinedChain%Count(self%RefinedChain%numRefinement)%verbose
1354 : end if
1355 :
1356 : ! regenerate the refined sample, this time with the user-specified sample size.
1357 :
1358 : call self%RefinedChain%get ( CFC = self%Chain &
1359 : , Err = self%Err &
1360 : , burninLoc = self%Stats%BurninLoc%compact &
1361 : , refinedChainSize = self%SpecBase%SampleSize%abs &
1362 : , sampleRefinementCount = self%SpecMCMC%SampleRefinementCount%val &
1363 : , sampleRefinementMethod = self%SpecMCMC%SampleRefinementMethod%val &
1364 85 : )
1365 85 : if (self%Err%occurred) then
1366 : ! LCOV_EXCL_START
1367 : self%Err%msg = PROCEDURE_NAME // self%Err%msg
1368 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
1369 : exit blockLeaderPostProcessing
1370 : return
1371 : ! LCOV_EXCL_STOP
1372 : end if
1373 :
1374 : end if
1375 :
1376 : ! open the output sample file
1377 :
1378 261 : self%SampleFile%unit = 3000001 + self%Image%id ! for some unknown reason, if newunit is used, GFortran opens the file as an internal file
1379 : open( unit = self%SampleFile%unit &
1380 : , file = self%SampleFile%Path%original &
1381 : , status = self%SampleFile%status &
1382 : , iostat = self%SampleFile%Err%stat &
1383 : #if defined INTEL_COMPILER_ENABLED && defined OS_IS_WINDOWS
1384 : , SHARED &
1385 : #endif
1386 261 : , position = self%SampleFile%Position%value )
1387 261 : self%Err = self%SampleFile%getOpenErr(self%SampleFile%Err%stat)
1388 261 : if (self%Err%occurred) then
1389 : ! LCOV_EXCL_START
1390 : self%Err%msg = PROCEDURE_NAME//": Error occurred while opening the "//self%name//" "//self%SampleFile%suffix//" file='"//self%SampleFile%Path%original//"'. "
1391 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
1392 : exit blockLeaderPostProcessing
1393 : return
1394 : ! LCOV_EXCL_STOP
1395 : end if
1396 :
1397 : ! determine the sample file's contents' format
1398 :
1399 261 : if (self%SpecBase%OutputColumnWidth%val>0) then
1400 2 : formatStr = "(*(E"//self%SpecBase%OutputColumnWidth%str//"."//self%SpecBase%OutputRealPrecision%str//"E3,:,'"//self%SpecBase%OutputDelimiter%val//"'))"
1401 : else
1402 259 : formatStr = self%SampleFile%format
1403 : end if
1404 :
1405 : ! write to the output sample file
1406 :
1407 261 : call self%RefinedChain%write(self%SampleFile%unit,self%SampleFile%format,formatStr)
1408 261 : close(self%SampleFile%unit)
1409 :
1410 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1411 : ! Compute the statistical properties of the refined sample
1412 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1413 :
1414 261 : if (self%Image%isFirst) then
1415 : call self%note ( prefix = self%brand &
1416 : , outputUnit = output_unit &
1417 : , newline = NLC &
1418 : , marginTop = 2_IK &
1419 225 : , msg = "Computing the statistical properties of the final refined sample..." )
1420 : end if
1421 :
1422 : call self%Decor%writeDecoratedText ( text = "\nThe statistical properties of the final refined sample \n" &
1423 : , marginTop = 1_IK &
1424 : , marginBot = 1_IK &
1425 : , newline = "\n" &
1426 261 : , outputUnit = self%LogFile%unit )
1427 :
1428 261 : self%Stats%Sample%count = self%RefinedChain%Count(self%RefinedChain%numRefinement)%verbose
1429 :
1430 : ! compute the covariance and correlation upper-triangle matrices
1431 :
1432 261 : if (allocated(self%Stats%Sample%Mean)) deallocate(self%Stats%Sample%Mean); allocate(self%Stats%Sample%Mean(ndim))
1433 261 : if (allocated(self%Stats%Sample%CovMat)) deallocate(self%Stats%Sample%CovMat); allocate(self%Stats%Sample%CovMat(ndim,ndim))
1434 261 : if (self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact > ndim + 1_IK) then
1435 122064 : ContiguousChain = self%RefinedChain%LogFuncState(1:ndim,1:self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact)
1436 504 : call getWeiSamCovUppMeanTrans ( np = self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact &
1437 504 : , sumWeight = self%RefinedChain%Count(self%RefinedChain%numRefinement)%verbose &
1438 : , nd = ndim & ! LCOV_EXCL_LINE
1439 : , Point = ContiguousChain & ! LCOV_EXCL_LINE
1440 : , Weight = self%RefinedChain%Weight(1:self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact) & ! LCOV_EXCL_LINE
1441 : , CovMatUpper = self%Stats%Sample%CovMat & ! LCOV_EXCL_LINE
1442 : , Mean = self%Stats%Sample%Mean & ! LCOV_EXCL_LINE
1443 252 : )
1444 2466 : self%Stats%Sample%CorMat = getCorMatUpperFromCovMatUpper(nd=ndim,CovMatUpper=self%Stats%Sample%CovMat)
1445 :
1446 : else
1447 : if (allocated(self%Stats%Sample%CovMat)) deallocate(self%Stats%Sample%CovMat) ! LCOV_EXCL_LINE
1448 59 : allocate(self%Stats%Sample%CovMat(ndim,ndim), source = POSINF_RK)
1449 59 : self%Stats%Sample%CorMat = getEye(ndim,ndim)
1450 : self%Stats%Sample%Mean = getMean ( nd = ndim &
1451 18 : , np = self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact &
1452 : , Point = self%RefinedChain%LogFuncState(1:ndim,1:self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact) & ! LCOV_EXCL_LINE
1453 54 : , Weight = self%RefinedChain%Weight(1:self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact) &
1454 58 : )
1455 : end if
1456 :
1457 : ! transpose the covariance and correlation matrices
1458 :
1459 679 : do i = 1, ndim
1460 418 : self%Stats%Sample%CorMat(i,i) = 1._RK
1461 575 : self%Stats%Sample%CorMat(i+1:ndim,i) = self%Stats%Sample%CorMat(i,i+1:ndim)
1462 836 : self%Stats%Sample%CovMat(i+1:ndim,i) = self%Stats%Sample%CovMat(i,i+1:ndim)
1463 : end do
1464 :
1465 : ! compute the quantiles
1466 :
1467 78034 : ContiguousChain = transpose(self%RefinedChain%LogFuncState(1:ndim,1:self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact)) ! avoid temporary array creation and the warning message in debug mode
1468 261 : if (allocated(self%Stats%Sample%Quantile)) deallocate(self%Stats%Sample%Quantile)
1469 261 : allocate(self%Stats%Sample%Quantile(QPROB%count,ndim))
1470 679 : do i = 1, ndim
1471 3344 : self%Stats%Sample%Quantile(1:QPROB%count,i) = getQuantile ( np = self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact &
1472 : , nq = QPROB%count & ! LCOV_EXCL_LINE
1473 : , SortedQuantileProbability = QPROB%Value & ! LCOV_EXCL_LINE
1474 : , Point = ContiguousChain(:,i) & ! LCOV_EXCL_LINE
1475 : , Weight = self%RefinedChain%Weight(1:self%RefinedChain%Count(self%RefinedChain%numRefinement)%compact) & ! LCOV_EXCL_LINE
1476 836 : , sumWeight = self%RefinedChain%Count(self%RefinedChain%numRefinement)%verbose &
1477 4441 : )
1478 : end do
1479 :
1480 : ! report the refined chain statistics
1481 :
1482 : !formatStr = "(1A" // self%LogFile%maxColWidth%str // ",*(E" // self%LogFile%maxColWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "))"
1483 :
1484 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1485 :
1486 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1487 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.length"
1488 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1489 261 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) self%Stats%Sample%count
1490 261 : msg = "This is the final output refined sample size. "
1491 261 : if (self%SpecBase%SampleSize%val/=-1_IK) then
1492 85 : if (abs(self%SpecBase%SampleSize%abs)<effectiveSampleSize) then
1493 : msg = msg // "The user-requested sample size ("// num2str(self%SpecBase%SampleSize%abs) // ") is smaller &
1494 : &than the potentially-optimal i.i.d. sample size (" // num2str(effectiveSampleSize) // "). The output sample &
1495 : &contains i.i.d. samples, however, the sample-size could have been larger if it had been set to the optimal size. &
1496 11 : &To get the optimal size in the future runs, set sampleSize = -1, or drop it from the input list."
1497 74 : elseif (abs(self%SpecBase%SampleSize%abs)>effectiveSampleSize) then
1498 : msg = msg // "The user-requested sample size ("// num2str(self%SpecBase%SampleSize%abs) // ") is larger than &
1499 : &the potentially-optimal i.i.d. sample size (" // num2str(effectiveSampleSize) // "). The resulting sample &
1500 : &likely contains duplicates and is not independently and identically distributed (i.i.d.).\nTo get the optimal &
1501 74 : &size in the future runs, set sampleSize = -1, or drop it from the input list."
1502 : else ! LCOV_EXCL_LINE
1503 : msg = msg // "How lucky that could be! The user-requested sample size (" // num2str(self%SpecBase%SampleSize%abs) // & ! LCOV_EXCL_LINE
1504 : ") is equal to the potentially-optimal i.i.d. sample size determined by the "//self%name//" sampler." ! LCOV_EXCL_LINE
1505 : end if
1506 : end if
1507 261 : call self%reportDesc(msg)
1508 :
1509 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1510 :
1511 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1512 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.avgStd"
1513 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1514 261 : write(self%LogFile%unit,self%LogFile%format) "variableName", "Mean", "Standard Deviation"
1515 679 : do i = 1, ndim
1516 679 : write(self%LogFile%unit,formatStrReal) trim(adjustl(self%SpecBase%VariableNameList%Val(i))), self%Stats%Sample%Mean(i), sqrt(self%Stats%Sample%CovMat(i,i))
1517 : end do
1518 261 : msg = "This is the Mean and standard deviation table of the final output refined sample."
1519 261 : call self%reportDesc(msg)
1520 :
1521 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1522 :
1523 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1524 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.covmat"
1525 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1526 679 : write(self%LogFile%unit,self%LogFile%format) "", (trim(adjustl(self%SpecBase%VariableNameList%Val(i))),i=1,ndim)
1527 679 : do i = 1, ndim
1528 679 : write(self%LogFile%unit,formatStrReal) trim(adjustl(self%SpecBase%VariableNameList%Val(i))), self%Stats%Sample%CovMat(1:ndim,i)
1529 : end do
1530 261 : msg = "This is the covariance matrix of the final output refined sample."
1531 261 : call self%reportDesc(msg)
1532 :
1533 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1534 :
1535 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1536 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.cormat"
1537 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1538 679 : write(self%LogFile%unit,self%LogFile%format) "", (trim(adjustl(self%SpecBase%VariableNameList%Val(i))),i=1,ndim)
1539 679 : do i = 1, ndim
1540 679 : write(self%LogFile%unit,formatStrReal) trim(adjustl(self%SpecBase%VariableNameList%Val(i))), self%Stats%Sample%CorMat(1:ndim,i)
1541 : end do
1542 261 : msg = "This is the correlation matrix of the final output refined sample."
1543 261 : call self%reportDesc(msg)
1544 :
1545 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1546 :
1547 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1548 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.quantile"
1549 261 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1550 679 : write(self%LogFile%unit,self%LogFile%format) "Quantile", (trim(adjustl(self%SpecBase%VariableNameList%Val(i))),i=1,ndim)
1551 2610 : do iq = 1, QPROB%count
1552 6372 : write(self%LogFile%unit,formatStrReal) trim(adjustl(QPROB%Name(iq))), (self%Stats%Sample%Quantile(iq,i),i=1,ndim)
1553 : end do
1554 261 : msg = "This is the quantiles table of the variables of the final output refined sample."
1555 261 : call self%reportDesc(msg)
1556 :
1557 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1558 : ! Begin inter-chain convergence test in multiChain parallelization mode
1559 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1560 :
1561 : #if defined CAF_ENABLED || defined MPI_ENABLED
1562 :
1563 261 : blockInterChainConvergence: if (self%SpecBase%ParallelizationModel%isMultiChain .and. self%Image%count>1_IK) then
1564 :
1565 : call self%note ( prefix = self%brand &
1566 : , outputUnit = self%LogFile%unit &
1567 : , newline = NLC &
1568 : , marginTop = 1_IK &
1569 : , marginBot = 1_IK &
1570 54 : , msg = "Computing the inter-chain convergence probabilities..." )
1571 54 : if (self%Image%isFirst) then ! print the message on stdout
1572 : call self%note ( prefix = self%brand &
1573 : , outputUnit = output_unit &
1574 : , newline = NLC &
1575 : , marginTop = 1_IK &
1576 : , marginBot = 1_IK &
1577 18 : , msg = "Computing the inter-chain convergence probabilities..." )
1578 : end if
1579 :
1580 54 : call self%Image%sync()
1581 :
1582 : ! read the sample files generated by other images
1583 :
1584 : multiChainConvergenceTest: block
1585 :
1586 : use Sort_mod, only: sortAscending
1587 : use String_mod, only: replaceStr
1588 : use Statistics_mod, only: doSortedKS2
1589 : use ParaMCMCRefinedChain_mod, only: readRefinedChain, RefinedChain_type
1590 54 : type(RefinedChain_type) :: RefinedChainThisImage, RefinedChainThatImage
1591 : integer(IK) :: imageID,indexMinProbKS,imageMinProbKS
1592 54 : real(RK) :: statKS, minProbKS
1593 : real(RK) , allocatable :: ProbKS(:)
1594 54 : character(:), allocatable :: inputSamplePath
1595 :
1596 54 : minProbKS = huge(minProbKS)
1597 : if (allocated(ProbKS)) deallocate(ProbKS) ! LCOV_EXCL_LINE
1598 54 : allocate(ProbKS(0:ndim))
1599 :
1600 : ! read the refined chain on the current image
1601 :
1602 54 : RefinedChainThisImage = readRefinedChain( sampleFilePath=self%SampleFile%Path%original, delimiter=self%SpecBase%OutputDelimiter%val, ndim=ndim, tenabled = .true. )
1603 54 : if (RefinedChainThisImage%Err%occurred) then
1604 : ! LCOV_EXCL_START
1605 : self%Err%occurred = .true.
1606 : self%Err%msg = PROCEDURE_NAME//RefinedChainThisImage%Err%msg
1607 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
1608 : exit blockLeaderPostProcessing
1609 : return
1610 : ! LCOV_EXCL_STOP
1611 : end if
1612 :
1613 : ! sort the refined chain on the current image
1614 :
1615 210 : do i = 0, ndim
1616 0 : call sortAscending ( np = RefinedChainThisImage%Count(RefinedChainThisImage%numRefinement)%verbose &
1617 0 : , Point = RefinedChainThisImage%LogFuncState(1:RefinedChainThisImage%Count(RefinedChainThisImage%numRefinement)%verbose,i) &
1618 : , Err = self%Err &
1619 156 : )
1620 210 : if (self%Err%occurred) then
1621 : ! LCOV_EXCL_START
1622 : self%Err%msg = PROCEDURE_NAME//self%Err%msg
1623 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
1624 : exit blockLeaderPostProcessing
1625 : return
1626 : ! LCOV_EXCL_STOP
1627 : end if
1628 : end do
1629 :
1630 : ! compute and report the KS convergence probabilities for all images
1631 :
1632 54 : formatStr = "(1I" // self%LogFile%maxColWidth%str // ",*(E" // self%LogFile%maxColWidth%str // "." // self%SpecBase%OutputRealPrecision%str // "E3))"
1633 :
1634 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1635 :
1636 54 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1637 54 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.kstest.prob"
1638 54 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1639 210 : write(self%LogFile%unit,self%LogFile%format) "ProcessID",("probKS("//trim(adjustl(self%RefinedChain%ColHeader(i)%record))//")",i=0,ndim)
1640 :
1641 216 : do imageID = 1, self%Image%count
1642 :
1643 216 : if (imageID/=self%Image%id) then
1644 :
1645 : ! read the refined chain on the other image
1646 :
1647 0 : inputSamplePath = replaceStr( string = self%SampleFile%Path%original &
1648 : , search = "process_"//num2str(self%Image%id) &
1649 108 : , substitute = "process_"//num2str(imageID) )
1650 :
1651 318 : RefinedChainThatImage = readRefinedChain( sampleFilePath=inputSamplePath, delimiter=self%SpecBase%OutputDelimiter%val, ndim=ndim, tenabled = .true. )
1652 108 : if (RefinedChainThatImage%Err%occurred) then
1653 : ! LCOV_EXCL_START
1654 : self%Err%occurred = .true.
1655 : self%Err%msg = PROCEDURE_NAME//RefinedChainThatImage%Err%msg
1656 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
1657 : exit blockLeaderPostProcessing
1658 : return
1659 : ! LCOV_EXCL_STOP
1660 : end if
1661 :
1662 420 : do i = 0, ndim
1663 :
1664 : ! sort the refined chain on the other image
1665 :
1666 0 : call sortAscending ( np = RefinedChainThatImage%Count(RefinedChainThatImage%numRefinement)%verbose &
1667 0 : , Point = RefinedChainThatImage%LogFuncState(1:RefinedChainThatImage%Count(RefinedChainThatImage%numRefinement)%verbose,i) &
1668 : , Err = self%Err &
1669 312 : )
1670 312 : if (self%Err%occurred) then
1671 : ! LCOV_EXCL_START
1672 : self%Err%msg = PROCEDURE_NAME//self%Err%msg
1673 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
1674 : exit blockLeaderPostProcessing
1675 : return
1676 : ! LCOV_EXCL_STOP
1677 : end if
1678 :
1679 : ! compute the inter-chain KS probability table
1680 :
1681 0 : call doSortedKS2( np1 = RefinedChainThisImage%Count(RefinedChainThisImage%numRefinement)%verbose &
1682 0 : , np2 = RefinedChainThatImage%Count(RefinedChainThatImage%numRefinement)%verbose &
1683 0 : , SortedPoint1 = RefinedChainThisImage%LogFuncState(1:RefinedChainThisImage%Count(RefinedChainThisImage%numRefinement)%verbose,i) &
1684 0 : , SortedPoint2 = RefinedChainThatImage%LogFuncState(1:RefinedChainThatImage%Count(RefinedChainThatImage%numRefinement)%verbose,i) &
1685 : , statKS = statKS &
1686 0 : , probKS = ProbKS(i) &
1687 312 : )
1688 :
1689 : ! determine the minimum KS
1690 :
1691 420 : if (ProbKS(i)<minProbKS) then
1692 139 : minProbKS = ProbKS(i)
1693 139 : indexMinProbKS = i
1694 139 : imageMinProbKS = imageID
1695 : end if
1696 :
1697 : end do
1698 :
1699 : ! write the inter-chain KS probability table row
1700 :
1701 420 : write(self%LogFile%unit, formatStr) imageID, (ProbKS(i), i = 0, ndim)
1702 :
1703 : end if
1704 :
1705 : end do
1706 : msg = "This is the table pairwise inter-chain Kolmogorov-Smirnov (KS) convergence (similarity) probabilities. &
1707 54 : &Higher KS probabilities are better, indicating less evidence for a lack of convergence."
1708 54 : call self%reportDesc(msg)
1709 :
1710 : ! write the smallest KS probabilities
1711 :
1712 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1713 :
1714 54 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1715 54 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT) "stats.chain.refined.kstest.prob.min"
1716 54 : write(self%LogFile%unit,GENERIC_OUTPUT_FORMAT)
1717 54 : write(self%LogFile%unit,GENERIC_TABBED_FORMAT) minProbKS
1718 108 : msg = "This is the smallest KS-test probability for the inter-chain sampling convergence, which has happened between "//self%RefinedChain%ColHeader(indexMinProbKS)%record// &
1719 54 : " on the chains generated by processes "//num2str(self%Image%id)//" and "//num2str(imageMinProbKS)//"."
1720 54 : call self%reportDesc(msg)
1721 :
1722 : ! report the smallest KS probabilities on stdout
1723 :
1724 54 : if (self%Image%isFirst) call self%note ( prefix = self%brand &
1725 : , outputUnit = output_unit &
1726 : , newline = NLC &
1727 : , marginTop = 2_IK &
1728 : , marginBot = 0_IK &
1729 18 : , msg = "The smallest KS probabilities for the inter-chain sampling convergence:" )
1730 :
1731 690 : do imageID = 1, self%Image%count
1732 162 : if (imageID==self%Image%id) then
1733 : call self%note ( prefix = self%brand &
1734 : , outputUnit = output_unit &
1735 : , newline = NLC &
1736 : , marginTop = 0_IK &
1737 : , marginBot = 0_IK &
1738 108 : , msg = num2str(minProbKS)//" for "//self%RefinedChain%ColHeader(indexMinProbKS)%record//&
1739 : " on the chains generated by processes "//num2str(self%Image%id)//&
1740 54 : " and "//num2str(imageMinProbKS)//"." )
1741 : end if
1742 : #if defined CAF_ENABLED || MPI_ENABLED
1743 162 : call execute_command_line(" ", cmdstat = self%Err%stat)
1744 162 : flush(output_unit)
1745 216 : call self%Image%sync()
1746 : #endif
1747 : end do
1748 :
1749 : end block multiChainConvergenceTest
1750 :
1751 : end if blockInterChainConvergence
1752 : #endif
1753 :
1754 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1755 : ! End inter-chain convergence test in multiChain parallelization mode
1756 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1757 :
1758 : end if blockSampleFileGeneration
1759 :
1760 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1761 : ! End of generating the i.i.d. sample statistics and output file (if requested)
1762 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1763 :
1764 263 : call self%Decor%write(self%LogFile%unit,1,1)
1765 :
1766 : ! Mission accomplished.
1767 :
1768 263 : call self%note( prefix = self%brand, outputUnit = self%LogFile%unit, newline = "\n", msg = "Mission Accomplished." )
1769 263 : if (output_unit/=self%LogFile%unit .and. self%Image%isFirst) then
1770 227 : call self%Decor%write(output_unit,1,1)
1771 227 : call self%note( prefix = self%brand, outputUnit = output_unit, newline = "\n", msg = "Mission Accomplished." )
1772 227 : call self%Decor%write(output_unit,1,1)
1773 : end if
1774 :
1775 263 : close(self%TimeFile%unit)
1776 263 : close(self%LogFile%unit)
1777 :
1778 : end if blockLeaderPostProcessing
1779 :
1780 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
1781 681 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
1782 : #endif
1783 : if (self%Err%occurred) return ! LCOV_EXCL_LINE
1784 :
1785 : !nullify(self%Proposal)
1786 : #if defined CAF_ENABLED || defined MPI_ENABLED
1787 681 : if (self%SpecBase%MpiFinalizeRequested%val) call self%Image%finalize()
1788 : #endif
1789 :
1790 2334 : end subroutine runSampler
1791 :
1792 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1793 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1794 :
1795 : !end submodule Setup_smod
1796 :
1797 : #undef ParaDXXXProposalAbstract_mod
1798 : #undef ParaDXXX_type
1799 : #undef ParaDXXX
|
