Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 1625463 : loopNextMove: do counterDRS = 0, spec%proposalDelayedRejectionCount%val
18 : #if OMP_ENABLED || CAFMPI_SINGLCHAIN_ENABLED
19 : ! stat%numFunCallAcceptedRejectedDelayedUnused is relevant only on the first image, despite being updated by all images
20 : stat%numFunCallAcceptedRejectedDelayedUnused = stat%numFunCallAcceptedRejectedDelayedUnused + spec%image%count
21 : #endif
22 : #if OMP_ENABLED
23 : do imageID = 1, spec%image%count
24 : call setProposalStateNew(spec, proposal, counterDRS, co_logFuncState(1 : ndim, imageID, counterDRS - 1), co_logFuncState(1 : ndim, imageID, counterDRS), spec%rng, err)
25 : if(err%occurred) then; err%msg = getFine(__FILE__, __LINE__)//PROCEDURE_NAME//SK_": "//err%msg; return; end if
26 : end do
27 : #else
28 1079152 : call setProposalStateNew(spec, proposal, counterDRS, co_logFuncState(1 : ndim, counterDRS - 1), co_logFuncState(1 : ndim, counterDRS), spec%rng, err)
29 1079152 : if(err%occurred) then; err%msg = getFine(__FILE__, __LINE__)//PROCEDURE_NAME//SK_": "//err%msg; return; end if
30 : #endif
31 : !!#if (MATLAB_ENABLED || PYTHON_ENABLED || R_ENABLED) && !CAF_ENABLED && !MPI_ENABLED
32 : !! if(err%occurred) then; err%msg = PROCEDURE_NAME//SK_": "//err%msg; return; end if
33 : !!#elif CODECOV_ENABLED || SAMPLER_TEST_ENABLED
34 : !#if CODECOV_ENABLED || SAMPLER_TEST_ENABLED
35 : ! ! This block is exclusively used to test the deterministic restart functionality of ParaXXXX samplers.
36 : ! ! This block must not be activated under any other circumstances.
37 : ! ! This block must be executed by all images.
38 : ! !spec%testSamplingCounter = spec%testSamplingCounter + 1_IK
39 : ! !err%occurred = err%occurred .or. spec%testSamplingCountTarget < spec%testSamplingCounter
40 : ! !if (err%occurred) then
41 : ! ! if (spec%testSamplingCounter > spec%testSamplingCountTarget) err%msg = "The simulation was interrupted at the requested sampling count for restart testing purposes."
42 : ! ! SET_CAFMPI(call isFailedImage(err%occurred))
43 : ! ! return
44 : ! !end if
45 : !#endif
46 : ! The following random function call is only needed in fresh runs to evaluate the acceptance of a proposal.
47 : ! However, it is taken out of the subsequent loops to achieve 100% deterministic reproducibility when
48 : ! a simulation is restarted.
49 1079152 : call setUnifRand(spec%rng, unifrnd)
50 1625450 : if (spec%run%is%new .or. numFunCallAcceptedPlusOne == stat%cfc%nsam) then
51 : #if OMP_ENABLED && (MATLAB_ENABLED || PYTHON_ENABLED || R_ENABLED)
52 : block
53 : real(RKD) :: avgTimePerFunCall, avgCommPerFunCall
54 : !avgTimePerFunCall = epsilon(0._RKD); avgCommPerFunCall = epsilon(0._RKD)
55 : !block
56 : !use pm_io
57 : !call disp%skip(count = 2)
58 : !call disp%show("before")
59 : !call disp%show("co_logFuncState(:, 1 : spec%image%count, counterDRS)")
60 : !call disp%show( co_logFuncState(:, 1 : spec%image%count, counterDRS) )
61 : !call disp%skip(count = 2)
62 : !end block
63 : mold = getLogFunc(co_logFuncState(:, 1 : spec%image%count, counterDRS), avgTimePerFunCall, avgCommPerFunCall)
64 : !block
65 : !use pm_io
66 : !call disp%skip(count = 2)
67 : !call disp%show("after")
68 : !call disp%show("co_logFuncState(:, 1 : spec%image%count, counterDRS)")
69 : !call disp%show( co_logFuncState(:, 1 : spec%image%count, counterDRS) )
70 : !call disp%skip(count = 2)
71 : !end block
72 : stat%avgTimePerFunCall = stat%avgTimePerFunCall + avgTimePerFunCall
73 : stat%avgCommPerFunCall = stat%avgCommPerFunCall + avgCommPerFunCall
74 : end block
75 : #elif OMP_ENABLED
76 : block
77 : use pm_timer, only: getTime => getTimeOMP
78 : real(RKD) :: start, avgTimePerFunCall
79 : avgTimePerFunCall = 0._RKD
80 : stat%timer%clock = stat%timer%time()
81 : ! preliminary tests indicate that false sharing is not a concern here.
82 : ! nevertheless, the current implementation attempts to avoid it by assuming a 128-byte cacheline.
83 : !!$omp parallel do reduction(+ : avgTimePerFunCall) num_threads(spec%image%count) default(none) shared(co_logFuncState, counterDRS, ndim, spec) private(imageID, start, logFuncState)
84 : !$omp parallel do reduction(+ : avgTimePerFunCall) num_threads(spec%image%count) default(none) shared(co_logFuncState, counterDRS, ndim, spec, logFuncState) private(imageID, start)
85 : do imageID = 1, spec%image%count
86 : start = getTime()
87 : logFuncState(1, imageID) = getLogFunc(co_logFuncState(1 : ndim, imageID, counterDRS))
88 : !avgTimePerFunCall(imageID * JUMP) = getTime() - start
89 : !co_logFuncState(0, imageID, counterDRS) = getLogFunc(co_logFuncState(1 : ndim, imageID, counterDRS))
90 : avgTimePerFunCall = avgTimePerFunCall + (getTime() - start)
91 : end do
92 : !$omp end parallel do
93 : co_logFuncState(0, :, counterDRS) = logFuncState(1, :)
94 : avgTimePerFunCall = avgTimePerFunCall / spec%image%count
95 : stat%avgTimePerFunCall = stat%avgTimePerFunCall + avgTimePerFunCall
96 : stat%avgCommPerFunCall = stat%avgCommPerFunCall + stat%timer%time() - stat%timer%clock - avgTimePerFunCall
97 : end block
98 : #else
99 1079152 : stat%timer%clock = stat%timer%time()
100 1079152 : co_logFuncState(0, counterDRS) = getLogFunc(co_logFuncState(1 : ndim, counterDRS))
101 1079152 : stat%avgTimePerFunCall = stat%avgTimePerFunCall + stat%timer%time() - stat%timer%clock
102 : #endif
103 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104 : ! accept or reject the proposed state
105 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 :
107 : #if OMP_ENABLED
108 : do imageID = 1, spec%image%count
109 : #endif
110 1079152 : if (CO_LOGFUNCSTATE(0, imageID, -1) <= CO_LOGFUNCSTATE(0, imageID, counterDRS)) then ! accept the proposed state.
111 142897 : CO_ACCR(imageID, counterDRS) = 1._RKC
112 142897 : CO_ACCR(imageID, -1) = real(counterDRS, RKC)
113 794861 : CO_LOGFUNCSTATE(0 : ndim, imageID, -1) = CO_LOGFUNCSTATE(0 : ndim, imageID, counterDRS)
114 : #if OMP_ENABLED || !CAFMPI_SINGLCHAIN_ENABLED
115 : exit loopNextMove
116 : #endif
117 936255 : elseif (CO_LOGFUNCSTATE(0, imageID, counterDRS) < GET_ELL(maxLogFuncRejectedProposal, imageID)) then ! reject the proposed state. This step should be reachable only when proposalDelayedRejectionCount > 0
118 209101 : CO_ACCR(imageID, counterDRS) = 0._RKC ! proposal rejected XXX is this correct? could `co_accr` be absolute zero?
119 : else ! accept with probability co_accr.
120 727154 : if (counterDRS == 0_IK) then ! This should be equivalent to GET_ELL(maxLogFuncRejectedProposal, imageID) == NEGBIG_RK
121 691831 : logFuncDiff = CO_LOGFUNCSTATE(0, imageID, counterDRS) - CO_LOGFUNCSTATE(0, imageID, -1)
122 691831 : if (logFuncDiff < log(tiny(0._RKC))) then ! xxx should the condition for LOGHUGE_RK be also added?
123 0 : CO_ACCR(imageID, counterDRS) = 0._RKC
124 : else
125 691831 : CO_ACCR(imageID, counterDRS) = exp(logFuncDiff)
126 : end if
127 : else ! ensure no arithmetic overflow/underflow. ATTN: GET_ELL(maxLogFuncRejectedProposal, imageID) < co_logFuncState(0, counterDRS) < co_logFuncState(0, -1)
128 : CO_ACCR(imageID, counterDRS) & ! LCOV_EXCL_LINE
129 : = getLogSubExp(smaller = GET_ELL(maxLogFuncRejectedProposal, imageID), larger = CO_LOGFUNCSTATE(0, imageID, counterDRS)) & ! LCOV_EXCL_LINE
130 35323 : - getLogSubExp(smaller = GET_ELL(maxLogFuncRejectedProposal, imageID), larger = CO_LOGFUNCSTATE(0, imageID, -1))
131 35323 : if (CO_ACCR(imageID, counterDRS) < log(tiny(0._RKC))) then
132 0 : CO_ACCR(imageID, counterDRS) = 0._RKC
133 : else
134 35323 : CO_ACCR(imageID, counterDRS) = exp(CO_ACCR(imageID, counterDRS))
135 : end if
136 : end if
137 727154 : if (GET_ELL(unifrnd, imageID) < CO_ACCR(imageID, counterDRS)) then ! accept the proposed state
138 790286 : CO_LOGFUNCSTATE(0 : ndim, imageID, -1) = CO_LOGFUNCSTATE(0 : ndim, imageID, counterDRS)
139 142140 : CO_ACCR(imageID, -1) = real(counterDRS, RKC)
140 : #if OMP_ENABLED || !CAFMPI_SINGLCHAIN_ENABLED
141 142140 : exit loopNextMove
142 : #endif
143 : end if
144 : end if
145 794115 : GET_ELL(maxLogFuncRejectedProposal, imageID) = max(GET_ELL(maxLogFuncRejectedProposal, imageID), CO_LOGFUNCSTATE(0, imageID, counterDRS))
146 : #if OMP_ENABLED
147 : end do
148 : #endif
149 : else ! if dryrun
150 : #if CAFMPI_SINGLCHAIN_ENABLED
151 : if (spec%image%id == stat%cfc%processID(numFunCallAcceptedPlusOne) .and. &
152 : currentStateWeight + spec%image%id - 1_IK == stat%cfc%sampleWeight(stat%numFunCallAccepted) .and. &
153 : counterDRS == stat%cfc%delayedRejectionStage(numFunCallAcceptedPlusOne)) then
154 : co_accr(-1) = real(counterDRS, RKC)
155 : co_logFuncState(0, counterDRS) = stat%cfc%sampleLogFunc(numFunCallAcceptedPlusOne)
156 : co_logFuncState(1 : ndim, counterDRS) = stat%cfc%sampleState(1 : ndim,numFunCallAcceptedPlusOne)
157 : co_logFuncState(0 : ndim, -1) = co_logFuncState(0 : ndim, counterDRS)
158 : end if
159 : #elif OMP_ENABLED
160 : imageID = stat%cfc%processID(numFunCallAcceptedPlusOne)
161 : if (currentStateWeight + imageID - 1_IK == stat%cfc%sampleWeight(stat%numFunCallAccepted) .and. counterDRS == stat%cfc%delayedRejectionStage(numFunCallAcceptedPlusOne)) then
162 : co_accr(imageID, -1) = real(counterDRS, RKC)
163 : co_logFuncState(0, imageID, counterDRS) = stat%cfc%sampleLogFunc(numFunCallAcceptedPlusOne)
164 : co_logFuncState(1 : ndim, imageID, counterDRS) = stat%cfc%sampleState(1 : ndim,numFunCallAcceptedPlusOne)
165 : co_logFuncState(0 : ndim, imageID, -1) = co_logFuncState(0 : ndim, imageID, counterDRS)
166 : exit loopNextMove
167 : end if
168 : #else
169 0 : if (currentStateWeight == stat%cfc%sampleWeight(stat%numFunCallAccepted) .and. counterDRS == stat%cfc%delayedRejectionStage(numFunCallAcceptedPlusOne)) then
170 0 : co_accr(-1) = real(counterDRS, RKC)
171 0 : co_logFuncState(0, -1) = stat%cfc%sampleLogFunc(numFunCallAcceptedPlusOne)
172 0 : co_logFuncState(1 : ndim, -1) = stat%cfc%sampleState(1 : ndim, numFunCallAcceptedPlusOne)
173 : exit loopNextMove
174 : end if
175 : #endif
176 : end if
177 : #if CAFMPI_SINGLCHAIN_ENABLED
178 : ! broadcast the sampling status from the all images to all images.
179 : ! This is needed in singleChain parallelism to avoid unnecessary delayed rejection if a proposal is already found.
180 : ! Note that this is not necessary when the delayed rejection is deactivated.
181 : if (0 < spec%proposalDelayedRejectionCount%val) then
182 : if (-0.5_RKC < co_accr(-1)) proposalFoundSinglChainMode = 1 ! -0.5 instead of -1. avoids possible real roundoff errors.
183 : stat%timer%clock = stat%timer%time()
184 : SET_CAF(call co_sum(proposalFoundSinglChainMode))
185 : ! send buffer, recv buffer, buffer size, datatype, mpi reduction operation, comm, ierr
186 : SET_MPI(call mpi_allreduce(proposalFoundSinglChainMode, proposalFoundSinglChainModeReduced, 1, mpi_integer, mpi_sum, mpi_comm_world, ierrMPI))
187 : SET_MPI(proposalFoundSinglChainMode = proposalFoundSinglChainModeReduced)
188 : stat%avgCommPerFunCall = stat%avgCommPerFunCall + stat%timer%time() - stat%timer%clock
189 : if (0_IK < proposalFoundSinglChainMode) exit loopNextMove
190 : end if
191 : #endif
192 : end do loopNextMove
|
