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 `Kernel_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_type ParaDRAM_type
54 : #define ParaDXXXProposalAbstract_mod ParaDRAMProposalAbstract_mod
55 :
56 : #elif defined PARADISE
57 :
58 : #define ParaDXXX_type ParaDISE_type
59 : #define ParaDXXXProposalAbstract_mod ParaDISEProposalAbstract_mod
60 :
61 : #else
62 :
63 : #error "Unrecognized sampler in ParaDXXX_mod.inc.f90"
64 :
65 : #endif
66 :
67 : use, intrinsic :: iso_fortran_env, only: output_unit
68 : !use Constants_mod, only: IK, RK ! gfortran 9.3 compile crashes with this line
69 : #if defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED || ( (defined MATLAB_ENABLED || defined PYTHON_ENABLED || defined R_ENABLED) && !defined CAF_ENABLED && !defined MPI_ENABLED )
70 : use ParaDXXXProposalAbstract_mod, only: ProposalErr
71 : #endif
72 :
73 : #if defined MPI_ENABLED
74 : use mpi
75 : #endif
76 :
77 : implicit none
78 :
79 : character(*), parameter :: SUBMODULE_NAME = MODULE_NAME // "@Kernel_smod"
80 :
81 : type :: SumAccRateSinceStart_type
82 : real(RK) :: acceptedRejected
83 : real(RK) :: acceptedRejectedDelayed
84 : end type SumAccRateSinceStart_type
85 :
86 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
87 :
88 : contains
89 :
90 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 :
92 : !> \brief
93 : !> This procedure is a method of [ParaDRAM_type](@ref paradram_type) and [ParaDISE_type](@ref paradise_type) classes.
94 : !> Run the sampler and return the sampled states.
95 : !>
96 : !> @param[inout] self : An object of class [ParaDRAM_type](@ref paradram_type) or [ParaDISE_type](@ref paradise_type).
97 : !> @param[in] getLogFunc : The target objective function that is to be sampled.
98 : !>
99 : !> \remark
100 : !> This procedure requires preprocessing.
101 687 : module subroutine runKernel ( self &
102 : , getLogFunc &
103 : )
104 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
105 : !DEC$ ATTRIBUTES DLLEXPORT :: runKernel
106 : #endif
107 : use Constants_mod, only: RK, IK, NEGINF_RK, NLC, LOGTINY_RK, NEGLOGINF_RK
108 : use ParaMonteLogFunc_mod, only: getLogFunc_proc
109 : use Math_mod, only: getLogSubExp
110 : use Decoration_mod, only: write
111 : use String_mod, only: num2str
112 :
113 : implicit none
114 :
115 : character(*), parameter :: PROCEDURE_NAME = SUBMODULE_NAME//"@runKernel()"
116 : integer(IK) , parameter :: CHAIN_RESTART_OFFSET = 2_IK
117 :
118 : class(ParaDXXX_type), intent(inout) :: self
119 : procedure(getLogFunc_proc) :: getLogFunc
120 :
121 : #if defined CAF_ENABLED
122 : real(RK) , allocatable :: co_AccRate(:)[:]
123 : real(RK) , allocatable :: co_LogFuncState(:,:)[:] ! (0:nd,-1:delayedRejectionCount), -1 corresponds to the current accepted state
124 : integer(IK) , save :: co_proposalFound_samplerUpdateOccurred(2)[*] ! merging these scalars would reduce the MPI communication overhead cost: co_proposalFound, co_samplerUpdateOccurred, co_counterDRS, 0 means false, 1 means true
125 : #else
126 : real(RK) , allocatable :: co_LogFuncState(:,:)
127 : real(RK) , allocatable :: co_AccRate(:)
128 : integer(IK) , save :: co_proposalFound_samplerUpdateOccurred(2) ! merging these scalars would reduce the MPI communication overhead cost: co_proposalFound, co_samplerUpdateOccurred, co_counterDRS, 0 means false, 1 means true
129 : #endif
130 : type(SumAccRateSinceStart_type) :: SumAccRateSinceStart ! used to figure out the average acceptance ratio for the entire chain.
131 : integer(IK) :: numFunCallAcceptedLastAdaptation ! number of function calls accepted at Last proposal adaptation occurrence
132 : integer(IK) :: counterAUP ! counter for adaptiveUpdatePeriod
133 : integer(IK) :: counterAUC ! counter for adaptiveUpdateCount
134 : integer(IK) :: counterPRP ! counter for progressReportPeriod
135 : integer(IK) :: counterDRS ! counter for Delayed Rejection Stages
136 : integer(IK) :: lastStateWeight ! This is used for passing the most recent verbose chain segment to the adaptive updater of the sampler
137 : integer(IK) :: currentStateWeight ! counter for SampleWeight, used only in in restart mode
138 : integer(IK) :: numFunCallAcceptedPlusOne ! counter for SampleWeight, used only in in restart mode
139 : integer(IK) :: numFunCallAcceptedRejectedLastReport ! used for progress-report: must be initialized to zero upon entry to the procedure
140 687 : real(RK) :: timeElapsedUntilLastReportInSeconds ! used for progress-report: must be initialized to zero upon entry to the procedure
141 687 : real(RK) :: inverseProgressReportPeriod ! used for progress-report: inverse of progressReportPeriod
142 687 : real(RK) :: sumAccRateLastReport ! used for progress-report: must be initialized to zero upon entry to the procedure
143 687 : real(RK) :: uniformRnd ! used for random number generation.
144 687 : real(RK) :: meanAccRateSinceStart ! used for restart file read
145 687 : real(RK) :: logFuncDiff ! The difference between the log of the old and the new states. Used to avoid underflow.
146 : !real(RK) :: adaptationMeasureDummy
147 687 : real(RK) :: maxLogFuncRejectedProposal
148 : logical :: samplerUpdateIsGreedy
149 : logical :: samplerUpdateSucceeded
150 : logical :: delayedRejectionRequested
151 : logical :: noDelayedRejectionRequested
152 : real(RK) , allocatable :: AdaptationMeasure(:), adaptationMeasurePlaceHolder(:)
153 : integer(IK) :: adaptationMeasureLen
154 : integer(IK) :: acceptedRejectedDelayedUnusedRestartMode
155 : integer(IK) :: imageID, dumint
156 : integer(IK) :: nd
157 : integer(IK), parameter :: STDOUT_SEGLEN = 25
158 : character(4*STDOUT_SEGLEN+2*3+1) :: txt
159 : #if defined CAF_ENABLED || defined MPI_ENABLED
160 : integer(IK) :: imageStartID, imageEndID
161 : integer(IK) :: proposalFoundSinglChainMode ! used in singlChain delayed rejection. zero if the proposal is not accepted. 1 if the proposal is accepted.
162 : #if defined MPI_ENABLED
163 : integer(IK) :: proposalFoundSinglChainModeReduced ! the reduced value by summing proposalFoundSinglChainModeReduced over all images.
164 : real(RK) , allocatable :: AccRateMatrix(:,:) ! matrix of size (-1:self%SpecDRAM%DelayedRejectionCount%val,1:self%Image%count)
165 : integer(IK) :: ndPlusOne
166 : integer(IK) :: ierrMPI
167 : integer(IK) :: delayedRejectionCountPlusTwo
168 : #endif
169 687 : self%Stats%avgCommTimePerFunCall = 0._RK ! Until the reporting time, this is in reality, sumCommTimePerFunCall. This is meaningful only in singlChain parallelism
170 687 : self%Stats%NumFunCall%acceptedRejectedDelayedUnused = self%Image%count ! used only in singlChain parallelism, and relevant only on the first image
171 : #endif
172 687 : acceptedRejectedDelayedUnusedRestartMode = 0_IK ! used to compute more accurate timings in the restart mode
173 687 : self%Stats%avgTimePerFunCalInSec = 0._RK
174 687 : numFunCallAcceptedRejectedLastReport = 0_IK
175 687 : timeElapsedUntilLastReportInSeconds = 0._RK
176 687 : inverseProgressReportPeriod = 1._RK/real(self%SpecBase%ProgressReportPeriod%val,kind=RK) ! this remains a constant except for the last the last report of the simulation
177 687 : sumAccRateLastReport = 0._RK
178 687 : nd = self%nd%val
179 :
180 681 : if (allocated(co_AccRate)) deallocate(co_AccRate)
181 687 : if (allocated(co_LogFuncState)) deallocate(co_LogFuncState)
182 : #if defined CAF_ENABLED
183 687 : allocate(co_LogFuncState(0:nd,-1:self%SpecDRAM%DelayedRejectionCount%val)[*])
184 687 : allocate(co_AccRate(-1:self%SpecDRAM%DelayedRejectionCount%val)[*]) ! the negative element will contain counterDRS
185 : #else
186 : allocate(co_LogFuncState(0:nd,-1:self%SpecDRAM%DelayedRejectionCount%val))
187 : allocate(co_AccRate(-1:self%SpecDRAM%DelayedRejectionCount%val)) ! the negative element will contain counterDRS
188 : #endif
189 687 : co_AccRate(-1) = 0._RK ! the real-value counterDRS, indicating the initial delayed rejection stage at which the first point is sampled
190 687 : co_AccRate(0) = 1._RK ! initial acceptance rate for the first zeroth DR stage.
191 849 : co_AccRate(1:self%SpecDRAM%DelayedRejectionCount%val) = 0._RK ! indicates the very first proposal acceptance on image 1
192 :
193 : #if defined MPI_ENABLED
194 : if (allocated(AccRateMatrix)) deallocate(AccRateMatrix)
195 : allocate(AccRateMatrix(-1:self%SpecDRAM%DelayedRejectionCount%val,1:self%Image%count)) ! the negative element will contain counterDRS
196 : AccRateMatrix = 0._RK ! -huge(1._RK) ! debug
197 : AccRateMatrix(0,1:self%Image%count) = 1._RK ! initial acceptance rate for the first zeroth DR stage.
198 : ndPlusOne = nd + 1_IK
199 : delayedRejectionCountPlusTwo = self%SpecDRAM%DelayedRejectionCount%val + 2_IK
200 : #endif
201 :
202 687 : delayedRejectionRequested = self%SpecDRAM%DelayedRejectionCount%val > 0_IK
203 687 : noDelayedRejectionRequested = .not. delayedRejectionRequested
204 687 : if (delayedRejectionRequested) then
205 48 : self%Stats%NumFunCall%acceptedRejectedDelayed = 0_IK ! Markov Chain counter
206 48 : SumAccRateSinceStart%acceptedRejectedDelayed = 0._RK ! sum of acceptance rate
207 : end if
208 :
209 : !adaptationMeasure = 0._RK ! needed for the first output
210 687 : SumAccRateSinceStart%acceptedRejected = 0._RK ! sum of acceptance rate
211 687 : self%Stats%NumFunCall%acceptedRejected = 0_IK ! Markov Chain counter
212 687 : counterAUC = 0_IK ! counter for padaptiveUpdateCount.
213 687 : counterPRP = 0_IK ! counter for progressReportPeriod.
214 687 : counterAUP = 0_IK ! counter for adaptiveUpdatePeriod.
215 687 : self%Stats%NumFunCall%accepted = 0_IK ! Markov Chain acceptance counter.
216 687 : samplerUpdateSucceeded = .true. ! needed to set up lastStateWeight and numFunCallAcceptedLastAdaptation for the first accepted proposal
217 687 : numFunCallAcceptedLastAdaptation = 0_IK
218 687 : lastStateWeight = -huge(lastStateWeight)
219 687 : meanAccRateSinceStart = 1._RK ! needed for the first restart output in fresh run.
220 :
221 687 : call self%Timer%tic()
222 :
223 687 : blockDryRunSetup: if (self%isFreshRun) then
224 :
225 627 : adaptationMeasureLen = 100_IK
226 :
227 627 : allocate(self%Chain%ProcessID ( self%SpecMCMC%ChainSize%val))
228 627 : allocate(self%Chain%DelRejStage ( self%SpecMCMC%ChainSize%val))
229 627 : allocate(self%Chain%MeanAccRate ( self%SpecMCMC%ChainSize%val))
230 627 : allocate(self%Chain%Adaptation ( self%SpecMCMC%ChainSize%val))
231 627 : allocate(self%Chain%BurninLoc ( self%SpecMCMC%ChainSize%val))
232 627 : allocate(self%Chain%Weight ( self%SpecMCMC%ChainSize%val))
233 627 : allocate(self%Chain%State (nd,self%SpecMCMC%ChainSize%val))
234 627 : allocate(self%Chain%LogFunc ( self%SpecMCMC%ChainSize%val))
235 :
236 : else blockDryRunSetup
237 :
238 : ! load the existing Chain file into self%Chain components
239 :
240 : call self%Chain%get ( chainFilePath = self%ChainFile%Path%original & ! LCOV_EXCL_LINE
241 : , chainFileForm = self%SpecBase%ChainFileFormat%val & ! LCOV_EXCL_LINE
242 : , Err = self%Err & ! LCOV_EXCL_LINE
243 : , targetChainSize = self%SpecMCMC%ChainSize%val & ! LCOV_EXCL_LINE
244 : , lenHeader = self%Chain%lenHeader & ! LCOV_EXCL_LINE
245 : , ndim = nd & ! LCOV_EXCL_LINE
246 : , delimiter = self%SpecBase%OutputDelimiter%val & ! LCOV_EXCL_LINE
247 60 : )
248 60 : if (self%Err%occurred) then
249 : ! LCOV_EXCL_START
250 : self%Err%msg = PROCEDURE_NAME//self%Err%msg
251 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
252 : return
253 : ! LCOV_EXCL_STOP
254 : end if
255 :
256 60 : if (self%Chain%Count%compact<=CHAIN_RESTART_OFFSET) then
257 0 : self%isFreshRun = .true.
258 0 : self%isDryRun = .not. self%isFreshRun
259 : end if
260 :
261 60 : call self%Image%sync()
262 :
263 60 : blockLeaderSetup: if (self%Image%isLeader) then
264 :
265 : ! set up the chain file
266 :
267 : block
268 :
269 : use System_mod, only: RandomFileName_type, copyFile, removeFile
270 32 : type(RandomFileName_type) :: RFN
271 :
272 : ! create a copy of the chain file, just for the sake of not losing the simulation results
273 :
274 32 : RFN = RandomFileName_type(dir = "", key = self%ChainFile%Path%original//".rst", ext="") ! temporary_restart_copy
275 32 : call copyFile(pathOld=self%ChainFile%Path%original,pathNew=RFN%path,isUnixShell=self%OS%Shell%isUnix,Err=self%err)
276 32 : if (self%Err%occurred) then
277 : ! LCOV_EXCL_START
278 : self%Err%msg = PROCEDURE_NAME//self%Err%msg
279 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
280 : exit blockLeaderSetup
281 : return
282 : ! LCOV_EXCL_STOP
283 : end if
284 :
285 : ! reopen the chain file to resume the simulation
286 :
287 : open( newunit = self%ChainFile%unit & ! LCOV_EXCL_LINE
288 : , file = self%ChainFile%Path%original & ! LCOV_EXCL_LINE
289 : , form = self%ChainFile%Form%value & ! LCOV_EXCL_LINE
290 : , status = self%ChainFile%status & ! LCOV_EXCL_LINE
291 : , iostat = self%ChainFile%Err%stat & ! LCOV_EXCL_LINE
292 : #if defined INTEL_COMPILER_ENABLED && defined OS_IS_WINDOWS
293 : , SHARED & ! LCOV_EXCL_LINE
294 : #endif
295 32 : , position = self%ChainFile%Position%value )
296 32 : self%Err = self%ChainFile%getOpenErr(self%ChainFile%Err%stat)
297 32 : if (self%Err%occurred) then
298 : ! LCOV_EXCL_START
299 : self%Err%msg = PROCEDURE_NAME//": Error occurred while opening "//self%name//" "//self%ChainFile%suffix//" file='"//self%ChainFile%Path%original//"'."
300 : call self%abort( Err = self%Err, prefix = self%brand, newline = NLC, outputUnit = self%LogFile%unit )
301 : exit blockLeaderSetup
302 : return
303 : ! LCOV_EXCL_STOP
304 : end if
305 :
306 : ! rewrite the chain file
307 :
308 : call self%Chain%writeChainFile ( ndim = nd & ! LCOV_EXCL_LINE
309 : , compactStartIndex = 1_IK & ! LCOV_EXCL_LINE
310 : , compactEndIndex = self%Chain%Count%compact-CHAIN_RESTART_OFFSET & ! LCOV_EXCL_LINE
311 : , chainFileUnit = self%ChainFile%unit & ! LCOV_EXCL_LINE
312 : , chainFileForm = self%SpecBase%ChainFileFormat%val & ! LCOV_EXCL_LINE
313 : , chainFileFormat = self%ChainFile%format & ! LCOV_EXCL_LINE
314 : , adaptiveUpdatePeriod = self%SpecDRAM%AdaptiveUpdatePeriod%val & ! LCOV_EXCL_LINE
315 32 : )
316 :
317 : ! remove the temporary copy of the chain file
318 :
319 64 : call removeFile(RFN%path, self%Err) ! Passing the optional Err argument, handles exceptions should any occur.
320 :
321 : end block
322 :
323 18368 : adaptationMeasureLen = maxval(self%Chain%Weight(1:self%Chain%Count%compact-CHAIN_RESTART_OFFSET))
324 :
325 : end if blockLeaderSetup
326 :
327 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
328 60 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
329 : #endif
330 60 : if (self%Err%occurred) return
331 :
332 : end if blockDryRunSetup
333 :
334 687 : if (self%Image%isLeader) then
335 269 : if (allocated(AdaptationMeasure)) deallocate(AdaptationMeasure); allocate(AdaptationMeasure(adaptationMeasureLen))
336 : end if
337 :
338 687 : if (self%isFreshRun) then ! this must be done separately from the above blockDryRunSetup
339 :
340 627 : self%Chain%BurninLoc(1) = 1_IK
341 :
342 1569 : co_LogFuncState(1:nd,0) = self%SpecMCMC%StartPointVec%Val ! proposal state
343 627 : call self%Timer%toc()
344 627 : co_LogFuncState(0,0) = getLogFunc( nd, co_LogFuncState(1:nd,0) ) ! proposal logFunc
345 627 : call self%Timer%toc(); self%Stats%avgTimePerFunCalInSec = self%Stats%avgTimePerFunCalInSec + self%Timer%Time%delta
346 :
347 : else
348 :
349 180 : co_LogFuncState(1:nd,0) = self%Chain%State(1:nd,1) ! proposal logFunc
350 60 : co_LogFuncState(0,0) = self%Chain%LogFunc(1) ! proposal logFunc
351 :
352 : end if
353 :
354 2436 : co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,0) ! set current logFunc and State equal to the first proposal
355 687 : self%Stats%LogFuncMode%val = -huge(self%Stats%LogFuncMode%val)
356 687 : self%Stats%LogFuncMode%Loc%compact = 0_IK
357 :
358 687 : if (self%Image%isFirst) then
359 : ! LCOV_EXCL_START
360 : txt = repeat(" ",STDOUT_SEGLEN) &
361 : // "Accepted/Total Func. Call " &
362 : // "Dynamic/Overall Acc. Rate " &
363 : // "Elapsed/Remained Time [s] "
364 : call write(string=txt)
365 : txt = repeat(" ",STDOUT_SEGLEN) &
366 : // repeat("=",STDOUT_SEGLEN) // " " &
367 : // repeat("=",STDOUT_SEGLEN) // " " &
368 : // repeat("=",STDOUT_SEGLEN) // " "
369 : call write(string=txt)
370 : !call execute_command_line(" ")
371 : flush(output_unit)
372 : ! LCOV_EXCL_STOP
373 : end if
374 :
375 : #if defined CAF_ENABLED || defined MPI_ENABLED
376 687 : if (self%SpecBase%ParallelizationModel%isSinglChain) then
377 627 : imageStartID = 1_IK
378 627 : imageEndID = self%Image%count
379 : else ! isMultiChain
380 60 : imageStartID = self%Image%id
381 60 : imageEndID = self%Image%id
382 : end if
383 : #else
384 : imageID = 1_IK ! needed even in the case of serial run to assign a proper value to self%Chain%ProcessID
385 : #endif
386 :
387 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% start of loopMarkovChain %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
388 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
389 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
390 :
391 383724 : loopMarkovChain: do
392 :
393 384411 : co_proposalFound_samplerUpdateOccurred(2) = 0_IK ! at each iteration assume no samplerUpdateOccurred, unless it occurs
394 :
395 : #if defined CAF_ENABLED || defined MPI_ENABLED
396 384411 : blockLeaderImage: if (self%Image%isLeader) then
397 : #endif
398 :
399 194091 : co_proposalFound_samplerUpdateOccurred(1) = 0_IK ! co_proposalFound = .false.
400 194091 : samplerUpdateIsGreedy = counterAUC < self%SpecDRAM%GreedyAdaptationCount%val
401 :
402 : #if defined CAF_ENABLED || defined MPI_ENABLED
403 383848 : loopOverImages: do imageID = imageStartID, imageEndID
404 : #if defined CAF_ENABLED
405 283911 : call self%Timer%toc()
406 283911 : if (imageID/=self%Image%id) co_AccRate(-1:self%SpecDRAM%DelayedRejectionCount%val) = co_AccRate(-1:self%SpecDRAM%DelayedRejectionCount%val)[imageID] ! happens only in isSinglChain=TRUE
407 283911 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
408 : #elif defined MPI_ENABLED
409 : if (imageID/=self%Image%id) co_AccRate(-1:self%SpecDRAM%DelayedRejectionCount%val) = AccRateMatrix(-1:self%SpecDRAM%DelayedRejectionCount%val,imageID) ! happens only in isSinglChain=TRUE
410 : #endif
411 : #endif
412 :
413 283911 : counterDRS = nint(co_AccRate(-1),kind=IK)
414 283911 : if (counterDRS > -1_IK) co_proposalFound_samplerUpdateOccurred(1) = 1_IK ! co_proposalFound = .true.
415 :
416 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% blockProposalAccepted %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
417 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
418 :
419 : ! On the very first iteration, this block is (and must be) executed for imageID==1,
420 : ! since it is for the first (starting) point, which is assumed to have been accepted
421 : ! as the first point by the first coarray imageID.
422 :
423 283911 : blockProposalAccepted: if ( co_proposalFound_samplerUpdateOccurred(1) == 1_IK ) then ! co_proposalAccepted = .true.
424 :
425 94150 : currentStateWeight = 0_IK
426 :
427 : ! communicate the accepted logFunc and State from the winning image to the leader/all images: co_LogFuncState
428 :
429 : #if defined MPI_ENABLED
430 : if (self%SpecBase%ParallelizationModel%isSinglChain) then
431 : call self%Timer%toc()
432 : ! broadcast winning image to all processes
433 : call mpi_bcast ( imageID & ! buffer
434 : , 1 & ! count
435 : , mpi_integer & ! datatype
436 : , 0 & ! root: broadcasting rank
437 : , mpi_comm_world & ! comm
438 : , ierrMPI & ! ierr
439 : )
440 : ! broadcast co_LogFuncState from the winning image to all others
441 : call mpi_bcast ( co_LogFuncState(0:nd,-1) & ! buffer
442 : , ndPlusOne & ! count
443 : , mpi_double_precision & ! datatype
444 : , imageID - 1_IK & ! root: broadcasting rank
445 : , mpi_comm_world & ! comm
446 : , ierrMPI & ! ierr
447 : )
448 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
449 : end if
450 : #elif defined CAF_ENABLED
451 94150 : if (imageID/=self%Image%id) then ! Avoid remote connection for something that is available locally.
452 26243 : call self%Timer%toc()
453 26243 : co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,-1)[imageID]
454 26243 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
455 : end if
456 : #endif
457 :
458 : ! Note: after every adaptive update of the sampler, counterAUP is reset to 0.
459 94150 : if (counterAUP==0_IK .and. samplerUpdateSucceeded) then
460 14732 : numFunCallAcceptedLastAdaptation = numFunCallAcceptedLastAdaptation + 1_IK
461 14732 : lastStateWeight = 0_IK
462 : end if
463 :
464 94150 : blockFreshDryRun: if (self%isFreshRun) then
465 :
466 75750 : call writeOutput()
467 :
468 75750 : self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
469 :
470 75750 : self%Chain%ProcessID(self%Stats%NumFunCall%accepted) = imageID
471 75750 : self%Chain%DelRejStage(self%Stats%NumFunCall%accepted) = counterDRS
472 75750 : self%Chain%Adaptation(self%Stats%NumFunCall%accepted) = 0._RK ! adaptationMeasure
473 75750 : self%Chain%Weight(self%Stats%NumFunCall%accepted) = 0_IK
474 75750 : self%Chain%LogFunc(self%Stats%NumFunCall%accepted) = co_LogFuncState(0,-1)
475 196818 : self%Chain%State(1:nd,self%Stats%NumFunCall%accepted) = co_LogFuncState(1:nd,-1)
476 :
477 : ! find the burnin point
478 :
479 151500 : self%Chain%BurninLoc(self%Stats%NumFunCall%accepted) = getBurninLoc ( lenLogFunc = self%Stats%NumFunCall%accepted &
480 : , refLogFunc = self%Stats%LogFuncMode%val &
481 303000 : , LogFunc = self%Chain%LogFunc(1:self%Stats%NumFunCall%accepted) &
482 75750 : )
483 :
484 : else blockFreshDryRun ! in restart mode: determine the correct value of co_proposalFound_samplerUpdateOccurred(1)
485 :
486 18400 : numFunCallAcceptedPlusOne = self%Stats%NumFunCall%accepted + 1_IK
487 18400 : if (numFunCallAcceptedPlusOne==self%Chain%Count%compact) then
488 32 : self%isFreshRun = .true.
489 32 : call writeOutput()
490 32 : self%isDryRun = .not. self%isFreshRun
491 32 : self%Chain%Weight(numFunCallAcceptedPlusOne) = 0_IK
492 32 : SumAccRateSinceStart%acceptedRejected = self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted) * real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
493 32 : if (delayedRejectionRequested) SumAccRateSinceStart%acceptedRejectedDelayed = self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted) * real(self%Stats%NumFunCall%acceptedRejectedDelayed,kind=RK)
494 : end if
495 18400 : self%Stats%NumFunCall%accepted = numFunCallAcceptedPlusOne
496 18400 : numFunCallAcceptedPlusOne = self%Stats%NumFunCall%accepted + 1_IK
497 :
498 : end if blockFreshDryRun
499 :
500 94150 : if (self%Stats%LogFuncMode%val < self%Chain%LogFunc(self%Stats%NumFunCall%accepted)) then
501 445 : self%Stats%LogFuncMode%val = self%Chain%LogFunc(self%Stats%NumFunCall%accepted)
502 445 : self%Stats%LogFuncMode%Loc%compact = self%Stats%NumFunCall%accepted
503 : end if
504 :
505 94150 : SumAccRateSinceStart%acceptedRejected = SumAccRateSinceStart%acceptedRejected + co_AccRate(counterDRS)
506 :
507 : else blockProposalAccepted
508 :
509 189761 : counterDRS = self%SpecDRAM%DelayedRejectionCount%val
510 189761 : SumAccRateSinceStart%acceptedRejected = SumAccRateSinceStart%acceptedRejected + co_AccRate(counterDRS)
511 :
512 : end if blockProposalAccepted
513 :
514 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
515 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% blockProposalAccepted %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
516 :
517 283911 : counterAUP = counterAUP + 1_IK
518 283911 : counterPRP = counterPRP + 1_IK
519 :
520 : ! It is critical for this if block to occur before updating `self%Stats%NumFunCall%acceptedRejected` otherwise,
521 : ! `numFunCallAcceptedRejectedLastReport` will be updated to `self%Stats%NumFunCall%acceptedRejected`
522 : ! which can lead to "Floating-point exception - erroneous arithmetic operation" in the computation
523 : ! of `inverseProgressReportPeriod` when `blockLastSample` happens to be activated.
524 283911 : if (counterPRP == self%SpecBase%ProgressReportPeriod%val) then
525 145 : counterPRP = 0_IK
526 145 : call reportProgress()
527 : end if
528 :
529 283911 : currentStateWeight = currentStateWeight + 1_IK
530 283911 : self%Stats%NumFunCall%acceptedRejected = self%Stats%NumFunCall%acceptedRejected + 1_IK
531 :
532 283911 : if (delayedRejectionRequested) then
533 88526 : SumAccRateSinceStart%acceptedRejectedDelayed = SumAccRateSinceStart%acceptedRejectedDelayed + sum(co_AccRate(0:counterDRS))
534 18941 : self%Stats%NumFunCall%acceptedRejectedDelayed = self%Stats%NumFunCall%acceptedRejectedDelayed + counterDRS + 1_IK
535 : end if
536 :
537 283911 : if (self%isFreshRun) then ! these are used for adaptive proposal updating, so they have to be set on every accepted or rejected iteration (excluding delayed rejections)
538 222935 : self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted) = SumAccRateSinceStart%acceptedRejected / real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
539 222935 : self%Chain%Weight(self%Stats%NumFunCall%accepted) = self%Chain%Weight(self%Stats%NumFunCall%accepted) + 1_IK
540 222935 : if (adaptationMeasureLen<self%Chain%Weight(self%Stats%NumFunCall%accepted)) then
541 8 : allocate(adaptationMeasurePlaceHolder(2*adaptationMeasureLen))
542 12808 : adaptationMeasurePlaceHolder(1:adaptationMeasureLen) = AdaptationMeasure
543 8 : call move_alloc(from=adaptationMeasurePlaceHolder,to=AdaptationMeasure)
544 8 : adaptationMeasureLen = 2_IK * adaptationMeasureLen
545 : end if
546 : else
547 : #if defined CAF_ENABLED || defined MPI_ENABLED
548 60976 : acceptedRejectedDelayedUnusedRestartMode = self%Stats%NumFunCall%acceptedRejectedDelayedUnused
549 : #else
550 : acceptedRejectedDelayedUnusedRestartMode = self%Stats%NumFunCall%acceptedRejectedDelayed
551 : #endif
552 : end if
553 :
554 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% last output write %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
555 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
556 :
557 : ! in paradise mode, it is imperative to finish the simulation before any further redundant sampler updates occurs.
558 : ! This is the reason why blockLastSample appears before blockSamplerAdaptation.
559 283911 : blockLastSample: if (self%Stats%NumFunCall%accepted==self%SpecMCMC%ChainSize%val) then !co_missionAccomplished = .true.
560 :
561 : ! on 3 images Windows, substituting co_missionAccomplished with the following leads to 10% less communication overhead for 1D Gaussian example
562 : ! on 3 images Linux , substituting co_missionAccomplished with the following leads to 16% less communication overhead for 1D Gaussian example
563 : ! on 5 images Linux , substituting co_missionAccomplished with the following leads to 11% less communication overhead for 1D Gaussian example
564 :
565 263 : co_proposalFound_samplerUpdateOccurred(1) = -1_IK ! equivalent to co_missionAccomplished = .true.
566 263 : inverseProgressReportPeriod = 1._RK / (self%Stats%NumFunCall%acceptedRejected-numFunCallAcceptedRejectedLastReport)
567 :
568 263 : if (self%isFreshRun) then
569 263 : call writeOutput()
570 263 : flush(self%ChainFile%unit)
571 : end if
572 :
573 263 : call reportProgress()
574 :
575 : #if defined CAF_ENABLED || defined MPI_ENABLED
576 263 : exit loopOverImages
577 : #endif
578 : end if blockLastSample
579 :
580 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
581 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% last output write %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
582 :
583 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Proposal Adaptation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
584 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
585 :
586 283648 : blockSamplerAdaptation: if ( counterAUC < self%SpecDRAM%AdaptiveUpdateCount%val .and. counterAUP == self%SpecDRAM%AdaptiveUpdatePeriod%val ) then
587 :
588 124416 : co_proposalFound_samplerUpdateOccurred(2) = 1_IK ! istart = numFunCallAcceptedLastAdaptation ! = max( numFunCallAcceptedLastAdaptation , self%Chain%BurninLoc(self%Stats%NumFunCall%accepted) ) ! this is experimental
589 :
590 : ! the order in the following two MUST be preserved as occasionally self%Stats%NumFunCall%accepted = numFunCallAcceptedLastAdaptation
591 :
592 124416 : dumint = self%Chain%Weight(self%Stats%NumFunCall%accepted) ! needed for the restart mode, not needed in the fresh run
593 124416 : if (self%Stats%NumFunCall%accepted==numFunCallAcceptedLastAdaptation) then ! no new point has been accepted since last time
594 49522 : self%Chain%Weight(numFunCallAcceptedLastAdaptation) = currentStateWeight - lastStateWeight
595 : #if (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED || defined DEBUG_ENABLED || defined TESTING_ENABLED) && !defined CAF_ENABLED && !defined MPI_ENABLED
596 : if (mod(self%Chain%Weight(numFunCallAcceptedLastAdaptation),self%SpecDRAM%AdaptiveUpdatePeriod%val)/=0) then
597 : write(output_unit,"(*(g0,:,' '))" ) PROCEDURE_NAME//": Internal error occurred: ", self%SpecDRAM%AdaptiveUpdatePeriod%val &
598 : , self%Chain%Weight(numFunCallAcceptedLastAdaptation), currentStateWeight, lastStateWeight
599 : !call execute_command_line(" ")
600 : flush(output_unit)
601 : error stop
602 : end if
603 : #endif
604 : else
605 74894 : self%Chain%Weight(numFunCallAcceptedLastAdaptation) = self%Chain%Weight(numFunCallAcceptedLastAdaptation) - lastStateWeight
606 74894 : self%Chain%Weight(self%Stats%NumFunCall%accepted) = currentStateWeight ! needed for the restart mode, not needed in the fresh run
607 : end if
608 :
609 124416 : meanAccRateSinceStart = self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted) ! used only in fresh run, but not worth putting it in a conditional block.
610 : call self%Proposal%doAdaptation ( nd = nd & ! LCOV_EXCL_LINE
611 : , chainSize = self%Stats%NumFunCall%accepted - numFunCallAcceptedLastAdaptation + 1_IK & ! LCOV_EXCL_LINE
612 : , Chain = self%Chain%State(1:nd,numFunCallAcceptedLastAdaptation:self%Stats%NumFunCall%accepted) & ! LCOV_EXCL_LINE
613 : , ChainWeight = self%Chain%Weight(numFunCallAcceptedLastAdaptation:self%Stats%NumFunCall%accepted) & ! LCOV_EXCL_LINE
614 : , isFreshRun = self%isFreshRun & ! LCOV_EXCL_LINE
615 : , samplerUpdateIsGreedy = samplerUpdateIsGreedy & ! LCOV_EXCL_LINE
616 : , meanAccRateSinceStart = meanAccRateSinceStart & ! LCOV_EXCL_LINE
617 : , samplerUpdateSucceeded = samplerUpdateSucceeded & ! LCOV_EXCL_LINE
618 : , adaptationMeasure = AdaptationMeasure(dumint) & ! LCOV_EXCL_LINE
619 124416 : )
620 : #if defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED || ( (defined MATLAB_ENABLED || defined PYTHON_ENABLED || defined R_ENABLED) && !defined CAF_ENABLED && !defined MPI_ENABLED )
621 124416 : if(ProposalErr%occurred) then; self%Err%occurred = .true.; self%Err%msg = ProposalErr%msg; exit loopMarkovChain; return; end if
622 : #endif
623 124410 : if (self%isDryRun) SumAccRateSinceStart%acceptedRejected = meanAccRateSinceStart * real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
624 :
625 124410 : self%Chain%Weight(self%Stats%NumFunCall%accepted) = dumint ! needed for the restart mode, not needed in the fresh run
626 124410 : if (self%Stats%NumFunCall%accepted==numFunCallAcceptedLastAdaptation) then
627 : !adaptationMeasure = adaptationMeasure + adaptationMeasureDummy ! this is the worst-case upper-bound
628 49521 : self%Chain%Adaptation(self%Stats%NumFunCall%accepted) = min(1._RK, self%Chain%Adaptation(self%Stats%NumFunCall%accepted) + AdaptationMeasure(dumint)) ! this is the worst-case upper-bound
629 : else
630 : !adaptationMeasure = adaptationMeasureDummy
631 74889 : self%Chain%Adaptation(self%Stats%NumFunCall%accepted) = AdaptationMeasure(dumint)
632 74889 : self%Chain%Weight(numFunCallAcceptedLastAdaptation) = self%Chain%Weight(numFunCallAcceptedLastAdaptation) + lastStateWeight
633 : end if
634 124410 : if (samplerUpdateSucceeded) then
635 30423 : lastStateWeight = currentStateWeight ! self%Chain%Weight(self%Stats%NumFunCall%accepted) ! informative, do not remove
636 30423 : numFunCallAcceptedLastAdaptation = self%Stats%NumFunCall%accepted
637 : end if
638 :
639 124410 : counterAUP = 0_IK
640 124410 : counterAUC = counterAUC + 1_IK
641 : !if (counterAUC==self%SpecDRAM%AdaptiveUpdateCount%val) adaptationMeasure = 0._RK
642 :
643 : else blockSamplerAdaptation
644 :
645 159232 : AdaptationMeasure(self%Chain%Weight(self%Stats%NumFunCall%accepted)) = 0._RK
646 :
647 : end if blockSamplerAdaptation
648 :
649 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
650 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Proposal Adaptation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
651 :
652 : #if defined CAF_ENABLED || defined MPI_ENABLED
653 383579 : if (co_proposalFound_samplerUpdateOccurred(1)==1_IK) exit loopOverImages
654 :
655 : end do loopOverImages
656 : #endif
657 :
658 : #if defined CAF_ENABLED
659 :
660 194085 : if (self%SpecBase%ParallelizationModel%isSinglChain) then
661 95160 : call self%Timer%toc()
662 95160 : sync images(*)
663 95160 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
664 : end if
665 :
666 : else blockLeaderImage ! ATTN: This block should be executed only when singlChain parallelizationModel is requested
667 :
668 190320 : sync images(1)
669 :
670 : ! get the accepted proposal from the first image
671 :
672 190320 : call self%Timer%toc()
673 190320 : co_proposalFound_samplerUpdateOccurred(1:2) = co_proposalFound_samplerUpdateOccurred(1:2)[1]
674 190320 : if (co_proposalFound_samplerUpdateOccurred(1)==1_IK) co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,-1)[1]
675 190320 : if (co_proposalFound_samplerUpdateOccurred(2)==1_IK) call self%Proposal%bcastAdaptation()
676 190320 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
677 :
678 190320 : if (self%isDryRun) then
679 17018 : if (co_proposalFound_samplerUpdateOccurred(1)==1_IK) then
680 13972 : self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
681 13972 : numFunCallAcceptedPlusOne = self%Stats%NumFunCall%accepted + 1_IK
682 13972 : currentStateWeight = 1_IK
683 13972 : if (self%Stats%NumFunCall%accepted==self%Chain%Count%compact) then
684 0 : self%isFreshRun = .true.
685 0 : self%isDryRun = .false.
686 : end if
687 : else
688 3046 : currentStateWeight = currentStateWeight + self%Image%count
689 : end if
690 : #if defined DEBUG_ENABLED
691 173302 : elseif (co_proposalFound_samplerUpdateOccurred(1)==1_IK) then
692 121050 : self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
693 : #endif
694 : end if
695 :
696 : end if blockLeaderImage
697 :
698 : #elif defined MPI_ENABLED
699 :
700 : if (self%SpecBase%ParallelizationModel%isSinglChain .and. co_proposalFound_samplerUpdateOccurred(1)==0_IK) then
701 : ! LCOV_EXCL_START
702 : imageID = 0_IK ! broadcast rank # 0 to all processes, indicating unsuccessful sampling
703 : call mpi_bcast ( imageID & ! buffer
704 : , 1 & ! count
705 : , mpi_integer & ! datatype
706 : , 0 & ! root: broadcasting rank
707 : , mpi_comm_world & ! comm
708 : , ierrMPI & ! ierr
709 : )
710 : ! LCOV_EXCL_STOP
711 : end if
712 :
713 : else blockLeaderImage ! This block should be executed only when singlChain parallelizationModel is requested
714 :
715 : ! fetch the winning rank from the main process
716 :
717 : call self%Timer%toc()
718 : call mpi_bcast ( imageID & ! LCOV_EXCL_LINE ! buffer
719 : , 1 & ! LCOV_EXCL_LINE ! count
720 : , mpi_integer & ! LCOV_EXCL_LINE ! datatype
721 : , 0 & ! LCOV_EXCL_LINE ! root: broadcasting rank
722 : , mpi_comm_world & ! LCOV_EXCL_LINE ! comm
723 : , ierrMPI & ! LCOV_EXCL_LINE ! ierr
724 : )
725 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
726 :
727 : if (imageID>0_IK) then ! co_ProposalFound = .true., sampling successful
728 :
729 : ! broadcast co_LogFuncState from the winning image to all others
730 :
731 : call self%Timer%toc()
732 : call mpi_bcast ( co_LogFuncState(0:nd,-1) & ! LCOV_EXCL_LINE ! buffer
733 : , ndPlusOne & ! LCOV_EXCL_LINE ! count
734 : , mpi_double_precision & ! LCOV_EXCL_LINE ! datatype
735 : , imageID - 1 & ! LCOV_EXCL_LINE ! root: broadcasting rank
736 : , mpi_comm_world & ! LCOV_EXCL_LINE ! comm
737 : , ierrMPI & ! LCOV_EXCL_LINE ! ierr
738 : )
739 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
740 :
741 : end if
742 :
743 : if (self%isDryRun) then
744 : if (imageID>0_IK) then ! equivalent to co_proposalFound_samplerUpdateOccurred(1)==1_IK
745 : self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
746 : numFunCallAcceptedPlusOne = self%Stats%NumFunCall%accepted + 1_IK
747 : currentStateWeight = 1_IK
748 : if (self%Stats%NumFunCall%accepted==self%Chain%Count%compact) then
749 : self%isFreshRun = .true.
750 : self%isDryRun = .false.
751 : end if
752 : else
753 : currentStateWeight = currentStateWeight + self%Image%count
754 : end if
755 : #if defined DEBUG_ENABLED
756 : elseif (co_proposalFound_samplerUpdateOccurred(1)==1_IK) then
757 : self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
758 : #endif
759 : end if
760 :
761 : end if blockLeaderImage
762 :
763 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
764 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin Common block between all images %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
765 :
766 : ! broadcast samplerUpdateOccurred from the root process to all others, and broadcast proposal adaptations if needed
767 :
768 : if (self%SpecBase%ParallelizationModel%isSinglChain) then
769 : call self%Timer%toc()
770 : ! buffer: XXX: first element of co_proposalFound_samplerUpdateOccurred not needed except to end the simulation.
771 : ! This could perhaps be enhanced in the further to only pass one elemtn.
772 : call mpi_bcast ( co_proposalFound_samplerUpdateOccurred & ! buffer
773 : , 2 & ! count
774 : , mpi_integer & ! datatype
775 : , 0 & ! root: broadcasting rank
776 : , mpi_comm_world & ! comm
777 : , ierrMPI & ! ierr
778 : )
779 : if (co_proposalFound_samplerUpdateOccurred(2)==1_IK) call self%Proposal%bcastAdaptation()
780 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
781 : end if
782 :
783 : #endif
784 :
785 384405 : if (co_proposalFound_samplerUpdateOccurred(1) == -1_IK) exit loopMarkovChain ! we are done: co_missionAccomplished = .true.
786 :
787 383724 : co_AccRate(-1) = -1._RK ! counterDRS at which new proposal is accepted. This initialization is essential for all serial and parallel modes
788 383724 : maxLogFuncRejectedProposal = NEGINF_RK
789 :
790 : #if defined CAF_ENABLED || defined MPI_ENABLED
791 383724 : blockSingleChainParallelism: if (self%SpecBase%ParallelizationModel%isSinglChain) then
792 : #define SINGLCHAIN_PARALLELISM
793 : #include "ParaDXXX_mod@Kernel_smod@nextMove.inc.f90"
794 : #undef SINGLCHAIN_PARALLELISM
795 : else blockSingleChainParallelism
796 : #endif
797 : #include "ParaDXXX_mod@Kernel_smod@nextMove.inc.f90"
798 : #if defined CAF_ENABLED || defined MPI_ENABLED
799 : end if blockSingleChainParallelism
800 : #endif
801 :
802 384405 : cycle loopMarkovChain
803 :
804 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end Common block between all images %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
805 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
806 :
807 : end do loopMarkovChain
808 :
809 : #if (defined MPI_ENABLED || defined CAF_ENABLED) && (defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED)
810 687 : block; use Err_mod, only: bcastErr; call bcastErr(ProposalErr); end block
811 : #endif
812 687 : if (self%Err%occurred) return
813 :
814 681 : call self%Timer%toc()
815 :
816 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
817 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
818 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%* end of loopMarkovChain %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
819 :
820 681 : self%chain%Count%target = self%Stats%NumFunCall%accepted
821 681 : self%chain%Count%compact = self%Stats%NumFunCall%accepted
822 681 : self%chain%Count%verbose = self%Stats%NumFunCall%acceptedRejected
823 :
824 681 : if (noDelayedRejectionRequested) then
825 633 : self%Stats%NumFunCall%acceptedRejectedDelayed = self%Stats%NumFunCall%acceptedRejected
826 633 : SumAccRateSinceStart%acceptedRejectedDelayed = SumAccRateSinceStart%acceptedRejected
827 : end if
828 :
829 :
830 681 : if (self%Image%isLeader) then
831 : block
832 : integer(IK) :: i
833 263 : if (self%SpecDRAM%BurninAdaptationMeasure%val>0.999999_RK) then
834 261 : self%Stats%AdaptationBurninLoc%compact = 1_IK
835 261 : self%Stats%AdaptationBurninLoc%verbose = 1_IK
836 : else
837 2 : self%Stats%AdaptationBurninLoc%compact = self%Stats%NumFunCall%accepted
838 2 : self%Stats%AdaptationBurninLoc%verbose = self%Stats%NumFunCall%acceptedRejected - self%Chain%Weight(self%Stats%NumFunCall%accepted) + 1_IK
839 598 : loopAdaptationBurnin: do i = self%Stats%NumFunCall%accepted-1, 1, -1
840 597 : if (self%Chain%Adaptation(i)>self%SpecDRAM%BurninAdaptationMeasure%val) exit loopAdaptationBurnin
841 596 : self%Stats%AdaptationBurninLoc%compact = self%Stats%AdaptationBurninLoc%compact - 1_IK
842 598 : self%Stats%AdaptationBurninLoc%verbose = self%Stats%AdaptationBurninLoc%verbose - self%Chain%Weight(i)
843 : end do loopAdaptationBurnin
844 : end if
845 : end block
846 263 : self%Stats%BurninLoc%compact = self%Chain%BurninLoc(self%Stats%NumFunCall%accepted)
847 373 : self%Stats%BurninLoc%verbose = sum(self%Chain%Weight(1:self%Stats%BurninLoc%compact-1)) + 1_IK
848 697 : self%Stats%LogFuncMode%Crd = self%Chain%State(1:nd,self%Stats%LogFuncMode%Loc%compact)
849 1492 : self%Stats%LogFuncMode%Loc%verbose = sum(self%Chain%Weight(1:self%Stats%LogFuncMode%Loc%compact-1)) + 1_IK
850 263 : close(self%RestartFile%unit)
851 263 : close(self%ChainFile%unit)
852 : endif
853 :
854 681 : if (self%Image%isFirst) then
855 : ! LCOV_EXCL_START
856 : call write()
857 : !call execute_command_line(" ")
858 : flush(output_unit)
859 : ! LCOV_EXCL_STOP
860 : end if
861 :
862 : #if defined CAF_ENABLED || defined MPI_ENABLED
863 1368 : if (self%SpecBase%ParallelizationModel%isMultiChain) then
864 : #endif
865 54 : self%Stats%avgCommTimePerFunCall = 0._RK
866 54 : self%Stats%NumFunCall%acceptedRejectedDelayedUnused = self%Stats%NumFunCall%acceptedRejectedDelayed
867 54 : dumint = self%Stats%NumFunCall%acceptedRejectedDelayedUnused - acceptedRejectedDelayedUnusedRestartMode ! this is needed to avoid division-by-zero undefined behavior
868 54 : if (dumint/=0_IK) self%Stats%avgTimePerFunCalInSec = self%Stats%avgTimePerFunCalInSec / dumint
869 : #if defined CAF_ENABLED || defined MPI_ENABLED
870 627 : elseif(self%Image%isFirst) then
871 : ! LCOV_EXCL_START
872 : self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall / self%Stats%NumFunCall%acceptedRejectedDelayed
873 : dumint = self%Stats%NumFunCall%acceptedRejectedDelayedUnused - acceptedRejectedDelayedUnusedRestartMode ! this is needed to avoid division-by-zero undefined behavior
874 : if (dumint/=0_IK) self%Stats%avgTimePerFunCalInSec = (self%Stats%avgTimePerFunCalInSec / dumint) * self%Image%count
875 : ! LCOV_EXCL_STOP
876 : end if
877 : #endif
878 :
879 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
880 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
881 :
882 : contains
883 :
884 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
885 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
886 :
887 : !> \brief
888 : !> Write the latest simulated state to the output file.
889 76045 : subroutine writeOutput()
890 :
891 : implicit none
892 : integer(IK) :: j
893 :
894 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output write %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
895 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
896 :
897 : ! if new point has been sampled, write the previous sampled point to output file
898 :
899 76045 : blockOutputWrite: if (self%Stats%NumFunCall%accepted>0_IK) then
900 75808 : if (self%SpecBase%ChainFileFormat%isCompact) then
901 62176 : write(self%ChainFile%unit,self%ChainFile%format ) self%Chain%ProcessID(self%Stats%NumFunCall%accepted) &
902 62176 : , self%Chain%DelRejStage(self%Stats%NumFunCall%accepted) &
903 62176 : , self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted) &
904 62176 : , self%Chain%Adaptation(self%Stats%NumFunCall%accepted) &
905 62176 : , self%Chain%BurninLoc(self%Stats%NumFunCall%accepted) &
906 62176 : , self%Chain%Weight(self%Stats%NumFunCall%accepted) &
907 62176 : , self%Chain%LogFunc(self%Stats%NumFunCall%accepted) &
908 124352 : , self%Chain%State(1:nd,self%Stats%NumFunCall%accepted)
909 13632 : elseif (self%SpecBase%ChainFileFormat%isBinary) then
910 6816 : write(self%ChainFile%unit ) self%Chain%ProcessID(self%Stats%NumFunCall%accepted) &
911 6816 : , self%Chain%DelRejStage(self%Stats%NumFunCall%accepted) &
912 6816 : , self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted) &
913 6816 : , self%Chain%Adaptation(self%Stats%NumFunCall%accepted) &
914 6816 : , self%Chain%BurninLoc(self%Stats%NumFunCall%accepted) &
915 6816 : , self%Chain%Weight(self%Stats%NumFunCall%accepted) &
916 6816 : , self%Chain%LogFunc(self%Stats%NumFunCall%accepted) &
917 13632 : , self%Chain%State(1:nd,self%Stats%NumFunCall%accepted)
918 6816 : elseif (self%SpecBase%ChainFileFormat%isVerbose) then
919 30381 : do j = 1, self%Chain%Weight(self%Stats%NumFunCall%accepted)
920 23565 : write(self%ChainFile%unit,self%ChainFile%format ) self%Chain%ProcessID(self%Stats%NumFunCall%accepted) &
921 23565 : , self%Chain%DelRejStage(self%Stats%NumFunCall%accepted) &
922 23565 : , self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted) &
923 0 : , AdaptationMeasure(j) &
924 : !, self%Chain%Adaptation(self%Stats%NumFunCall%accepted) &
925 23565 : , self%Chain%BurninLoc(self%Stats%NumFunCall%accepted) &
926 23565 : , 1_IK &
927 23565 : , self%Chain%LogFunc(self%Stats%NumFunCall%accepted) &
928 53946 : , self%Chain%State(1:nd,self%Stats%NumFunCall%accepted)
929 : end do
930 : end if
931 : end if blockOutputWrite
932 :
933 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
934 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output write %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
935 :
936 693 : end subroutine writeOutput
937 :
938 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
939 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
940 :
941 : !> \brief
942 : !> Update the user with the simulation progress information.
943 : !>
944 : !> \remark
945 : !> Objects that change state in this subroutine are: `Timer`, `timeElapsedUntilLastReportInSeconds`, `sumAccRateLastReport`.
946 408 : subroutine reportProgress()
947 :
948 76045 : use Constants_mod, only: CARRIAGE_RETURN
949 : implicit none
950 :
951 408 : real(RK) :: meanAccRateSinceStart
952 408 : real(RK) :: meanAccRateSinceLastReport
953 408 : real(RK) :: estimatedTimeToFinishInSeconds
954 408 : real(RK) :: timeElapsedSinceLastReportInSeconds
955 : integer(IK) :: numFunCallAccepted_dummy ! used in restart mode
956 :
957 :
958 408 : if (self%isFreshRun) then
959 :
960 363 : call self%Timer%toc()
961 363 : timeElapsedSinceLastReportInSeconds = self%Timer%Time%total - timeElapsedUntilLastReportInSeconds
962 363 : timeElapsedUntilLastReportInSeconds = self%Timer%Time%total
963 363 : meanAccRateSinceStart = SumAccRateSinceStart%acceptedRejected / real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
964 363 : meanAccRateSinceLastReport = (SumAccRateSinceStart%acceptedRejected-sumAccRateLastReport) * inverseProgressReportPeriod
965 363 : estimatedTimeToFinishInSeconds = getRemainingSimulationFraction() * self%Timer%Time%total
966 :
967 363 : write( self%TimeFile%unit, self%TimeFile%format ) self%Stats%NumFunCall%acceptedRejected &
968 363 : , self%Stats%NumFunCall%accepted &
969 : , meanAccRateSinceStart &
970 363 : , meanAccRateSinceLastReport &
971 363 : , timeElapsedSinceLastReportInSeconds &
972 363 : , self%Timer%Time%total & ! timeElapsedUntilLastReportInSeconds
973 726 : , estimatedTimeToFinishInSeconds
974 363 : flush(self%TimeFile%unit)
975 :
976 : else
977 :
978 810 : block
979 : use String_mod, only: String_type
980 45 : type(String_type) :: Record
981 45 : allocate( character(600) :: Record%value )
982 45 : read(self%TimeFile%unit, "(A)" ) Record%value
983 1035 : Record%Parts = Record%split(trim(adjustl(Record%value)),self%SpecBase%OutputDelimiter%val,Record%nPart)
984 45 : read(Record%Parts(1)%record,*) numFunCallAcceptedRejectedLastReport
985 45 : read(Record%Parts(2)%record,*) numFunCallAccepted_dummy
986 45 : read(Record%Parts(3)%record,*) meanAccRateSinceStart
987 45 : read(Record%Parts(4)%record,*) meanAccRateSinceLastReport
988 45 : read(Record%Parts(5)%record,*) timeElapsedSinceLastReportInSeconds
989 45 : read(Record%Parts(6)%record,*) timeElapsedUntilLastReportInSeconds
990 855 : read(Record%Parts(7)%record,*) estimatedTimeToFinishInSeconds
991 : end block
992 45 : SumAccRateSinceStart%acceptedRejected = meanAccRateSinceStart * numFunCallAcceptedRejectedLastReport
993 :
994 : end if
995 :
996 : ! report progress in the standard output
997 408 : if (self%Image%isFirst) then
998 : ! LCOV_EXCL_START
999 : write( &
1000 : #if defined MEXPRINT_ENABLED
1001 : txt, &
1002 : #else
1003 : output_unit, advance = "no", &
1004 : #endif
1005 : fmt = "(A,25X,*(A25,3X))" ) &
1006 : CARRIAGE_RETURN, &
1007 : num2str(self%Stats%NumFunCall%accepted)//" / "//num2str(self%Stats%NumFunCall%acceptedRejected,formatIn="(1I10)"), &
1008 : num2str(meanAccRateSinceLastReport,"(1F11.4)")//" / "//num2str(SumAccRateSinceStart%acceptedRejected / real(self%Stats%NumFunCall%acceptedRejected,kind=RK),"(1F11.4)"), &
1009 : num2str(timeElapsedUntilLastReportInSeconds,"(1F11.4)")//" / "//num2str(estimatedTimeToFinishInSeconds,"(1F11.4)")
1010 : #if defined MEXPRINT_ENABLED
1011 : call write(string=txt,advance=.false.)
1012 : #else
1013 : !call execute_command_line(" ")
1014 : flush(output_unit)
1015 : #endif
1016 : ! LCOV_EXCL_STOP
1017 : end if
1018 :
1019 408 : numFunCallAcceptedRejectedLastReport = self%Stats%NumFunCall%acceptedRejected
1020 408 : sumAccRateLastReport = SumAccRateSinceStart%acceptedRejected
1021 :
1022 408 : end subroutine reportProgress
1023 :
1024 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1025 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1026 :
1027 : !> \brief
1028 : !> Return the predicted remaining fraction of the simulation for the purpose of updating the user.
1029 : !>
1030 : !> \return
1031 : !> `remainingSimulationFraction` : The predicted remaining fraction of the simulation.
1032 363 : pure function getRemainingSimulationFraction() result(remainingSimulationFraction)
1033 408 : use Constants_mod, only: IK, RK
1034 : implicit none
1035 : real(RK) :: remainingSimulationFraction
1036 363 : remainingSimulationFraction = real(self%SpecMCMC%ChainSize%val-self%Stats%NumFunCall%accepted,kind=RK) / self%Stats%NumFunCall%accepted
1037 363 : end function getRemainingSimulationFraction
1038 :
1039 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1040 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1041 :
1042 : end subroutine runKernel
1043 :
1044 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1045 :
1046 : !> \brief
1047 : !> Return the predicted burnin location within the currently sampled chain.
1048 : !>
1049 : !> @param[in] lenLogFunc : The length of the input `lenLogFunc`.
1050 : !> @param[in] refLogFunc : The reference logFunc value with respect to which the relative importance of the points is gauged.
1051 : !> @param[in] LogFunc : The vector of length `lenLogFunc` of log(function) values.
1052 : !>
1053 : !> \return
1054 : !> `burninLoc` : The location of burnin in the input `LogFunc` vector.
1055 75750 : pure function getBurninLoc(lenLogFunc,refLogFunc,LogFunc) result(burninLoc)
1056 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1057 : !DEC$ ATTRIBUTES DLLEXPORT :: getBurninLoc
1058 : #endif
1059 : use Constants_mod, only: IK, RK
1060 : implicit none
1061 : integer(IK), intent(in) :: lenLogFunc
1062 : real(RK) , intent(in) :: refLogFunc,LogFunc(lenLogFunc)
1063 75750 : real(RK) :: negLogIncidenceProb
1064 : integer(IK) :: burninLoc
1065 75750 : negLogIncidenceProb = log( real(lenLogFunc,kind=RK) )
1066 75750 : burninLoc = 0_IK
1067 : do ! LCOV_EXCL_LINE
1068 98672 : burninLoc = burninLoc + 1_IK
1069 98672 : if ( burninLoc<lenLogFunc .and. refLogFunc-LogFunc(burninLoc)>negLogIncidenceProb ) cycle
1070 : !burninLoc = burninLoc - 1
1071 98672 : exit
1072 : end do
1073 75750 : end function getBurninLoc
1074 :
1075 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1076 :
1077 : #undef ParaDXXX_type
1078 : #undef ParaDXXXProposalAbstract_mod
|
