The ParaMonte Documentation Website
Current view: top level - kernel - ParaDXXX_mod@Kernel_smod.inc.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 329 333 98.8 %
Date: 2021-01-08 13:07:16 Functions: 7 7 100.0 %
Legend: Lines: hit not hit

          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             :         allocate(co_LogFuncState(0:nd,-1:self%SpecDRAM%DelayedRejectionCount%val)[*])
     184             :         allocate(co_AccRate(-1:self%SpecDRAM%DelayedRejectionCount%val)[*])                         ! the negative element will contain counterDRS
     185             : #else
     186         687 :         allocate(co_LogFuncState(0:nd,-1:self%SpecDRAM%DelayedRejectionCount%val))
     187         687 :         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         687 :         if (allocated(AccRateMatrix)) deallocate(AccRateMatrix)
     195         687 :         allocate(AccRateMatrix(-1:self%SpecDRAM%DelayedRejectionCount%val,1:self%Image%count))      ! the negative element will contain counterDRS
     196        7356 :         AccRateMatrix = 0._RK                                                                       ! -huge(1._RK)  ! debug
     197        2748 :         AccRateMatrix(0,1:self%Image%count) = 1._RK                                                 ! initial acceptance rate for the first zeroth DR stage.
     198         687 :         ndPlusOne = nd + 1_IK
     199         687 :         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      378659 :         loopMarkovChain: do
     392             : 
     393      379346 :             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      379346 :             blockLeaderImage: if (self%Image%isLeader) then
     397             : #endif
     398             : 
     399      191332 :                 co_proposalFound_samplerUpdateOccurred(1) = 0_IK ! co_proposalFound = .false.
     400      191332 :                 samplerUpdateIsGreedy = counterAUC < self%SpecDRAM%GreedyAdaptationCount%val
     401             : 
     402             : #if defined CAF_ENABLED || defined MPI_ENABLED
     403      376056 :                 loopOverImages: do imageID = imageStartID, imageEndID
     404             : #if defined CAF_ENABLED
     405             :                     call self%Timer%toc()
     406             :                     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             :                     call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     408             : #elif defined MPI_ENABLED
     409      491234 :                     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      278878 :                     counterDRS = nint(co_AccRate(-1),kind=IK)
     414      278878 :                     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      278878 :                     blockProposalAccepted: if ( co_proposalFound_samplerUpdateOccurred(1) == 1_IK ) then ! co_proposalAccepted = .true.
     424             : 
     425       94149 :                         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       94149 :                         if (self%SpecBase%ParallelizationModel%isSinglChain) then
     431       67720 :                             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       67720 :                                             )
     440             :                             ! broadcast co_LogFuncState from the winning image to all others
     441           0 :                             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       67720 :                                             )
     448       67720 :                             call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     449             :                         end if
     450             : #elif defined CAF_ENABLED
     451             :                         if (imageID/=self%Image%id) then ! Avoid remote connection for something that is available locally.
     452             :                             call self%Timer%toc()
     453             :                             co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,-1)[imageID]
     454             :                             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       94149 :                         if (counterAUP==0_IK .and. samplerUpdateSucceeded) then
     460       14670 :                             numFunCallAcceptedLastAdaptation = numFunCallAcceptedLastAdaptation + 1_IK
     461       14670 :                             lastStateWeight = 0_IK
     462             :                         end if
     463             : 
     464       94149 :                         blockFreshDryRun: if (self%isFreshRun) then
     465             : 
     466       75749 :                             call writeOutput()
     467             : 
     468       75749 :                             self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
     469             : 
     470       75749 :                             self%Chain%ProcessID(self%Stats%NumFunCall%accepted)    = imageID
     471       75749 :                             self%Chain%DelRejStage(self%Stats%NumFunCall%accepted)  = counterDRS
     472       75749 :                             self%Chain%Adaptation(self%Stats%NumFunCall%accepted)   = 0._RK ! adaptationMeasure
     473       75749 :                             self%Chain%Weight(self%Stats%NumFunCall%accepted)       = 0_IK
     474       75749 :                             self%Chain%LogFunc(self%Stats%NumFunCall%accepted)      = co_LogFuncState(0,-1)
     475      196816 :                             self%Chain%State(1:nd,self%Stats%NumFunCall%accepted)   = co_LogFuncState(1:nd,-1)
     476             : 
     477             :                             ! find the burnin point
     478             : 
     479      151498 :                             self%Chain%BurninLoc(self%Stats%NumFunCall%accepted) = getBurninLoc ( lenLogFunc    = self%Stats%NumFunCall%accepted                        &
     480             :                                                                                                 , refLogFunc    = self%Stats%LogFuncMode%val                            &
     481      302996 :                                                                                                 , LogFunc       = self%Chain%LogFunc(1:self%Stats%NumFunCall%accepted)  &
     482       75749 :                                                                                                 )
     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       94149 :                         if (self%Stats%LogFuncMode%val < self%Chain%LogFunc(self%Stats%NumFunCall%accepted)) then
     501         464 :                             self%Stats%LogFuncMode%val = self%Chain%LogFunc(self%Stats%NumFunCall%accepted)
     502         464 :                             self%Stats%LogFuncMode%Loc%compact = self%Stats%NumFunCall%accepted
     503             :                         end if
     504             : 
     505       94149 :                         SumAccRateSinceStart%acceptedRejected = SumAccRateSinceStart%acceptedRejected + co_AccRate(counterDRS)
     506             : 
     507             :                     else blockProposalAccepted
     508             : 
     509      184729 :                         counterDRS = self%SpecDRAM%DelayedRejectionCount%val
     510      184729 :                         SumAccRateSinceStart%acceptedRejected = SumAccRateSinceStart%acceptedRejected + co_AccRate(counterDRS)
     511             : 
     512             :                     end if blockProposalAccepted
     513             : 
     514             :                     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     515             :                     !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% blockProposalAccepted %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     516             : 
     517      278878 :                     counterAUP = counterAUP + 1_IK
     518      278878 :                     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      278878 :                     if (counterPRP == self%SpecBase%ProgressReportPeriod%val) then
     525         144 :                         counterPRP = 0_IK
     526         144 :                         call reportProgress()
     527             :                     end if
     528             : 
     529      278878 :                     currentStateWeight = currentStateWeight + 1_IK
     530      278878 :                     self%Stats%NumFunCall%acceptedRejected = self%Stats%NumFunCall%acceptedRejected + 1_IK
     531             : 
     532      278878 :                     if (delayedRejectionRequested) then
     533       87246 :                         SumAccRateSinceStart%acceptedRejectedDelayed = SumAccRateSinceStart%acceptedRejectedDelayed + sum(co_AccRate(0:counterDRS))
     534       18750 :                         self%Stats%NumFunCall%acceptedRejectedDelayed = self%Stats%NumFunCall%acceptedRejectedDelayed + counterDRS + 1_IK
     535             :                     end if
     536             : 
     537      278878 :                     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      218593 :                         self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted)  = SumAccRateSinceStart%acceptedRejected / real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
     539      218593 :                         self%Chain%Weight(self%Stats%NumFunCall%accepted)       = self%Chain%Weight(self%Stats%NumFunCall%accepted) + 1_IK
     540      218593 :                         if (adaptationMeasureLen<self%Chain%Weight(self%Stats%NumFunCall%accepted)) then
     541          10 :                             allocate(adaptationMeasurePlaceHolder(2*adaptationMeasureLen))
     542        6910 :                             adaptationMeasurePlaceHolder(1:adaptationMeasureLen) = AdaptationMeasure
     543          10 :                             call move_alloc(from=adaptationMeasurePlaceHolder,to=AdaptationMeasure)
     544          10 :                             adaptationMeasureLen = 2_IK * adaptationMeasureLen
     545             :                         end if
     546             :                     else
     547             : #if defined CAF_ENABLED || defined MPI_ENABLED
     548       60285 :                         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      278878 :                     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      278615 :                     blockSamplerAdaptation: if ( counterAUC < self%SpecDRAM%AdaptiveUpdateCount%val .and. counterAUP == self%SpecDRAM%AdaptiveUpdatePeriod%val ) then
     587             : 
     588      122107 :                         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      122107 :                         dumint = self%Chain%Weight(self%Stats%NumFunCall%accepted) ! needed for the restart mode, not needed in the fresh run
     593      122107 :                         if (self%Stats%NumFunCall%accepted==numFunCallAcceptedLastAdaptation) then    ! no new point has been accepted since last time
     594       46669 :                             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       75438 :                             self%Chain%Weight(numFunCallAcceptedLastAdaptation) = self%Chain%Weight(numFunCallAcceptedLastAdaptation) - lastStateWeight
     606       75438 :                             self%Chain%Weight(self%Stats%NumFunCall%accepted) = currentStateWeight ! needed for the restart mode, not needed in the fresh run
     607             :                         end if
     608             : 
     609      122107 :                         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      122107 :                                                         )
     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      122107 :                         if(ProposalErr%occurred) then; self%Err%occurred = .true.; self%Err%msg = ProposalErr%msg; exit loopMarkovChain; return; end if
     622             : #endif
     623      122101 :                         if (self%isDryRun) SumAccRateSinceStart%acceptedRejected = meanAccRateSinceStart * real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
     624             : 
     625      122101 :                         self%Chain%Weight(self%Stats%NumFunCall%accepted) = dumint   ! needed for the restart mode, not needed in the fresh run
     626      122101 :                         if (self%Stats%NumFunCall%accepted==numFunCallAcceptedLastAdaptation) then
     627             :                             !adaptationMeasure = adaptationMeasure + adaptationMeasureDummy ! this is the worst-case upper-bound
     628       46667 :                             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       75434 :                             self%Chain%Adaptation(self%Stats%NumFunCall%accepted) = AdaptationMeasure(dumint)
     632       75434 :                             self%Chain%Weight(numFunCallAcceptedLastAdaptation) = self%Chain%Weight(numFunCallAcceptedLastAdaptation) + lastStateWeight
     633             :                         end if
     634      122101 :                         if (samplerUpdateSucceeded) then
     635       30425 :                             lastStateWeight = currentStateWeight ! self%Chain%Weight(self%Stats%NumFunCall%accepted) ! informative, do not remove
     636       30425 :                             numFunCallAcceptedLastAdaptation = self%Stats%NumFunCall%accepted
     637             :                         end if
     638             : 
     639      122101 :                         counterAUP = 0_IK
     640      122101 :                         counterAUC = counterAUC + 1_IK
     641             :                         !if (counterAUC==self%SpecDRAM%AdaptiveUpdateCount%val) adaptationMeasure = 0._RK
     642             : 
     643             :                     else blockSamplerAdaptation
     644             : 
     645      156508 :                         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      375787 :                     if (co_proposalFound_samplerUpdateOccurred(1)==1_IK) exit loopOverImages
     654             : 
     655             :                 end do loopOverImages
     656             : #endif
     657             : 
     658             : #if defined CAF_ENABLED
     659             : 
     660             :                 if (self%SpecBase%ParallelizationModel%isSinglChain) then
     661             :                     call self%Timer%toc()
     662             :                     sync images(*)
     663             :                     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             :                 sync images(1)
     669             : 
     670             :                 ! get the accepted proposal from the first image
     671             : 
     672             :                 call self%Timer%toc()
     673             :                 co_proposalFound_samplerUpdateOccurred(1:2) = co_proposalFound_samplerUpdateOccurred(1:2)[1]
     674             :                 if (co_proposalFound_samplerUpdateOccurred(1)==1_IK) co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,-1)[1]
     675             :                 if (co_proposalFound_samplerUpdateOccurred(2)==1_IK) call self%Proposal%bcastAdaptation()
     676             :                 call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     677             : 
     678             :                 if (self%isDryRun) then
     679             :                     if (co_proposalFound_samplerUpdateOccurred(1)==1_IK) then
     680             :                         self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
     681             :                         numFunCallAcceptedPlusOne = self%Stats%NumFunCall%accepted + 1_IK
     682             :                         currentStateWeight = 1_IK
     683             :                         if (self%Stats%NumFunCall%accepted==self%Chain%Count%compact) then
     684             :                             self%isFreshRun = .true.
     685             :                             self%isDryRun = .false.
     686             :                         end if
     687             :                     else
     688             :                         currentStateWeight = currentStateWeight + self%Image%count
     689             :                     end if
     690             : #if defined DEBUG_ENABLED
     691             :                 elseif (co_proposalFound_samplerUpdateOccurred(1)==1_IK) then
     692             :                     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      191326 :                 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      188014 :                 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      188014 :                                 )
     725      188014 :                 call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     726             : 
     727      188014 :                 if (imageID>0_IK) then ! co_ProposalFound = .true., sampling successful
     728             : 
     729             :                     ! broadcast co_LogFuncState from the winning image to all others
     730             : 
     731      135440 :                     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      135440 :                                     )
     739      135440 :                     call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     740             : 
     741             :                 end if
     742             : 
     743      188014 :                 if (self%isDryRun) then
     744       17070 :                     if (imageID>0_IK) then ! equivalent to co_proposalFound_samplerUpdateOccurred(1)==1_IK
     745       14000 :                         self%Stats%NumFunCall%accepted = self%Stats%NumFunCall%accepted + 1_IK
     746       14000 :                         numFunCallAcceptedPlusOne = self%Stats%NumFunCall%accepted + 1_IK
     747       14000 :                         currentStateWeight = 1_IK
     748       14000 :                         if (self%Stats%NumFunCall%accepted==self%Chain%Count%compact) then
     749          28 :                             self%isFreshRun = .true.
     750          28 :                             self%isDryRun = .false.
     751             :                         end if
     752             :                     else
     753        3070 :                         currentStateWeight = currentStateWeight + self%Image%count
     754             :                     end if
     755             : #if defined DEBUG_ENABLED
     756      170944 :                 elseif (co_proposalFound_samplerUpdateOccurred(1)==1_IK) then
     757      121050 :                     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      379340 :             if (self%SpecBase%ParallelizationModel%isSinglChain) then
     769      282021 :                 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      282021 :                                 )
     779      282021 :                 if (co_proposalFound_samplerUpdateOccurred(2)==1_IK) call self%Proposal%bcastAdaptation()
     780      282021 :                 call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
     781             :             end if
     782             : 
     783             : #endif
     784             : 
     785      379340 :             if (co_proposalFound_samplerUpdateOccurred(1) == -1_IK) exit loopMarkovChain   ! we are done: co_missionAccomplished = .true.
     786             : 
     787      378659 :             co_AccRate(-1) = -1._RK ! counterDRS at which new proposal is accepted. This initialization is essential for all serial and parallel modes
     788      378659 :             maxLogFuncRejectedProposal = NEGINF_RK
     789             : 
     790             : #if defined CAF_ENABLED || defined MPI_ENABLED
     791      378659 :             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      378659 :             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         597 :                         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         383 :             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        1600 :             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       76044 :         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       76044 :             blockOutputWrite: if (self%Stats%NumFunCall%accepted>0_IK) then
     900       75807 :                 if (self%SpecBase%ChainFileFormat%isCompact) then
     901       62175 :                     write(self%ChainFile%unit,self%ChainFile%format ) self%Chain%ProcessID(self%Stats%NumFunCall%accepted)      &
     902       62175 :                                                                     , self%Chain%DelRejStage(self%Stats%NumFunCall%accepted)    &
     903       62175 :                                                                     , self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted)    &
     904       62175 :                                                                     , self%Chain%Adaptation(self%Stats%NumFunCall%accepted)     &
     905       62175 :                                                                     , self%Chain%BurninLoc(self%Stats%NumFunCall%accepted)      &
     906       62175 :                                                                     , self%Chain%Weight(self%Stats%NumFunCall%accepted)         &
     907       62175 :                                                                     , self%Chain%LogFunc(self%Stats%NumFunCall%accepted)        &
     908      124350 :                                                                     , 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       30037 :                     do j = 1, self%Chain%Weight(self%Stats%NumFunCall%accepted)
     920       23221 :                     write(self%ChainFile%unit,self%ChainFile%format ) self%Chain%ProcessID(self%Stats%NumFunCall%accepted)      &
     921       23221 :                                                                     , self%Chain%DelRejStage(self%Stats%NumFunCall%accepted)    &
     922       23221 :                                                                     , self%Chain%MeanAccRate(self%Stats%NumFunCall%accepted)    &
     923           0 :                                                                     , AdaptationMeasure(j)                                      &
     924             :                                                                    !, self%Chain%Adaptation(self%Stats%NumFunCall%accepted)     &
     925       23221 :                                                                     , self%Chain%BurninLoc(self%Stats%NumFunCall%accepted)      &
     926       23221 :                                                                     , 1_IK                                                      &
     927       23221 :                                                                     , self%Chain%LogFunc(self%Stats%NumFunCall%accepted)        &
     928       53258 :                                                                     , 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         687 :         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         407 :         subroutine reportProgress()
     947             : 
     948       76044 :             use Constants_mod, only: CARRIAGE_RETURN
     949             :             implicit none
     950             : 
     951         407 :             real(RK)                            :: meanAccRateSinceStart
     952         407 :             real(RK)                            :: meanAccRateSinceLastReport
     953         407 :             real(RK)                            :: estimatedTimeToFinishInSeconds
     954         407 :             real(RK)                            :: timeElapsedSinceLastReportInSeconds
     955             :             integer(IK)                         :: numFunCallAccepted_dummy ! used in restart mode
     956             : 
     957             : 
     958         407 :             if (self%isFreshRun) then
     959             : 
     960         361 :                 call self%Timer%toc()
     961         361 :                 timeElapsedSinceLastReportInSeconds = self%Timer%Time%total - timeElapsedUntilLastReportInSeconds
     962         361 :                 timeElapsedUntilLastReportInSeconds = self%Timer%Time%total
     963         361 :                 meanAccRateSinceStart = SumAccRateSinceStart%acceptedRejected / real(self%Stats%NumFunCall%acceptedRejected,kind=RK)
     964         361 :                 meanAccRateSinceLastReport = (SumAccRateSinceStart%acceptedRejected-sumAccRateLastReport) * inverseProgressReportPeriod
     965         361 :                 estimatedTimeToFinishInSeconds = getRemainingSimulationFraction() * self%Timer%Time%total
     966             : 
     967         361 :                 write( self%TimeFile%unit, self%TimeFile%format ) self%Stats%NumFunCall%acceptedRejected    &
     968         361 :                                                                 , self%Stats%NumFunCall%accepted            &
     969             :                                                                 , meanAccRateSinceStart                     &
     970         361 :                                                                 , meanAccRateSinceLastReport                &
     971         361 :                                                                 , timeElapsedSinceLastReportInSeconds       &
     972         361 :                                                                 , self%Timer%Time%total                     & ! timeElapsedUntilLastReportInSeconds
     973         722 :                                                                 , estimatedTimeToFinishInSeconds
     974         361 :                 flush(self%TimeFile%unit)
     975             : 
     976             :             else
     977             : 
     978         828 :                 block
     979             :                     use String_mod, only: String_type
     980          46 :                     type(String_type) :: Record
     981          46 :                     allocate( character(600) :: Record%value )
     982          46 :                     read(self%TimeFile%unit, "(A)" ) Record%value
     983         736 :                     Record%Parts = Record%split(trim(adjustl(Record%value)),self%SpecBase%OutputDelimiter%val,Record%nPart)
     984          46 :                     read(Record%Parts(1)%record,*) numFunCallAcceptedRejectedLastReport
     985          46 :                     read(Record%Parts(2)%record,*) numFunCallAccepted_dummy
     986          46 :                     read(Record%Parts(3)%record,*) meanAccRateSinceStart
     987          46 :                     read(Record%Parts(4)%record,*) meanAccRateSinceLastReport
     988          46 :                     read(Record%Parts(5)%record,*) timeElapsedSinceLastReportInSeconds
     989          46 :                     read(Record%Parts(6)%record,*) timeElapsedUntilLastReportInSeconds
     990         874 :                     read(Record%Parts(7)%record,*) estimatedTimeToFinishInSeconds
     991             :                 end block
     992          46 :                 SumAccRateSinceStart%acceptedRejected = meanAccRateSinceStart * numFunCallAcceptedRejectedLastReport
     993             : 
     994             :             end if
     995             : 
     996             :             ! report progress in the standard output
     997         407 :             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         407 :             numFunCallAcceptedRejectedLastReport = self%Stats%NumFunCall%acceptedRejected
    1020         407 :             sumAccRateLastReport = SumAccRateSinceStart%acceptedRejected
    1021             : 
    1022         407 :         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         361 :         pure function getRemainingSimulationFraction() result(remainingSimulationFraction)
    1033         407 :             use Constants_mod, only: IK, RK
    1034             :             implicit none
    1035             :             real(RK) :: remainingSimulationFraction
    1036         361 :             remainingSimulationFraction = real(self%SpecMCMC%ChainSize%val-self%Stats%NumFunCall%accepted,kind=RK) / self%Stats%NumFunCall%accepted
    1037         361 :         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       75749 :     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       75749 :         real(RK)                :: negLogIncidenceProb
    1064             :         integer(IK)             :: burninLoc
    1065       75749 :         negLogIncidenceProb = log( real(lenLogFunc,kind=RK) )
    1066       75749 :         burninLoc = 0_IK
    1067             :         do ! LCOV_EXCL_LINE
    1068       98478 :             burninLoc = burninLoc + 1_IK
    1069       98478 :             if ( burninLoc<lenLogFunc .and. refLogFunc-LogFunc(burninLoc)>negLogIncidenceProb ) cycle
    1070             :             !burninLoc = burninLoc - 1
    1071       75749 :             exit
    1072             :         end do
    1073       75749 :     end function getBurninLoc
    1074             : 
    1075             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1076             : 
    1077             : #undef ParaDXXX_type
    1078             : #undef ParaDXXXProposalAbstract_mod

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021