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 Proposal modules of the `ParaDRAM` and `ParaDISE` samplers.
45 : !>
46 : !> \remark
47 : !> This module requires preprocessing, prior to compilation.
48 : !>
49 : !> \author Amir Shahmoradi
50 :
51 : #if defined UNIFORM
52 :
53 : #define GET_RANDOM_PROPOSAL getRandMVU
54 :
55 : #elif defined NORMAL
56 :
57 : #define GET_RANDOM_PROPOSAL getRandMVN
58 :
59 : #else
60 :
61 : #error "Unknown Proposal model in ParaDXXXProposal_mod.inc.f90"
62 :
63 : #endif
64 :
65 :
66 : #if defined PARADRAM
67 :
68 : #define ParaDXXX ParaDRAM
69 : use ParaDRAMProposalAbstract_mod, only: ProposalAbstract_type, ProposalErr
70 :
71 : #elif defined PARADISE
72 :
73 : #define ParaDXXX ParaDISE
74 : use ParaDISEProposalAbstract_mod, only: ProposalAbstract_type, ProposalErr
75 :
76 : #endif
77 :
78 : use ParaMonte_mod, only: Image_type
79 : use Constants_mod, only: IK, RK, PMSM
80 : use String_mod, only: IntStr_type
81 :
82 : implicit none
83 :
84 : private
85 : public :: Proposal_type
86 :
87 : #if defined NORMAL
88 : character(*), parameter :: MODULE_NAME = "@"//PMSM%ParaDXXX//"ProposalNormal_mod"
89 : #elif defined UNIFORM
90 : character(*), parameter :: MODULE_NAME = "@"//PMSM%ParaDXXX//"ProposalUniform_mod"
91 : #endif
92 :
93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94 :
95 : type :: AccRate_type
96 : real(RK) :: sumUpToLastUpdate
97 : real(RK) :: target
98 : end type AccRate_type
99 :
100 : !> The `Proposal_type` class.
101 : type, extends(ProposalAbstract_type) :: Proposal_type
102 : !type(AccRate_type) :: AccRate
103 : contains
104 : procedure , nopass :: getNew
105 : #if defined PARADISE
106 : procedure , nopass :: getLogProb
107 : #endif
108 : procedure , nopass :: doAdaptation
109 : !procedure , nopass :: readRestartFile
110 : !procedure , nopass :: writeRestartFile
111 : #if defined CAF_ENABLED || defined MPI_ENABLED
112 : procedure , nopass :: bcastAdaptation
113 : #endif
114 : end type Proposal_type
115 :
116 : interface Proposal_type
117 : module procedure :: constructProposalSymmetric
118 : end interface Proposal_type
119 :
120 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 :
122 : ! Covariance Matrix of the proposal distribution. Last index belongs to delayed rejection
123 :
124 : #if defined CAF_ENABLED
125 : real(RK) , save , allocatable :: comv_CholDiagLower(:,:,:)[:]
126 : #else
127 : real(RK) , save , allocatable :: comv_CholDiagLower(:,:,:)
128 : #endif
129 : real(RK) , save , allocatable :: mv_logSqrtDetInvCovMat(:)
130 : real(RK) , save , allocatable :: mv_InvCovMat(:,:,:)
131 :
132 : #if defined MPI_ENABLED
133 : integer(IK) , save :: mc_ndimSqPlusNdim
134 : #endif
135 : type(Image_type), save :: mc_Image
136 : integer(IK) , save :: mc_ndim
137 : integer(IK) , save :: mc_logFileUnit
138 : integer(IK) , save :: mc_restartFileUnit
139 : logical , save :: mc_scalingRequested
140 : real(RK) , save :: mc_defaultScaleFactorSq
141 : integer(IK) , save :: mc_DelayedRejectionCount
142 : integer(IK) , save :: mc_MaxNumDomainCheckToWarn
143 : integer(IK) , save :: mc_MaxNumDomainCheckToStop
144 : logical , save :: mc_delayedRejectionRequested
145 : real(RK) , save :: mc_ndimInverse
146 : real(RK) , save :: mc_targetAcceptanceRate
147 : real(RK) , save :: mc_TargetAcceptanceRateLimit(2)
148 : !real(RK) , save :: mc_maxScaleFactor != 2._RK
149 : !real(RK) , save :: mc_maxScaleFactorSq != mc_maxScaleFactor**2
150 : real(RK) , save , allocatable :: mc_DelayedRejectionScaleFactorVec(:)
151 : real(RK) , save , allocatable :: mc_DomainLowerLimitVec(:)
152 : real(RK) , save , allocatable :: mc_DomainUpperLimitVec(:)
153 : logical , save , allocatable :: mc_isAsciiRestartFileFormat
154 : logical , save , allocatable :: mc_isBinaryRestartFileFormat
155 : character(:) , save , allocatable :: mc_MaxNumDomainCheckToWarnMsg
156 : character(:) , save , allocatable :: mc_MaxNumDomainCheckToStopMsg
157 : character(:) , save , allocatable :: mc_negativeTotalVariationMsg
158 : character(:) , save , allocatable :: mc_restartFileFormat
159 : character(:) , save , allocatable :: mc_methodBrand
160 : character(:) , save , allocatable :: mc_methodName
161 : #if defined UNIFORM
162 : real(RK) , save , allocatable :: mc_negLogVolUnitBall
163 : #endif
164 :
165 : ! the following had to be defined globally for the sake of restart file generation
166 :
167 : real(RK) , save , allocatable :: mv_MeanOld_save(:)
168 : real(RK) , save :: mv_logSqrtDetOld_save
169 : real(RK) , save :: mv_adaptiveScaleFactorSq_save ! = 1._RK
170 : integer(IK) , save :: mv_sampleSizeOld_save ! = 0_IK
171 :
172 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
173 :
174 : contains
175 :
176 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
177 :
178 : !> This is the constructor of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
179 : !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
180 : !> Return an object of [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
181 : !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes given the input simulation characteristics.
182 : !>
183 : !> \remark
184 : !> This interface madness is a result of the internal compiler bug in GFortran as of Jan 2020, which diagnoses a `ParaDXXX_type`
185 : !> argument as circular dependency due to this constructor appearing in the type-bound setup procedure of `ParaDXXX_type`.
186 : !> Intel does not complain. Until GFortran comes up with a fix, we have to live with this interface.
187 687 : function constructProposalSymmetric ( ndim &
188 : , SpecBase &
189 : , SpecMCMC &
190 : , SpecDRAM &
191 : , Image &
192 : , name &
193 : , brand &
194 : , LogFile &
195 : , RestartFile &
196 : , isFreshRun &
197 : ) result(Proposal)
198 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
199 : !DEC$ ATTRIBUTES DLLEXPORT :: constructProposalSymmetric
200 : #endif
201 : use Constants_mod, only: IK, RK, NULL_RK
202 : use ParaMonte_mod, only: Image_type
203 : use ParaMonte_mod, only: LogFile_type
204 : use ParaMonte_mod, only: RestartFile_type
205 : use SpecBase_mod, only: SpecBase_type
206 : use SpecMCMC_mod, only: SpecMCMC_type
207 : use SpecDRAM_mod, only: SpecDRAM_type
208 : use Matrix_mod, only: getCholeskyFactor
209 : use String_mod, only: num2str
210 : use Err_mod, only: abort
211 : #if defined UNIFORM
212 : use Math_mod, only: getLogVolUnitBall
213 : #endif
214 : implicit none
215 :
216 : integer(IK) , intent(in) :: ndim
217 : type(SpecBase_type) , intent(in) :: SpecBase
218 : type(SpecMCMC_type) , intent(in) :: SpecMCMC
219 : type(SpecDRAM_type) , intent(in) :: SpecDRAM
220 : type(Image_type) , intent(in) :: Image
221 : character(*) , intent(in) :: name
222 : character(*) , intent(in) :: brand
223 : type(LogFile_type) , intent(in) :: LogFile
224 : type(RestartFile_type) , intent(in) :: RestartFile
225 : logical :: isFreshRun
226 :
227 : type(Proposal_type) :: Proposal
228 :
229 : character(*), parameter :: PROCEDURE_NAME = MODULE_NAME // "@constructProposalSymmetric()"
230 : integer :: i, j
231 :
232 687 : ProposalErr%occurred = .false.
233 :
234 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
235 : ! setup sampler update global save variables
236 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
237 :
238 687 : if (allocated(mv_MeanOld_save)) deallocate(mv_MeanOld_save); allocate(mv_MeanOld_save(ndim))
239 1749 : mv_MeanOld_save(1:ndim) = SpecMCMC%StartPointVec%Val
240 687 : mv_logSqrtDetOld_save = NULL_RK
241 687 : mv_sampleSizeOld_save = 1_IK
242 687 : mv_adaptiveScaleFactorSq_save = 1._RK
243 :
244 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
245 : ! setup general proposal specifications
246 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
247 :
248 : #if defined MPI_ENABLED
249 : mc_ndimSqPlusNdim = ndim*(ndim+1_IK)
250 : #endif
251 687 : mc_ndim = ndim
252 2319 : mc_DomainLowerLimitVec = SpecBase%DomainLowerLimitVec%Val
253 2319 : mc_DomainUpperLimitVec = SpecBase%DomainUpperLimitVec%Val
254 1473 : mc_DelayedRejectionScaleFactorVec = SpecDRAM%delayedRejectionScaleFactorVec%Val
255 : !mc_isNormal = SpecMCMC%ProposalModel%isNormal
256 687 : mc_Image = Image
257 687 : mc_methodName = name
258 687 : mc_methodBrand = brand
259 687 : mc_logFileUnit = LogFile%unit
260 687 : mc_restartFileUnit = RestartFile%unit
261 687 : mc_restartFileFormat = RestartFile%format
262 687 : mc_isBinaryRestartFileFormat = SpecBase%RestartFileFormat%isBinary
263 687 : mc_isAsciiRestartFileFormat = SpecBase%RestartFileFormat%isAscii
264 687 : mc_defaultScaleFactorSq = SpecMCMC%ScaleFactor%val**2
265 : !Proposal%AccRate%sumUpToLastUpdate = 0._RK
266 687 : mc_maxNumDomainCheckToWarn = SpecBase%MaxNumDomainCheckToWarn%val
267 687 : mc_maxNumDomainCheckToStop = SpecBase%MaxNumDomainCheckToStop%val
268 687 : mc_delayedRejectionCount = SpecDRAM%DelayedRejectionCount%val
269 687 : mc_delayedRejectionRequested = mc_DelayedRejectionCount > 0_IK
270 687 : mc_scalingRequested = SpecBase%TargetAcceptanceRate%scalingRequested
271 687 : mc_ndimInverse = 1._RK/real(ndim,kind=RK)
272 2061 : mc_TargetAcceptanceRateLimit = SpecBase%TargetAcceptanceRate%Val
273 687 : mc_targetAcceptanceRate = 0.234_RK
274 795 : if (mc_scalingRequested) mc_targetAcceptanceRate = sum(mc_TargetAcceptanceRateLimit) / 2._RK
275 : !mc_maxScaleFactorSq = 4._RK**mc_ndimInverse
276 : !mc_maxScaleFactor = sqrt(mc_maxScaleFactorSq)
277 :
278 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
279 : ! setup ProposalSymmetric specifications
280 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281 :
282 : ! setup covariance matrix
283 :
284 687 : if (allocated(mv_InvCovMat)) deallocate(mv_InvCovMat)
285 687 : allocate( mv_InvCovMat(ndim,0:ndim,0:mc_DelayedRejectionCount) )
286 687 : if (allocated(mv_logSqrtDetInvCovMat)) deallocate(mv_logSqrtDetInvCovMat)
287 687 : allocate( mv_logSqrtDetInvCovMat(0:mc_DelayedRejectionCount) )
288 :
289 687 : if (allocated(comv_CholDiagLower)) deallocate(comv_CholDiagLower)
290 : #if defined CAF_ENABLED
291 : ! on the second dimension, the zeroth index refers to the Diagonal elements of the Cholesky lower triangular matrix
292 : ! This rearrangement was done for more efficient communication of the matrix across processes.
293 687 : allocate( comv_CholDiagLower(ndim,0:ndim,0:mc_DelayedRejectionCount)[*] )
294 : #else
295 : allocate( comv_CholDiagLower(ndim,0:ndim,0:mc_DelayedRejectionCount) )
296 : #endif
297 :
298 3561 : comv_CholDiagLower(1:ndim,1:ndim,0) = SpecMCMC%ProposalStartCovMat%Val
299 :
300 : ! Now scale the covariance matrix
301 :
302 1749 : do j = 1, ndim
303 3186 : do i = 1, j
304 2499 : comv_CholDiagLower(i,j,0) = comv_CholDiagLower(i,j,0) * mc_defaultScaleFactorSq
305 : end do
306 : end do
307 :
308 : ! Now get the Cholesky Factor of the Covariance Matrix. Lower comv_CholDiagLower will be the CholFac
309 :
310 687 : call getCholeskyFactor( ndim, comv_CholDiagLower(:,1:ndim,0), comv_CholDiagLower(1:ndim,0,0) ) ! The `:` instead of `1:ndim` avoids temporary array creation.
311 687 : if (comv_CholDiagLower(1,0,0)<0._RK) then
312 : ! LCOV_EXCL_START
313 : ProposalErr%msg = mc_Image%name // PROCEDURE_NAME // ": Singular input covariance matrix by user was detected. This is strange.\nCovariance matrix lower triangle:"
314 : do j = 1, ndim
315 : do i = 1, j
316 : ProposalErr%msg = ProposalErr%msg // "\n" // num2str(comv_CholDiagLower(1:i,j,0))
317 : end do
318 : end do
319 : ProposalErr%occurred = .true.
320 : call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
321 : return
322 : ! LCOV_EXCL_STOP
323 : end if
324 :
325 687 : if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
326 687 : call getInvCovMat()
327 1749 : mv_logSqrtDetOld_save = sum(log( comv_CholDiagLower(1:ndim,0,0) ))
328 :
329 : ! Scale the higher-stage delayed-rejection Cholesky Lower matrices
330 :
331 687 : if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
332 :
333 : ! This will be used for Domain boundary checking during the simulation
334 :
335 687 : mc_negativeTotalVariationMsg = MODULE_NAME//"@doAdaptation(): Non-positive adaptation measure detected, possibly due to round-off error: "
336 : mc_MaxNumDomainCheckToWarnMsg = MODULE_NAME//"@getNew(): "//num2str(mc_MaxNumDomainCheckToWarn)//&
337 687 : " proposals were drawn out of the objective function's Domain without any acceptance."
338 : mc_MaxNumDomainCheckToStopMsg = MODULE_NAME//"@getNew(): "//num2str(mc_MaxNumDomainCheckToStop)//&
339 : " proposals were drawn out of the objective function's Domain. As per the value set for the&
340 687 : & simulation specification variable 'maxNumDomainCheckToStop', "//mc_methodName//" will abort now."
341 :
342 : #if defined UNIFORM
343 120 : if (ndim>1_IK) then
344 108 : mc_negLogVolUnitBall = -getLogVolUnitBall(ndim)
345 : else
346 12 : mc_negLogVolUnitBall = 0._RK
347 : end if
348 : #endif
349 :
350 : ! read/write the first entry of the restart file
351 :
352 687 : if (mc_Image%isLeader) then
353 : block
354 269 : real(RK) :: meanAccRateSinceStart
355 269 : if (isFreshRun) then
356 237 : call writeRestartFile(meanAccRateSinceStart=1._RK)
357 237 : call writeRestartFile()
358 : else
359 32 : call readRestartFile(meanAccRateSinceStart)
360 32 : call readRestartFile()
361 : end if
362 : end block
363 : end if
364 :
365 687 : end function constructProposalSymmetric
366 :
367 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
368 :
369 : !> \brief
370 : !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
371 : !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
372 : !> Return a newly sampled state to the kernel routine for acceptance check.
373 : !>
374 : !>
375 : !> @param[in] nd : The number of dimensions of the objective function.
376 : !> @param[in] counterDRS : The Delayed Rejection Stage counter.
377 : !> @param[in] StateOld : The old sampled state.
378 : !>
379 : !> \return
380 : !> `StateNew` : The newly sampled state.
381 422253 : function getNew ( nd &
382 : , counterDRS &
383 422253 : , StateOld &
384 1590378 : ) result (StateNew)
385 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
386 : !DEC$ ATTRIBUTES DLLEXPORT :: getNew
387 : #endif
388 687 : use Statistics_mod, only: GET_RANDOM_PROPOSAL
389 : use Constants_mod, only: IK, RK
390 : use Err_mod, only: warn, abort
391 : use ParaMonteLogFunc_mod, only: getLogFunc_proc
392 :
393 : implicit none
394 :
395 : character(*), parameter :: PROCEDURE_NAME = MODULE_NAME // "@getNew()"
396 :
397 : integer(IK), intent(in) :: nd
398 : integer(IK), intent(in) :: counterDRS
399 : real(RK) , intent(in) :: StateOld(nd)
400 : real(RK) :: StateNew(nd)
401 : integer(IK) :: domainCheckCounter
402 :
403 422253 : domainCheckCounter = 0_IK
404 :
405 62428 : loopBoundaryCheck: do ! Check for the support Region consistency:
406 : #if defined UNIFORM || defined NORMAL
407 : StateNew(1:nd) = GET_RANDOM_PROPOSAL( nd & ! LCOV_EXCL_LINE
408 : , StateOld & ! LCOV_EXCL_LINE
409 : ! ATTN: The colon index in place of 1:nd below avoids the temporary array creation
410 : , comv_CholDiagLower(:,1:nd,counterDRS) & ! LCOV_EXCL_LINE
411 : , comv_CholDiagLower(1:nd,0,counterDRS) & ! LCOV_EXCL_LINE
412 484681 : )
413 : #endif
414 2148107 : if ( any(StateNew(1:nd)<=mc_DomainLowerLimitVec) .or. any(StateNew(1:nd)>=mc_DomainUpperLimitVec) ) then
415 62446 : domainCheckCounter = domainCheckCounter + 1
416 62473 : if (domainCheckCounter==mc_MaxNumDomainCheckToWarn) then
417 27 : call warn( prefix = mc_methodBrand, outputUnit = mc_logFileUnit, msg = mc_MaxNumDomainCheckToWarnMsg )
418 : end if
419 62446 : if (domainCheckCounter==mc_MaxNumDomainCheckToStop) then
420 18 : ProposalErr%occurred = .true.
421 18 : ProposalErr%msg = mc_MaxNumDomainCheckToStopMsg
422 : #if !defined CODECOV_ENABLED && ((!defined MATLAB_ENABLED && !defined PYTHON_ENABLED && !defined R_ENABLED) || defined CAF_ENABLED && defined MPI_ENABLED )
423 : call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
424 : #endif
425 18 : return
426 : end if
427 62428 : cycle loopBoundaryCheck
428 : end if
429 422235 : exit loopBoundaryCheck
430 : end do loopBoundaryCheck
431 :
432 844506 : end function getNew
433 :
434 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
435 :
436 : #if defined PARADISE
437 : !> \brief
438 : !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
439 : !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
440 : !> Return the natural logarithm of the probability of acceptance of the new state given the old state.
441 : !>
442 : !> @param[in] nd : The number of dimensions of the objective function.
443 : !> @param[in] counterDRS : The Delayed Rejection Stage counter.
444 : !> @param[in] StateOld : The old sampled state.
445 : !> @param[in] StateNew : The new sampled state.
446 : !>
447 : !> \return
448 : !> `logProb` : The log probability of obtaining obtaining the new sample given the old sample.
449 : ! LCOV_EXCL_START
450 : pure function getLogProb( nd &
451 : , counterDRS &
452 : , StateOld &
453 : , StateNew &
454 : ) result(logProb)
455 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
456 : !DEC$ ATTRIBUTES DLLEXPORT :: getLogProb
457 : #endif
458 : #if defined NORMAL
459 : use Statistics_mod, only: getLogProbMVN
460 : #elif defined UNIFORM
461 : use Statistics_mod, only: isInsideEllipsoid
462 : #endif
463 : use Constants_mod, only: IK, RK, NEGINF_RK
464 : implicit none
465 : integer(IK), intent(in) :: nd
466 : integer(IK), intent(in) :: counterDRS
467 : real(RK) , intent(in) :: StateOld(nd)
468 : real(RK) , intent(in) :: StateNew(nd)
469 : real(RK) :: logProb
470 : #if defined UNIFORM
471 : if (isInsideEllipsoid(nd,StateNew-StateOld,mv_InvCovMat(1:mc_ndim,1:mc_ndim,counterDRS))) then
472 : logProb = mc_negLogVolUnitBall + mv_logSqrtDetInvCovMat(counterDRS)
473 : else
474 : logProb = NEGINF_RK
475 : end if
476 : #elif defined NORMAL
477 : logProb = getLogProbMVN ( nd = nd &
478 : , MeanVec = StateOld &
479 : , InvCovMat = mv_InvCovMat(1:nd,1:nd,counterDRS) &
480 : , logSqrtDetInvCovMat = mv_logSqrtDetInvCovMat(counterDRS) &
481 : , Point = StateNew &
482 : )
483 : #endif
484 : end function getLogProb
485 : ! LCOV_EXCL_STOP
486 : #endif
487 :
488 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
489 :
490 : !> \brief
491 : !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
492 : !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
493 : !> Perform adaptation of the proposal distribution of the MCMC sampler.
494 : !>
495 : !> @param[in] nd : The number of dimensions of the objective function.
496 : !> @param[in] chainSize : The size of the input chain in compact format.
497 : !> @param[in] Chain : The two dimensional `(nd, chainSize)` array of accepted states.
498 : !> @param[in] ChainWeight : The one dimensional integer array of the weights of the sampled states in the chain.
499 : !> @param[in] isFreshRun : The logical flag indicating whether the simulation is in regular (non-restart) mode.
500 : !> @param[in] samplerUpdateIsGreedy : The logical flag indicating whether the adaptation is in greedy mode..
501 : !> @param[inout] meanAccRateSinceStart : The mean acceptance rate since the start of the sampling.
502 : !> @param[out] samplerUpdateSucceeded : The output logical flag indicating whether the sampling succeeded.
503 : !> @param[out] adaptationMeasure : The output real number in the range `[0,1]` indicating the amount of adaptation,
504 : !> with zero indicating no adaptation and one indicating extreme adaptation to the extent
505 : !> that the new adapted proposal distribution is completely different from the previous proposal.
506 : !> \warning
507 : !> This routine must be exclusively called by the leader images.
508 : !>
509 : !> \remark
510 : !> For information on the meaning of `adaptationMeasure`, see the paper by Shahmoradi and Bagheri, 2020, whose PDF is available at:
511 : !> [https://www.cdslab.org/paramonte/notes/overview/preface/#the-paradram-sampler](https://www.cdslab.org/paramonte/notes/overview/preface/#the-paradram-sampler)
512 124416 : subroutine doAdaptation ( nd &
513 : , chainSize &
514 124416 : , Chain &
515 124416 : , ChainWeight &
516 : , isFreshRun &
517 : , samplerUpdateIsGreedy &
518 : , meanAccRateSinceStart &
519 : , samplerUpdateSucceeded &
520 : , adaptationMeasure &
521 : )
522 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
523 : !DEC$ ATTRIBUTES DLLEXPORT :: doAdaptation
524 : #endif
525 :
526 : use Statistics_mod, only: getSamCovUpperMeanTrans, getWeiSamCovUppMeanTrans, mergeMeanCovUpper ! LCOV_EXCL_LINE
527 : use Matrix_mod, only: getCholeskyFactor, getLogSqrtDetPosDefMat
528 : use Constants_mod, only: RK, IK ! , EPS_RK
529 : use String_mod, only: num2str
530 : use Err_mod, only: abort, warn
531 : implicit none
532 :
533 : character(*), parameter :: PROCEDURE_NAME = MODULE_NAME // "@doAdaptation()"
534 :
535 : integer(IK), intent(in) :: nd
536 : integer(IK), intent(in) :: chainSize
537 : real(RK) , intent(in) :: Chain(nd,chainSize)
538 : integer(IK), intent(in) :: ChainWeight(chainSize)
539 : logical , intent(in) :: isFreshRun
540 : logical , intent(in) :: samplerUpdateIsGreedy
541 : real(RK) , intent(inout) :: meanAccRateSinceStart ! is intent(out) in restart mode, intent(in) in fresh mode.
542 : logical , intent(out) :: samplerUpdateSucceeded
543 : real(RK) , intent(out) :: adaptationMeasure
544 :
545 : integer(IK) :: i, j
546 359381 : real(RK) :: MeanNew(nd)
547 483797 : real(RK) :: MeanCurrent(nd)
548 939860 : real(RK) :: CovMatUpperOld(nd,nd)
549 939860 : real(RK) :: CovMatUpperCurrent(nd,nd)
550 124416 : real(RK) :: logSqrtDetNew
551 124416 : real(RK) :: logSqrtDetSum
552 124416 : real(RK) :: adaptiveScaleFactor
553 : logical :: adaptationMeasureComputationNeeded ! only used to avoid redundant affinity computation, if no update occurs
554 : logical :: singularityOccurred, scalingNeeded
555 : integer(IK) :: sampleSizeOld, sampleSizeCurrent
556 :
557 124416 : scalingNeeded = .false.
558 124416 : sampleSizeOld = mv_sampleSizeOld_save ! this is kept only for restoration of mv_sampleSizeOld_save, if needed.
559 :
560 : ! read/write meanAccRateSinceStart from/to restart file
561 :
562 124416 : if (mc_Image%isLeader) then
563 124416 : if (isFreshRun) then
564 75846 : call writeRestartFile(meanAccRateSinceStart)
565 : else
566 48570 : call readRestartFile(meanAccRateSinceStart)
567 : end if
568 : end if
569 :
570 : ! First if there are less than nd+1 points for new covariance computation, then just scale the covariance and return
571 :
572 124416 : blockSufficientSampleSizeCheck: if (chainSize>nd) then
573 :
574 : ! get the upper covariance matrix and Mean of the new sample
575 :
576 30428 : if (samplerUpdateIsGreedy) then
577 42 : sampleSizeCurrent = chainSize
578 42 : call getSamCovUpperMeanTrans(sampleSizeCurrent,nd,Chain,CovMatUpperCurrent,MeanCurrent)
579 : else
580 138478 : sampleSizeCurrent = sum(ChainWeight)
581 30386 : call getWeiSamCovUppMeanTrans(chainSize,sampleSizeCurrent,nd,Chain,ChainWeight,CovMatUpperCurrent,MeanCurrent)
582 : end if
583 :
584 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
585 :
586 : ! combine old and new covariance matrices if both exist
587 :
588 30428 : blockMergeCovMat: if (mv_sampleSizeOld_save==1_IK) then
589 :
590 : ! There is no prior old Covariance matrix to combine with the new one from the new chain
591 :
592 667 : mv_MeanOld_save(1:nd) = MeanCurrent
593 258 : mv_sampleSizeOld_save = sampleSizeCurrent
594 :
595 : ! copy and then scale the new covariance matrix by the default scale factor, which will be then used to get the Cholesky Factor
596 :
597 667 : do j = 1, nd
598 1227 : do i = 1, j
599 560 : CovMatUpperOld(i,j) = comv_CholDiagLower(i,j,0) ! This will be used to recover the old covariance in case of update failure, and to compute the adaptation measure
600 969 : comv_CholDiagLower(i,j,0) = CovMatUpperCurrent(i,j) * mc_defaultScaleFactorSq
601 : end do
602 : end do
603 :
604 : else blockMergeCovMat
605 :
606 : ! first scale the new covariance matrix by the default scale factor, which will be then used to get the Cholesky Factor
607 :
608 78896 : do j = 1, nd
609 146178 : do i = 1, j
610 67282 : CovMatUpperOld(i,j) = comv_CholDiagLower(i,j,0) ! This will be used to recover the old covariance in case of update failure, and to compute the adaptation measure
611 116008 : CovMatUpperCurrent(i,j) = CovMatUpperCurrent(i,j) * mc_defaultScaleFactorSq
612 : end do
613 : end do
614 :
615 : ! now combine it with the old covariance matrix.
616 : ! Do not set the full boundaries' range `(1:nd)` for `comv_CholDiagLower` in the following subroutine call.
617 : ! Setting the boundaries forces the compiler to generate a temporary array.
618 :
619 : call mergeMeanCovUpper ( nd = nd & ! LCOV_EXCL_LINE
620 : , npA = mv_sampleSizeOld_save & ! LCOV_EXCL_LINE
621 : , MeanVecA = mv_MeanOld_save & ! LCOV_EXCL_LINE
622 : , CovMatUpperA = CovMatUpperOld & ! LCOV_EXCL_LINE
623 : , npB = sampleSizeCurrent & ! LCOV_EXCL_LINE
624 : , MeanVecB = MeanCurrent & ! LCOV_EXCL_LINE
625 : , CovMatUpperB = CovMatUpperCurrent & ! LCOV_EXCL_LINE
626 : , MeanVecAB = MeanNew & ! LCOV_EXCL_LINE
627 : , CovMatUpperAB = comv_CholDiagLower(:,1:nd,0) & ! LCOV_EXCL_LINE
628 30170 : )
629 78896 : mv_MeanOld_save(1:nd) = MeanNew
630 :
631 : end if blockMergeCovMat
632 :
633 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
634 :
635 : ! now get the Cholesky factorization
636 :
637 : ! WARNING: Do not set the full boundaries' range `(1:nd)` for the first index of `comv_CholDiagLower` in the following subroutine call.
638 : ! WARNING: Setting the boundaries forces the compiler to generate a temporary array.
639 :
640 30428 : call getCholeskyFactor( nd, comv_CholDiagLower(:,1:nd,0), comv_CholDiagLower(1:nd,0,0) )
641 :
642 30428 : blockPosDefCheck: if (comv_CholDiagLower(1,0,0)>0._RK) then
643 :
644 : !singularityOccurred = .false.
645 30428 : samplerUpdateSucceeded = .true.
646 30428 : adaptationMeasureComputationNeeded = .true.
647 30428 : mv_sampleSizeOld_save = mv_sampleSizeOld_save + sampleSizeCurrent
648 :
649 : ! LCOV_EXCL_START
650 : else blockPosDefCheck
651 :
652 : adaptationMeasure = 0._RK
653 : !singularityOccurred = .true.
654 : samplerUpdateSucceeded = .false.
655 : adaptationMeasureComputationNeeded = .false.
656 : mv_sampleSizeOld_save = sampleSizeOld
657 :
658 : ! it may be a good idea to add a warning message printed out here for the singularity occurrence
659 :
660 : call warn ( prefix = mc_methodBrand & ! LCOV_EXCL_LINE
661 : , outputUnit = mc_logFileUnit & ! LCOV_EXCL_LINE
662 : , marginTop = 0_IK & ! LCOV_EXCL_LINE
663 : , marginBot = 0_IK & ! LCOV_EXCL_LINE
664 : , msg = "Singularity occurred while updating the proposal distribution's covariance matrix." & ! LCOV_EXCL_LINE
665 : )
666 :
667 : ! recover the old upper covariance matrix
668 :
669 : do j = 1, nd
670 : do i = 1, j
671 : comv_CholDiagLower(i,j,0) = CovMatUpperOld(i,j)
672 : end do
673 : end do
674 :
675 : ! ensure the old Cholesky factorization can be recovered
676 :
677 : ! WARNING: Do not set the full boundaries' range `(1:nd)` for the first index of `comv_CholDiagLower` in the following subroutine call.
678 : ! WARNING: Setting the boundaries forces the compiler to generate a temporary array.
679 :
680 : call getCholeskyFactor( nd, comv_CholDiagLower(:,1:nd,0), comv_CholDiagLower(1:nd,0,0) ) ! avoid temporary array creation by using : indexer
681 : if (comv_CholDiagLower(1,0,0)<0._RK) then
682 : write(mc_logFileUnit,"(A)")
683 : write(mc_logFileUnit,"(A)") "Singular covariance matrix detected:"
684 : write(mc_logFileUnit,"(A)")
685 : do j = 1, nd
686 : write(mc_logFileUnit,"(*(E25.15))") comv_CholDiagLower(1:j,j,0)
687 : end do
688 : write(mc_logFileUnit,"(A)")
689 : ProposalErr%occurred = .true.
690 : ProposalErr%msg = PROCEDURE_NAME // &
691 : ": Error occurred while attempting to compute the Cholesky factorization of the &
692 : &covariance matrix of the proposal distribution of " // mc_methodName // "'s sampler. &
693 : &This is highly unusual, and can be indicative of some major underlying problems.\n&
694 : &It may also be due to a runtime computational glitch, in particular, for high-dimensional simulations. &
695 : &In such case, consider increasing the value of the input variable adaptiveUpdatePeriod.\n&
696 : &It may also be that your input objective function has been incorrectly implemented.\n&
697 : &For example, ensure that you are passing a correct value of ndim to the ParaMonte sampler routine,\n&
698 : &the same value that is expected as input to your objective function's implementation.\n&
699 : &Otherwise, restarting the simulation might resolve the error."
700 : call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
701 : return
702 : end if
703 :
704 : end if blockPosDefCheck
705 : ! LCOV_EXCL_STOP
706 :
707 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
708 :
709 : !! perform global adaptive scaling is requested
710 : !if (mc_scalingRequested) then
711 : ! if (meanAccRateSinceStart<mc_targetAcceptanceRate) then
712 : ! mv_adaptiveScaleFactorSq_save = mc_maxScaleFactorSq**((mc_targetAcceptanceRate-meanAccRateSinceStart)/mc_targetAcceptanceRate)
713 : ! else
714 : ! mv_adaptiveScaleFactorSq_save = mc_maxScaleFactorSq**((mc_targetAcceptanceRate-meanAccRateSinceStart)/(1._RK-mc_targetAcceptanceRate))
715 : ! end if
716 : ! mv_adaptiveScaleFactorSq_save = mv_adaptiveScaleFactorSq_save * (meanAccRateSinceStart/mc_targetAcceptanceRate)**mc_ndimInverse
717 : !end if
718 30428 : if (mc_scalingRequested) scalingNeeded = .true.
719 :
720 : else blockSufficientSampleSizeCheck
721 :
722 : ! singularity has occurred. If the first covariance merging has not occurred yet, set the scaling factor appropriately to shrink the covariance matrix.
723 :
724 93988 : samplerUpdateSucceeded = .false.
725 93988 : if (mv_sampleSizeOld_save==1_IK .or. mc_scalingRequested) then
726 84696 : scalingNeeded = .true.
727 84696 : adaptationMeasureComputationNeeded = .true.
728 : ! save the old covariance matrix for the computation of the adaptation measure
729 254014 : do j = 1, nd
730 507954 : do i = 1,j
731 423258 : CovMatUpperOld(i,j) = comv_CholDiagLower(i,j,0)
732 : end do
733 : end do
734 : else
735 9292 : adaptationMeasure = 0._RK
736 9292 : adaptationMeasureComputationNeeded = .false.
737 : end if
738 :
739 :
740 : end if blockSufficientSampleSizeCheck
741 :
742 : ! adjust the scale of the covariance matrix and the Cholesky factor, if needed
743 :
744 : !scalingNeeded = .true.
745 124416 : if (scalingNeeded) then
746 93793 : if ( meanAccRateSinceStart < mc_TargetAcceptanceRateLimit(1) .or. meanAccRateSinceStart > mc_TargetAcceptanceRateLimit(2) ) then
747 9919 : mv_adaptiveScaleFactorSq_save = (meanAccRateSinceStart/mc_targetAcceptanceRate)**mc_ndimInverse
748 : !block
749 : ! use Statistics_mod, only: getRandUniform
750 : ! integer, save :: counter = 0_IK
751 : ! counter = counter - 1
752 : ! mv_adaptiveScaleFactorSq_save = mv_adaptiveScaleFactorSq_save * exp(-counter*getRandUniform(-1.e0_RK,1.e0_RK)/1.e4_RK)
753 : ! !use Statistics_mod, only: getRandInt
754 : ! !mv_adaptiveScaleFactorSq_save = mv_adaptiveScaleFactorSq_save * exp(real(getRandInt(-1_IK,1_IK),kind=RK))
755 : ! !write(*,*) counter, mv_adaptiveScaleFactorSq_save
756 : !end block
757 9919 : adaptiveScaleFactor = sqrt(mv_adaptiveScaleFactorSq_save)
758 29734 : do j = 1, nd
759 : ! update the Cholesky diagonal elements
760 19815 : comv_CholDiagLower(j,0,0) = comv_CholDiagLower(j,0,0) * adaptiveScaleFactor
761 : ! update covariance matrix
762 49526 : do i = 1,j
763 49526 : comv_CholDiagLower(i,j,0) = comv_CholDiagLower(i,j,0) * mv_adaptiveScaleFactorSq_save
764 : end do
765 : ! update the Cholesky factorization
766 39630 : do i = j+1, nd
767 29711 : comv_CholDiagLower(i,j,0) = comv_CholDiagLower(i,j,0) * adaptiveScaleFactor
768 : end do
769 : end do
770 : end if
771 : end if
772 :
773 : ! compute the adaptivity only if any updates has occurred
774 :
775 :
776 124416 : blockAdaptationMeasureComputation: if (adaptationMeasureComputationNeeded) then
777 :
778 333577 : logSqrtDetNew = sum(log( comv_CholDiagLower(1:nd,0,0) ))
779 :
780 : ! use the universal upper bound
781 :
782 333577 : do j = 1, nd
783 655359 : do i = 1,j
784 540235 : CovMatUpperCurrent(i,j) = 0.5_RK * ( comv_CholDiagLower(i,j,0) + CovMatUpperOld(i,j) ) ! dummy
785 : end do
786 : end do
787 115124 : call getLogSqrtDetPosDefMat(nd,CovMatUpperCurrent,logSqrtDetSum,singularityOccurred)
788 115124 : if (singularityOccurred) then
789 : ! LCOV_EXCL_START
790 : write(mc_logFileUnit,"(A)")
791 : write(mc_logFileUnit,"(A)") "Singular covariance matrix detected while computing the Adaptation measure:"
792 : write(mc_logFileUnit,"(A)")
793 : do j = 1, nd
794 : write(mc_logFileUnit,"(*(E25.15))") CovMatUpperCurrent(1:j,j)
795 : end do
796 : ProposalErr%occurred = .true.
797 : ProposalErr%msg = PROCEDURE_NAME // &
798 : ": Error occurred while computing the Cholesky factorization of &
799 : &a matrix needed for the computation of the Adaptation measure. &
800 : &Such error is highly unusual, and requires an in depth investigation of the case.\n&
801 : &It may also be due to a runtime computational glitch, in particular, for high-dimensional simulations. &
802 : &In such case, consider increasing the value of the input variable adaptiveUpdatePeriod.\n&
803 : &It may also be that your input objective function has been incorrectly implemented.\n&
804 : &For example, ensure that you are passing a correct value of ndim to the ParaMonte sampler routine,\n&
805 : &the same value that is expected as input to your objective function's implementation.\n&
806 : &Otherwise, restarting the simulation might resolve the error."
807 : call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
808 : return
809 : ! LCOV_EXCL_STOP
810 : end if
811 :
812 : !adaptationMeasure = 1._RK - exp( 0.5_RK*(mv_logSqrtDetOld_save+logSqrtDetNew) - logSqrtDetSum )
813 115124 : adaptationMeasure = 1._RK - exp( 0.5*(mv_logSqrtDetOld_save + logSqrtDetNew) - logSqrtDetSum ) ! totalVariationUpperBound
814 115124 : if (adaptationMeasure>0._RK) then
815 40324 : adaptationMeasure = sqrt(adaptationMeasure) ! totalVariationUpperBound
816 : ! LCOV_EXCL_START
817 : elseif (adaptationMeasure<0._RK) then
818 : call warn ( prefix = mc_methodBrand &
819 : , outputUnit = mc_logFileUnit &
820 : , msg = mc_negativeTotalVariationMsg//num2str(adaptationMeasure) )
821 : adaptationMeasure = 0._RK
822 : end if
823 : ! LCOV_EXCL_STOP
824 115124 : mv_logSqrtDetOld_save = logSqrtDetNew
825 :
826 : !block
827 : !integer, save :: counter = 0
828 : !counter = counter + 1
829 : !!if (counter==1) then
830 : !if (adaptationMeasure>1._RK) then
831 : !write(*,*)
832 : !write(*,*) mv_logSqrtDetOld_save
833 : !write(*,*) logSqrtDetNew
834 : !write(*,*) logSqrtDetSum
835 : !write(*,*) mv_logSqrtDetOld_save + logSqrtDetNew - 2_IK * logSqrtDetSum
836 : !write(*,*) exp( mv_logSqrtDetOld_save + logSqrtDetNew - 2_IK * logSqrtDetSum )
837 : !write(*,*)
838 : !end if
839 : !end block
840 :
841 : ! update the higher-stage delayed-rejection Cholesky Lower matrices
842 :
843 115124 : if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
844 115124 : call getInvCovMat()
845 :
846 : end if blockAdaptationMeasureComputation
847 :
848 : ! read/write restart file
849 :
850 124416 : if (mc_Image%isLeader) then
851 124416 : if (isFreshRun) then
852 75846 : call writeRestartFile()
853 : else
854 48570 : call readRestartFile()
855 : end if
856 : end if
857 :
858 124416 : end subroutine doAdaptation
859 :
860 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
861 :
862 : ! LCOV_EXCL_START
863 : ! ATTN: This routine needs further correction for the delayed rejection method
864 : #if false
865 : subroutine doAutoTune ( adaptationMeasure &
866 : , AutoTuneScaleSq &
867 : )
868 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
869 : !DEC$ ATTRIBUTES DLLEXPORT :: doAutoTune
870 : #endif
871 : use Matrix_mod, only: getLogSqrtDetPosDefMat
872 : use Constants_mod, only: RK, IK
873 : use Err_mod, only: abort
874 : implicit none
875 : character(*), parameter :: PROCEDURE_NAME = MODULE_NAME // "@doAutoTune()"
876 :
877 : real(RK) , intent(in) :: AutoTuneScaleSq(1)
878 : real(RK) , intent(inout) :: adaptationMeasure
879 : real(RK) :: logSqrtDetSum, logSqrtDetOld, logSqrtDetNew
880 : real(RK) :: CovMatUpperOld(1,1), CovMatUpperCurrent(1,1)
881 : logical :: singularityOccurred
882 :
883 : CovMatUpperOld = comv_CholDiagLower(1:mc_ndim,1:mc_ndim,0)
884 : logSqrtDetOld = sum(log( comv_CholDiagLower(1:mc_ndim,0,0) ))
885 :
886 : if (AutoTuneScaleSq(1)==0._RK) then
887 : comv_CholDiagLower(1,1,0) = 0.25_RK*comv_CholDiagLower(1,1,0)
888 : comv_CholDiagLower(1,0,0) = sqrt(comv_CholDiagLower(1,1,0))
889 : else
890 : comv_CholDiagLower(1,1,0) = AutoTuneScaleSq(1)
891 : comv_CholDiagLower(1,0,0) = sqrt(AutoTuneScaleSq(1))
892 : end if
893 :
894 : ! compute the adaptivity
895 :
896 : logSqrtDetNew = sum(log( comv_CholDiagLower(1:mc_ndim,0,0) ))
897 : CovMatUpperCurrent = 0.5_RK * ( comv_CholDiagLower(1:mc_ndim,1:mc_ndim,0) + CovMatUpperOld )
898 : call getLogSqrtDetPosDefMat(1_IK,CovMatUpperCurrent,logSqrtDetSum,singularityOccurred)
899 : if (singularityOccurred) then
900 : ProposalErr%occurred = .true.
901 : ProposalErr%msg = PROCEDURE_NAME // &
902 : ": Error occurred while computing the Cholesky factorization of &
903 : &a matrix needed for the computation of the proposal distribution's adaptation measure. &
904 : &Such error is highly unusual, and requires an in depth investigation of the case. &
905 : &It may also be that your input objective function has been incorrectly implemented.\n&
906 : &For example, ensure that you are passing a correct value of ndim to the ParaMonte sampler routine,\n&
907 : &the same value that is expected as input to your objective function's implementation.\n&
908 : &Otherwise, restarting the simulation might resolve the error."
909 : call abort( Err = ProposalErr, prefix = mc_methodBrand, newline = "\n", outputUnit = mc_logFileUnit )
910 : return
911 : end if
912 : adaptationMeasure = 1._RK - exp( 0.5_RK*(logSqrtDetOld+logSqrtDetNew) - logSqrtDetSum )
913 :
914 : end subroutine doAutoTune
915 : ! LCOV_EXCL_STOP
916 : #endif
917 :
918 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
919 :
920 : ! Note: based on some benchmarks with ndim = 1, the new design with merging cholesky diag and lower is faster than the original
921 : ! Note: double-communication, implementation. Here are some timings on 4 images:
922 : ! Note: new single communication:
923 : ! Note: image 2: avgTime = 6.960734060198531E-006
924 : ! Note: image 3: avgTime = 7.658279491640721E-006
925 : ! Note: image 4: avgTime = 9.261191328273417E-006
926 : ! Note: avg(avgTime): 7.960068293370891e-06
927 : ! Note: old double communication:
928 : ! Note: image 4: avgTime = 1.733505615104837E-005
929 : ! Note: image 3: avgTime = 1.442608268140151E-005
930 : ! Note: image 2: avgTime = 1.420345299036357E-005
931 : ! Note: avg(avgTime): 1.532153060760448e-05
932 : ! Note: avg(speedup): 1.924798889020109
933 : ! Note: One would expect this speed up to diminish as ndim goes to infinity,
934 : ! Note: since data transfer will dominate the communication overhead.
935 : !> \brief
936 : !> Broadcast adaptation to all images.
937 : !> \warning
938 : !> When CAF parallelism is used, this routine must be exclusively called by the rooter images.
939 : !> When MPI parallelism is used, this routine must be called by all images.
940 : #if defined CAF_ENABLED || MPI_ENABLED
941 67452 : subroutine bcastAdaptation()
942 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
943 : !DEC$ ATTRIBUTES DLLEXPORT :: bcastAdaptation
944 : #endif
945 : #if defined CAF_ENABLED
946 : implicit none
947 67452 : comv_CholDiagLower(1:mc_ndim,0:mc_ndim,0) = comv_CholDiagLower(1:mc_ndim,0:mc_ndim,0)[1]
948 67452 : if (mc_delayedRejectionRequested) call updateDelRejCholDiagLower() ! update the higher-stage delayed-rejection Cholesky Lower matrices
949 : #elif defined MPI_ENABLED
950 : use mpi ! LCOV_EXCL_LINE
951 : implicit none
952 : integer :: ierrMPI
953 : call mpi_bcast ( comv_CholDiagLower & ! LCOV_EXCL_LINE ! buffer: XXX: first element need not be shared. This may need a fix in future.
954 : , mc_ndimSqPlusNdim & ! LCOV_EXCL_LINE ! count
955 : , mpi_double_precision & ! LCOV_EXCL_LINE ! datatype
956 : , 0 & ! LCOV_EXCL_LINE ! root: broadcasting rank
957 : , mpi_comm_world & ! LCOV_EXCL_LINE ! comm
958 : , ierrMPI & ! LCOV_EXCL_LINE ! ierr
959 : )
960 : ! It is essential for the following to be exclusively done by the rooter images. The leaders have had their updates in `doAdaptation()`.
961 : if (mc_Image%isRooter .and. mc_delayedRejectionRequested) call updateDelRejCholDiagLower()
962 : #endif
963 67452 : call getInvCovMat()
964 191868 : end subroutine bcastAdaptation
965 : #endif
966 :
967 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
968 :
969 : !> \brief
970 : !> Return the inverse covariance matrix of the current covariance of the proposal distribution.
971 : !>
972 : !> \warning
973 : !> This procedure is `impure` and sets the values of multiple global variables upon exit.
974 : !>
975 : !> \remark
976 : !> This procedure is required for the proposal probabilities
977 : !>
978 : !> \todo
979 : !> The specification of the array bounds causes the compiler to generate a temporary copy of the array, despite being unnecessary.
980 : !> Right now, this is resolved by replacing the array bounds with `:`. A better solution is to add the `contiguous` attribute to
981 : !> the corresponding argument of [getInvMatFromCholFac](ref matrix_mod::getinvmatfromcholfac) to guarantee it to the compiler.
982 : !> More than improving performance, this would turn off the pesky compiler warnings about temporary array creation.
983 183263 : subroutine getInvCovMat()
984 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
985 : !DEC$ ATTRIBUTES DLLEXPORT :: getInvCovMat
986 : #endif
987 : use Matrix_mod, only: getInvMatFromCholFac ! LCOV_EXCL_LINE
988 : implicit none
989 : integer(IK) :: istage
990 : ! update the inverse covariance matrix of the proposal from the computed Cholesky factor
991 183263 : do concurrent(istage=0:mc_DelayedRejectionCount)
992 : ! WARNING: Do not set the full boundaries' range `(1:mc_ndim)` for the first index of `comv_CholDiagLower` in the following subroutine call.
993 : ! WARNING: Setting the boundaries forces the compiler to generate a temporary array.
994 : mv_InvCovMat(1:mc_ndim,1:mc_ndim,istage) = getInvMatFromCholFac ( nd = mc_ndim & ! LCOV_EXCL_LINE
995 : , CholeskyLower = comv_CholDiagLower(:,1:mc_ndim,istage) & ! LCOV_EXCL_LINE
996 : , Diagonal = comv_CholDiagLower(1:mc_ndim,0,istage) & ! LCOV_EXCL_LINE
997 1941215 : )
998 1045853 : mv_logSqrtDetInvCovMat(istage) = -sum(log( comv_CholDiagLower(1:mc_ndim,0,istage) ))
999 : end do
1000 : !if (mc_ndim==1 .and. abs(log(sqrt(mv_InvCovMat(1,1,0)))-mv_logSqrtDetInvCovMat(0))>1.e-13_RK) then
1001 : !write(*,"(*(g0,:,' '))") "log(sqrt(mv_InvCovMat(1,1,0))) /= mv_logSqrtDetInvCovMat(0)"
1002 : !write(*,"(*(g0,:,' '))") log(sqrt(mv_InvCovMat(1,1,0))), mv_logSqrtDetInvCovMat(0)
1003 : !write(*,"(*(g0,:,' '))") abs(-log(sqrt(mv_InvCovMat(1,1,0)))-mv_logSqrtDetInvCovMat(0))
1004 : !error stop
1005 : !endif
1006 183263 : end subroutine getInvCovMat
1007 :
1008 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1009 :
1010 : !> \brief
1011 : !> Update the Cholesky factorization of the higher-stage delayed rejection covariance matrices.
1012 : !>
1013 : !> \warning
1014 : !> This procedure is impure and sets the values of multiple global variables upon exit.
1015 : !>
1016 : !> \todo
1017 : !> The performance of this update could be improved by only updating the higher-stage covariance, only when needed.
1018 : !> However, the gain will be likely minimal, especially in low-dimensions.
1019 30546 : subroutine updateDelRejCholDiagLower()
1020 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1021 : !DEC$ ATTRIBUTES DLLEXPORT :: updateDelRejCholDiagLower
1022 : #endif
1023 : implicit none
1024 : integer(IK) :: j, istage
1025 : ! update the Cholesky factor of the delayed-rejection-stage proposal distributions
1026 149274 : do istage = 1, mc_DelayedRejectionCount
1027 350142 : comv_CholDiagLower(1:mc_ndim,0,istage) = comv_CholDiagLower(1:mc_ndim,0,istage-1) * mc_DelayedRejectionScaleFactorVec(istage)
1028 380688 : do j = 1, mc_ndim
1029 462828 : comv_CholDiagLower(j+1:mc_ndim,j,istage) = comv_CholDiagLower(j+1:mc_ndim,j,istage-1) * mc_DelayedRejectionScaleFactorVec(istage)
1030 : end do
1031 : end do
1032 : ! There is no need to check for positive-definiteness of the comv_CholDiagLower, it is already checked on the first image.
1033 213809 : end subroutine updateDelRejCholDiagLower
1034 :
1035 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1036 :
1037 : !> \brief
1038 : !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
1039 : !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
1040 : !> This procedure is called by the sampler kernel routines.\n
1041 : !> Write the restart information to the output file.
1042 : !>
1043 : !> @param[in] meanAccRateSinceStart : The current mean acceptance rate of the sampling (**optional**).
1044 152166 : subroutine writeRestartFile(meanAccRateSinceStart)
1045 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1046 : !DEC$ ATTRIBUTES DLLEXPORT :: writeRestartFile
1047 : #endif
1048 : implicit none
1049 : real(RK), intent(in), optional :: meanAccRateSinceStart
1050 : integer(IK) :: i, j
1051 152166 : if (present(meanAccRateSinceStart)) then
1052 76083 : if (mc_isBinaryRestartFileFormat) then
1053 74756 : write(mc_restartFileUnit) meanAccRateSinceStart
1054 : else
1055 1327 : write(mc_restartFileUnit,mc_restartFileFormat) "meanAcceptanceRateSinceStart", meanAccRateSinceStart
1056 : end if
1057 76083 : elseif (mc_isAsciiRestartFileFormat) then
1058 1327 : write( mc_restartFileUnit, mc_restartFileFormat ) "sampleSize" & ! sampleSizeOld
1059 : , mv_sampleSizeOld_save &
1060 1327 : , "logSqrtDeterminant" & ! logSqrtDetOld
1061 : , mv_logSqrtDetOld_save &
1062 1327 : , "adaptiveScaleFactorSquared" & ! adaptiveScaleFactorSq
1063 : , mv_adaptiveScaleFactorSq_save * mc_defaultScaleFactorSq &
1064 1327 : , "meanVec" & ! MeanOld(1:ndim)
1065 : , mv_MeanOld_save(1:mc_ndim) & ! LCOV_EXCL_LINE
1066 1327 : , "covMat" & ! CholDiagLower(1:ndim,0:ndim,0)
1067 8530 : , ((comv_CholDiagLower(i,j,0),i=1,j),j=1,mc_ndim)
1068 : !, (comv_CholDiagLower(1:mc_ndim,0:mc_ndim,0)
1069 : end if
1070 152166 : flush(mc_restartFileUnit)
1071 182712 : end subroutine writeRestartFile
1072 :
1073 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1074 :
1075 : !> \brief
1076 : !> This procedure is a static method of the [ParaDXXXProposalNormal_type](@ref paradxxxproposalnormal_type)
1077 : !> or [ParaDXXXProposalUniform_type](@ref paradxxxproposaluniform_type) classes.\n
1078 : !> This procedure is called by the sampler kernel routines.\n
1079 : !> Read the restart information from the restart file.
1080 : !>
1081 : !> @param[out] meanAccRateSinceStart : The current mean acceptance rate of the sampling (**optional**).
1082 97204 : subroutine readRestartFile(meanAccRateSinceStart)
1083 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
1084 : !DEC$ ATTRIBUTES DLLEXPORT :: readRestartFile
1085 : #endif
1086 : implicit none
1087 : real(RK), intent(out), optional :: meanAccRateSinceStart
1088 : integer(IK) :: i
1089 97204 : if (present(meanAccRateSinceStart)) then
1090 48602 : if (mc_isBinaryRestartFileFormat) then
1091 48065 : read(mc_restartFileUnit) meanAccRateSinceStart
1092 : else
1093 537 : read(mc_restartFileUnit,*)
1094 537 : read(mc_restartFileUnit,*) meanAccRateSinceStart
1095 : end if
1096 48602 : elseif (mc_isAsciiRestartFileFormat) then
1097 7518 : do i = 1, 8 + mc_ndim * (mc_ndim+3) / 2
1098 : !read( mc_restartFileUnit, mc_restartFileFormat )
1099 7518 : read( mc_restartFileUnit, * )
1100 : end do
1101 : end if
1102 249370 : end subroutine readRestartFile
1103 :
1104 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1105 :
1106 :
1107 : #undef ParaDXXX
1108 : #undef GET_RANDOM_PROPOSAL
1109 :
|
