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 : #if defined SINGLCHAIN_PARALLELISM
44 : #define LOOP_NEXT_MOVE loopNextMoveSinglChain
45 : #else
46 : #define LOOP_NEXT_MOVE loopNextMove
47 : #endif
48 :
49 : #if defined SINGLCHAIN_PARALLELISM
50 284853 : proposalFoundSinglChainMode = 0_IK
51 : #if defined CAF_ENABLED
52 : ! This syncing is necessary since the co_LogFuncState has to be fetched from the
53 : ! first image by all other images before it is updated again below by the first image.
54 284853 : if (self%Image%isLeader) then
55 94951 : call self%Timer%toc()
56 94951 : sync images(*)
57 94951 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
58 : else
59 189902 : sync images(1)
60 : end if
61 : #endif
62 : #endif
63 :
64 760455 : LOOP_NEXT_MOVE : do counterDRS = 0, self%SpecDRAM%DelayedRejectionCount%val
65 :
66 : #if defined SINGLCHAIN_PARALLELISM
67 : ! self%Stats%NumFunCall%acceptedRejectedDelayedUnused is relevant only on the first image, despite being updated by all images
68 323382 : self%Stats%NumFunCall%acceptedRejectedDelayedUnused = self%Stats%NumFunCall%acceptedRejectedDelayedUnused + self%Image%count
69 : #endif
70 :
71 : co_LogFuncState(1:nd,counterDRS) = self%Proposal%getNew ( nd = nd & ! LCOV_EXCL_LINE
72 : , counterDRS = counterDRS & ! LCOV_EXCL_LINE
73 : , StateOld = co_LogFuncState(1:nd,counterDRS-1) & ! LCOV_EXCL_LINE
74 1168125 : )
75 : #if (defined MATLAB_ENABLED || defined PYTHON_ENABLED || defined R_ENABLED) && !defined CAF_ENABLED && !defined MPI_ENABLED
76 : if(ProposalErr%occurred) then; self%Err%occurred = .true.; self%Err%msg = ProposalErr%msg; return; end if
77 : #elif defined CODECOV_ENABLED || defined SAMPLER_TEST_ENABLED
78 : ! This block is exclusively used to test the deterministic restart functionality of ParaDXXX samplers.
79 : ! This block must not be activated under any other circumstances.
80 : ! This block must be executed by all images.
81 422253 : self%testSamplingCounter = self%testSamplingCounter + 1_IK
82 422253 : self%Err%occurred = self%Err%occurred .or. self%testSamplingCounter > self%testSamplingCountTarget
83 422253 : if (self%Err%occurred) then
84 0 : if (self%testSamplingCounter > self%testSamplingCountTarget) self%Err%msg = "The simulation was interrupted at the requested sampling count for restart testing purposes."
85 : #if defined MPI_ENABLED || defined CAF_ENABLED
86 0 : block; use Err_mod, only: bcastErr; call bcastErr(self%Err); end block
87 : #endif
88 0 : return
89 : end if
90 : #endif
91 :
92 : ! The following random function call is only needed fresh runs to evaluate the acceptance of a proposal.
93 : ! However, it is taken out of the subsequent loops to achieve 100% deterministic reproducibility when
94 : ! a simulation is restarted.
95 :
96 422253 : call random_number(uniformRnd)
97 :
98 494754 : if (self%isFreshRun .or. numFunCallAcceptedPlusOne==self%Chain%Count%compact) then
99 :
100 332621 : call self%Timer%toc()
101 332621 : co_LogFuncState(0,counterDRS) = getLogFunc(nd,co_LogFuncState(1:nd,counterDRS))
102 332621 : call self%Timer%toc(); self%Stats%avgTimePerFunCalInSec = self%Stats%avgTimePerFunCalInSec + self%Timer%Time%delta
103 :
104 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105 : ! accept or reject the proposed state
106 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107 :
108 332621 : if ( co_LogFuncState(0,counterDRS) >= co_LogFuncState(0,-1) ) then ! accept the proposed state
109 :
110 63500 : co_AccRate(counterDRS) = 1._RK
111 63500 : co_AccRate(-1) = real(counterDRS,kind=RK)
112 225894 : co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,counterDRS)
113 : #if !defined SINGLCHAIN_PARALLELISM
114 7425 : exit LOOP_NEXT_MOVE
115 : #endif
116 :
117 269121 : elseif ( co_LogFuncState(0,counterDRS) < maxLogFuncRejectedProposal ) then ! reject the proposed state. This step should be reachable only when delayedRejectionCount > 0
118 :
119 15379 : co_AccRate(counterDRS) = 0._RK ! proposal rejected XXX is this correct? could co_AccRatebe absolute zero?
120 :
121 : else ! accept with probability co_AccRate
122 :
123 253742 : if ( counterDRS == 0_IK ) then ! This should be equivalent to maxLogFuncRejectedProposal == NEGINF_RK
124 249484 : logFuncDiff = co_LogFuncState(0,counterDRS) - co_LogFuncState(0,-1)
125 249484 : if (logFuncDiff < NEGLOGINF_RK) then ! xxx should the condition for LOGHUGE_RK be also added?
126 18757 : co_AccRate(counterDRS) = 0._RK
127 : else
128 230727 : co_AccRate(counterDRS) = exp( logFuncDiff )
129 : end if
130 : else ! ensure no arithmetic overflow/underflow. ATT: co_LogFuncState(0,-1) > co_LogFuncState(0,counterDRS) > maxLogFuncRejectedProposal
131 : co_AccRate(counterDRS) = getLogSubExp( co_LogFuncState(0,counterDRS) , maxLogFuncRejectedProposal ) & ! LCOV_EXCL_LINE
132 4258 : - getLogSubExp( co_LogFuncState(0,-1) , maxLogFuncRejectedProposal )
133 4258 : if (co_AccRate(counterDRS) < NEGLOGINF_RK) then
134 2 : co_AccRate(counterDRS) = 0._RK
135 : else
136 4256 : co_AccRate(counterDRS) = exp(co_AccRate(counterDRS))
137 : end if
138 : end if
139 :
140 253742 : if (uniformRnd<co_AccRate(counterDRS)) then ! accept the proposed state
141 62009 : co_AccRate(-1) = real(counterDRS,kind=RK)
142 220384 : co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,counterDRS)
143 : #if !defined SINGLCHAIN_PARALLELISM
144 7581 : exit LOOP_NEXT_MOVE
145 : #endif
146 : end if
147 :
148 : end if
149 :
150 317615 : maxLogFuncRejectedProposal = max( maxLogFuncRejectedProposal, co_LogFuncState(0,counterDRS) )
151 :
152 : else ! if dryrun
153 :
154 : #if defined SINGLCHAIN_PARALLELISM
155 87798 : if (self%Image%id==self%Chain%ProcessID(numFunCallAcceptedPlusOne) .and. &
156 87798 : currentStateWeight+self%Image%id-1_IK==self%Chain%Weight(self%Stats%NumFunCall%accepted) .and. &
157 43899 : counterDRS==self%Chain%DelRejStage(numFunCallAcceptedPlusOne)) then
158 6972 : co_AccRate(-1) = real(counterDRS,kind=RK)
159 6972 : co_LogFuncState( 0,counterDRS) = self%Chain%LogFunc (numFunCallAcceptedPlusOne)
160 20916 : co_LogFuncState(1:nd,counterDRS) = self%Chain%State(1:nd,numFunCallAcceptedPlusOne)
161 27888 : co_LogFuncState(0:nd,-1) = co_LogFuncState(0:nd,counterDRS)
162 : end if
163 : #else
164 45733 : if (currentStateWeight==self%Chain%Weight(self%Stats%NumFunCall%accepted) .and. counterDRS==self%Chain%DelRejStage(numFunCallAcceptedPlusOne)) then
165 11364 : co_AccRate(-1) = real(counterDRS,kind=RK)
166 11364 : co_LogFuncState( 0,-1) = self%Chain%LogFunc (numFunCallAcceptedPlusOne)
167 34092 : co_LogFuncState(1:nd,-1) = self%Chain%State(1:nd,numFunCallAcceptedPlusOne)
168 11364 : exit LOOP_NEXT_MOVE
169 : end if
170 : #endif
171 : end if
172 :
173 : #if defined SINGLCHAIN_PARALLELISM
174 : ! broadcast the sampling status from the all images to all images.
175 : ! This is needed in singlChain parallelism to avoid unnecessary delayed rejection if a proposal is already found.
176 : ! Note that this is not necessary when the delayed rejection is deactivated.
177 608235 : if (delayedRejectionRequested) then
178 65469 : if (co_AccRate(-1) > -0.5_RK) proposalFoundSinglChainMode = 1_IK ! -0.5 instead of -1. avoids possible real roundoff errors.
179 65469 : call self%Timer%toc()
180 : #if defined CAF_ENABLED
181 65469 : call co_sum(proposalFoundSinglChainMode)
182 : #elif defined MPI_ENABLED
183 : call mpi_allreduce ( proposalFoundSinglChainMode & ! send buffer
184 : , proposalFoundSinglChainModeReduced & ! recv buffer
185 : , 1 & ! buffer size
186 : , mpi_integer & ! datatype
187 : , mpi_sum & ! mpi reduction operation
188 : , mpi_comm_world & ! comm
189 : , ierrMPI & ! ierr
190 : )
191 : proposalFoundSinglChainMode = proposalFoundSinglChainModeReduced
192 : #endif
193 65469 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
194 65469 : if (proposalFoundSinglChainMode>0_IK) exit LOOP_NEXT_MOVE
195 : end if
196 : #endif
197 :
198 : end do LOOP_NEXT_MOVE
199 :
200 : #if defined SINGLCHAIN_PARALLELISM && defined MPI_ENABLED
201 : ! This SINGLCHAIN_PARALLELISM preprocessing avoids unnecessary communication in serial or multichain-parallel modes.
202 : ! gather all accr on the leader image.
203 : call self%Timer%toc()
204 : call mpi_gather ( co_AccRate(:) & ! send buffer
205 : , delayedRejectionCountPlusTwo & ! send count
206 : , mpi_double_precision & ! send datatype
207 : , AccRateMatrix(:,:) & ! recv buffer
208 : , delayedRejectionCountPlusTwo & ! recv count
209 : , mpi_double_precision & ! recv datatype
210 : , 0 & ! root rank
211 : , mpi_comm_world & ! comm
212 : , ierrMPI & ! mpi error flag
213 : )
214 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
215 : #elif defined SINGLCHAIN_PARALLELISM && defined CAF_ENABLED
216 284853 : if (self%Image%isLeader) then
217 94951 : call self%Timer%toc()
218 94951 : sync images(*)
219 94951 : call self%Timer%toc(); self%Stats%avgCommTimePerFunCall = self%Stats%avgCommTimePerFunCall + self%Timer%Time%delta
220 : else
221 189902 : sync images(1)
222 : end if
223 : #endif
224 :
225 : #undef LOOP_NEXT_MOVE
|
