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 : #if CAF_ENABLED
18 : #define SET_CAF(X)X
19 : #else
20 : #define SET_CAF(X)
21 : #endif
22 : ! \bug
23 : ! Avoid Intel ifort bug for too many `use` statements in a submodule by placing them all in the submodule header.
24 : #if CHECK_ENABLED
25 : use pm_err, only: setAsserted, getFine
26 : #define CHECK_ASSERTION(LINE,ASSERTION,MSG)call setAsserted(ASSERTION,getFine(__FILE__,LINE)//MODULE_NAME//MSG);
27 : #else
28 : #define CHECK_ASSERTION(LINE,ASSERTION,MSG) continue;
29 : #endif
30 : ! Define error handling.
31 : #define RETURN_IF_FAILED(LINE,FAILED,MSG) \
32 : if (FAILED) then; \
33 : err%occurred = .true._LK; \
34 : err%msg = PROCEDURE_NAME//getLine(LINE)//SKC_": "//MSG; \
35 : call spec%disp%stop%show(err%msg); \
36 : return; \
37 : end if;
38 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39 : #if ParaDISE_ENABLED || ParaDRAM_ENABLED
40 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41 :
42 : use pm_err, only: getLine
43 : use pm_err, only: err_type
44 : use pm_val2str, only: getStr
45 : use pm_kind, only: SKC => SK, SK, IK, LK
46 : use pm_paramonte, only: PARAMONTE_WEB_ISSUES
47 : #if ParaDISE_ENABLED
48 : use pm_sampling_dise, only: NL1, NL2
49 : use pm_sampling_dise, only: spec_type => specdise_type
50 : use pm_sampling_dise, only: stat_type => statdise_type
51 : use pm_sampling_scio, only: getErrChainRead => getErrChainReadDISE
52 : use pm_sampling_scio, only: getErrChainWrite => getErrChainWriteDISE
53 : use pm_sampling_scio, only: chainFileColName => chainFileColNameDISE
54 : use pm_sampling_scio, only: isFailedChainResize => isFailedChainResizeDISE
55 : implicit none
56 : character(*,SKC), parameter :: MODULE_NAME = SK_"@pm_sampling_proposal_dise"
57 : #elif ParaDRAM_ENABLED
58 : use pm_sampling_dram, only: NL1, NL2
59 : use pm_sampling_dram, only: spec_type => specdram_type
60 : use pm_sampling_dram, only: stat_type => statdram_type
61 : use pm_sampling_scio, only: getErrChainRead => getErrChainReadDRAM
62 : use pm_sampling_scio, only: getErrChainWrite => getErrChainWriteDRAM
63 : use pm_sampling_scio, only: chainFileColName => chainFileColNameDRAM
64 : use pm_sampling_scio, only: isFailedChainResize => isFailedChainResizeDRAM
65 : implicit none
66 : character(*,SKC), parameter :: MODULE_NAME = SK_"@pm_sampling_proposal_dram"
67 : #endif
68 : #if CAF_ENABLED
69 : real(RKC), save , allocatable :: comv_covLowChoUpp(:,:)[:] ! lower covariance, upper Cholesky.
70 : #endif
71 : integer(IK) , private :: ndimp1
72 : integer(IK) , private :: ndim
73 : type :: scaleFactorSq_type
74 : real(RKC) :: default
75 : real(RKC) :: running
76 : end type
77 : type :: proposal_type
78 : ! the following are made components for the sake of thread-safe save attribute and the restart file generation.
79 : character(:,SKC), allocatable :: format
80 : type(scaleFactorSq_type) :: scaleFactorSq
81 : logical(LK) :: delRejEnabled
82 : integer(IK) :: sampleSizeOld
83 : real(RKC) :: logSqrtDetOld
84 : real(RKC) , allocatable :: covLowChoUpp(:,:,:) ! The covariance Matrix of the proposal distribution. Last index belongs to delayed rejection.
85 : real(RKC) , allocatable :: meanOld(:)
86 : #if ParaDISE_ENABLED
87 : real(RKC) :: negLogVolUnitBall
88 : real(RKC) , allocatable :: logSqrtDetInvCov(:) ! (0 : spec%proposalDelayedRejectionCount%val)
89 : real(RKC) , allocatable :: invCov(:,:,:)
90 : #define SET_ParaDISE(X)X
91 : #else
92 : #define SET_ParaDISE(X)
93 : #endif
94 : end type
95 : interface proposal_type
96 : procedure :: constructProposal
97 : end interface
98 :
99 : contains
100 :
101 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
102 :
103 0 : subroutine killMeAlreadyCMake1_RK5(); use pm_sampling_scio_RK5, only: RKC; end subroutine
104 0 : subroutine killMeAlreadyCMake1_RK4(); use pm_sampling_scio_RK4, only: RKC; end subroutine
105 0 : subroutine killMeAlreadyCMake1_RK3(); use pm_sampling_scio_RK3, only: RKC; end subroutine
106 0 : subroutine killMeAlreadyCMake1_RK2(); use pm_sampling_scio_RK2, only: RKC; end subroutine
107 0 : subroutine killMeAlreadyCMake1_RK1(); use pm_sampling_scio_RK1, only: RKC; end subroutine
108 : #if ParaDISE_ENABLED
109 0 : subroutine killMeAlreadyCMake2_RK5(); use pm_sampling_dise_RK5, only: RKC; end subroutine
110 0 : subroutine killMeAlreadyCMake2_RK4(); use pm_sampling_dise_RK4, only: RKC; end subroutine
111 0 : subroutine killMeAlreadyCMake2_RK3(); use pm_sampling_dise_RK3, only: RKC; end subroutine
112 0 : subroutine killMeAlreadyCMake2_RK2(); use pm_sampling_dise_RK2, only: RKC; end subroutine
113 0 : subroutine killMeAlreadyCMake2_RK1(); use pm_sampling_dise_RK1, only: RKC; end subroutine
114 : #elif ParaDRAM_ENABLED
115 0 : subroutine killMeAlreadyCMake2_RK5(); use pm_sampling_dram_RK5, only: RKC; end subroutine
116 0 : subroutine killMeAlreadyCMake2_RK4(); use pm_sampling_dram_RK4, only: RKC; end subroutine
117 0 : subroutine killMeAlreadyCMake2_RK3(); use pm_sampling_dram_RK3, only: RKC; end subroutine
118 0 : subroutine killMeAlreadyCMake2_RK2(); use pm_sampling_dram_RK2, only: RKC; end subroutine
119 0 : subroutine killMeAlreadyCMake2_RK1(); use pm_sampling_dram_RK1, only: RKC; end subroutine
120 : #elif ParaNest_ENABLED
121 : subroutine killMeAlreadyCMake2_RK5(); use pm_sampling_nest_RK5, only: RKC; end subroutine
122 : subroutine killMeAlreadyCMake2_RK4(); use pm_sampling_nest_RK4, only: RKC; end subroutine
123 : subroutine killMeAlreadyCMake2_RK3(); use pm_sampling_nest_RK3, only: RKC; end subroutine
124 : subroutine killMeAlreadyCMake2_RK2(); use pm_sampling_nest_RK2, only: RKC; end subroutine
125 : subroutine killMeAlreadyCMake2_RK1(); use pm_sampling_nest_RK1, only: RKC; end subroutine
126 : #else
127 : #error "Unrecognized interface."
128 : #endif
129 :
130 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 :
132 : #if CAF_ENABLED || MPI_ENABLED
133 : !> \brief
134 : !> Broadcast adaptation to all images.
135 : !>
136 : !> \warning
137 : !> When CAF parallelism is used, this routine must be first called by the leader image, then exclusively called by the rooter images.
138 : !> When MPI parallelism is used, this routine must be called by all images.
139 : !>
140 : !> \devnote
141 : !> based on some benchmarks with ndim = 1, the new design with merging cholesky diag and lower is faster than the original
142 : !> double-communication, implementation. Here are some timings on 4 images:
143 : !> new single communication:
144 : !> image 2: avgTime = 6.960734060198531E-006
145 : !> image 3: avgTime = 7.658279491640721E-006
146 : !> image 4: avgTime = 9.261191328273417E-006
147 : !> avg(avgTime): 7.960068293370891e-06
148 : !> old double communication:
149 : !> image 4: avgTime = 1.733505615104837E-005
150 : !> image 3: avgTime = 1.442608268140151E-005
151 : !> image 2: avgTime = 1.420345299036357E-005
152 : !> avg(avgTime): 1.532153060760448e-05
153 : !> avg(speedup): 1.924798889020109
154 : !> One would expect this speed up to diminish as ndim goes to infinity,
155 : !> since data transfer will dominate the communication overhead.
156 : subroutine bcastProposalAdaptation(spec, proposal)
157 : #if CAF_ENABLED
158 : type(proposal_type), intent(inout) :: proposal
159 : type(spec_type), intent(in) :: spec
160 : if (spec%image%is%leader) then
161 : comv_covLowChoUpp(1 : ndim, 1 : ndimp1) = proposal%covLowChoUpp(1: ndim, 1 : ndimp1, 0)
162 : else
163 : proposal%covLowChoUpp(1: ndim, 1 : ndimp1, 0) = comv_covLowChoUpp(1 : ndim, 1 : ndimp1)[1]
164 : if (proposal%delRejEnabled) call setProposalDelRejChol(spec, proposal) ! update the higher-stage delayed-rejection Cholesky Upper matrices.
165 : SET_ParaDISE(call setProposalInvCov(proposal))
166 : end if
167 : #elif MPI_ENABLED
168 : use mpi !mpi_f08, only: mpi_bcast, mpi_double_precision, mpi_comm_world ! LCOV_EXCL_LINE
169 : use, intrinsic :: iso_c_binding, only: c_sizeof
170 : integer, parameter :: NBRK = c_sizeof(0._RKC) ! Number of bytes in the real kind.
171 : type(spec_type), intent(in) :: spec
172 : type(proposal_type), intent(inout) :: proposal
173 : integer :: ierrMPI
174 : call mpi_bcast ( proposal%covLowChoUpp & ! LCOV_EXCL_LINE
175 : !, size(proposal%covLowChoUpp) & ! LCOV_EXCL_LINE ! count
176 : !, mpi_double_precision & ! LCOV_EXCL_LINE ! datatype
177 : , size(proposal%covLowChoUpp) * NBRK & ! LCOV_EXCL_LINE ! count
178 : , mpi_byte & ! LCOV_EXCL_LINE ! datatype
179 : , 0 & ! LCOV_EXCL_LINE ! root: broadcasting rank
180 : , mpi_comm_world & ! LCOV_EXCL_LINE ! comm
181 : , ierrMPI & ! LCOV_EXCL_LINE ! ierr
182 : )
183 : ! It is essential for the following to be exclusively done by the rooter images. The leaders have had their updates in `setProposalAdapted()`.
184 : if (spec%image%is%rooter .and. proposal%delRejEnabled) call setProposalDelRejChol(spec, proposal)
185 : SET_ParaDISE(call setProposalInvCov(spec, proposal))
186 : #endif
187 : end subroutine bcastProposalAdaptation
188 : #endif
189 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190 :
191 : #if ParaDISE_ENABLED
192 : !> \brief
193 : !> Return the inverse covariance matrix of the current covariance of the proposal distribution.
194 : !>
195 : !> \remark
196 : !> This procedure is required for the proposal probabilities
197 : !>
198 : !> \todo
199 : !> The specification of the array bounds causes the compiler to generate a temporary copy of the array, despite being unnecessary.
200 : !> Right now, this is resolved by replacing the array bounds with `:`. A better solution is to add the `contiguous` attribute to
201 : !> the corresponding argument of [getMatInvFromChoLow](ref pm_matrixInv::getMatInvfromcholow) to guarantee it to the compiler.
202 : !> More than improving performance, this would turn off the pesky compiler warnings about temporary array creation.
203 0 : PURE subroutine setProposalInvCov(spec, proposal)
204 : use pm_matrixInv, only: setMatInv, choUpp
205 : type(proposal_type), intent(inout) :: proposal
206 : type(spec_type), intent(in) :: spec
207 : integer(IK) :: istage
208 : ! update the inverse covariance matrix of the proposal from the computed Cholesky factor.
209 : istage = 0
210 0 : call setMatInv(proposal%invCov(:, 1 : ndim, istage), proposal%covLowChoUpp(:, 2 : ndimp1, istage), choUpp)
211 0 : proposal%logSqrtDetInvCov(istage) = -proposal%logSqrtDetOld
212 0 : do istage = 0, spec%proposalDelayedRejectionCount%val
213 0 : call setMatInv(proposal%invCov(:, 1 : ndim, istage), proposal%covLowChoUpp(:, 2 : ndimp1, istage), choUpp)
214 : !proposal%logSqrtDetInvCov(istage) = .5_RKC * getMatMulTraceLog(proposal%covLowChoUpp(:, 2 : ndimp1, istage))
215 0 : proposal%logSqrtDetInvCov(istage) = proposal%logSqrtDetInvCov(istage) - ndim * log(spec%proposalDelayedRejectionScaleFactor%val(istage))
216 : end do
217 0 : end subroutine setProposalInvCov
218 :
219 : !> \brief
220 : !> Generate and return the natural logarithm of the probability of acceptance of the new state given the old state.
221 0 : PURE function getProposalJumpLogProb(spec, proposal, delayedRejectionStage, origin, destin) result(logPDF)
222 : use pm_distMultiNorm, only: getMultiNormLogPDFNF
223 : use pm_distMultiNorm, only: getMultiNormLogPDF
224 : use pm_ellipsoid, only: isMemberEll
225 : type(spec_type), intent(in) :: spec
226 : type(proposal_type), intent(in) :: proposal
227 : integer(IK), intent(in) :: delayedRejectionStage
228 : real(RKC), intent(in), contiguous :: origin(:), destin(:)
229 : character(*,SKC), parameter :: MSG = SKC_"@getProposalJumpLogProb()/: Internal library error; The specified uniform proposal jump is not within the proposal domain. "
230 : character(*,SKC), parameter :: ERRMSG = MSG//SKC_"ndim, delayedRejectionStage, proposal%invCov(:,1 : ndim, delayedRejectionStage), origin, destin = "
231 : real(RKC) :: logPDF
232 0 : if (spec%proposal%is%normal) then
233 0 : logPDF = getMultiNormLogPDF(destin, origin, proposal%invCov(:,1 : ndim, delayedRejectionStage), logPDFNF = getMultiNormLogPDFNF(ndim, proposal%logSqrtDetInvCov(delayedRejectionStage)))
234 0 : elseif (spec%proposal%is%uniform) then
235 0 : logPDF = proposal%negLogVolUnitBall + proposal%logSqrtDetInvCov(delayedRejectionStage)
236 0 : CHECK_ASSERTION(__LINE__, isMemberEll(proposal%invCov(:,1 : ndim, delayedRejectionStage), origin, destin), \
237 : ERRMSG//getStr([ndim, delayedRejectionStage])//SKC_", "//getStr([proposal%invCov(:,1 : ndim, delayedRejectionStage), origin, destin]))
238 : end if
239 0 : end function getProposalJumpLogProb
240 : #endif
241 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
242 :
243 13 : function constructProposal(spec, err) result(proposal)
244 : use pm_matrixChol, only: setMatChol, lowDia, transHerm
245 : use pm_matrixTrace, only: getMatMulTraceLog
246 : use pm_ellipsoid, only: getLogVolUnitBall ! for uniform proposal.
247 : use pm_arrayRebind, only: setRebound
248 : type(err_type), intent(inout) :: err
249 : type(spec_type), intent(in) :: spec
250 : type(proposal_type) :: proposal
251 : character(*, SK), parameter :: PROCEDURE_NAME = MODULE_NAME//"@constructProposal()"
252 : integer(IK) :: info
253 : !call setNAN(proposal%logSqrtDetOld)
254 13 : ndim = spec%ndim%val
255 13 : ndimp1 = ndim + 1_IK
256 13 : proposal%format = SKC_"("//spec%ndim%str//SKC_"(g0,:,', '),"//new_line(SKC_"a")//SKC_")"
257 66 : proposal%meanOld = spec%proposalStart%val
258 13 : proposal%sampleSizeOld = 1_IK ! This initialization is crucial important for proposal adaptation.
259 13 : proposal%scaleFactorSq%running = 1._RKC
260 13 : proposal%scaleFactorSq%default = spec%proposalScaleFactor%val**2
261 13 : proposal%delRejEnabled = 0_IK < spec%proposalDelayedRejectionCount%val
262 : !proposal%delayedRejection%scaleFactor%val = real(spec%proposalDelayedRejectionScaleFactor%val, RKC)
263 : !proposal%delayedRejection%scaleFactor%log = log(proposal%delayedRejection%scaleFactor%val)
264 : !proposal%domainCubeLimitLower = real(spec%domainCubeLimitLower%val, RKC)
265 : !proposal%domainCubeLimitUpper = real(spec%domainCubeLimitUpper%val, RKC)
266 : !if (spec%domain%isFinite) proposal%domainBallCovMatInv = real(spec%domainBallCovMat%inv, RKC)
267 : !proposal%domainBallCenter = real(spec%domainBallCenter%val, RKC)
268 : ! setup covariance matrix
269 0 : SET_ParaDISE(call setRebound(proposal%logSqrtDetInvCov, 0_IK, spec%proposalDelayedRejectionCount%val))
270 0 : SET_ParaDISE(call setRebound(proposal%invCov, [1_IK, 1_IK, 0_IK], [ndim, ndimp1, spec%proposalDelayedRejectionCount%val]))
271 52 : call setRebound(proposal%covLowChoUpp, [1_IK, 1_IK, 0_IK], [ndim, ndimp1, spec%proposalDelayedRejectionCount%val])
272 : SET_CAF(if (allocated(comv_covLowChoUpp)) deallocate(comv_covLowChoUpp))
273 : SET_CAF(allocate(comv_covLowChoUpp(ndim, 1 : ndimp1)[*]))
274 119 : proposal%covLowChoUpp(1 : ndim, 1 : ndim, 0) = spec%proposalCovMat%val * proposal%scaleFactorSq%default ! Scale the initial covariance matrix
275 13 : call setMatChol(proposal%covLowChoUpp(:, 1 : ndim, 0), lowDia, info, proposal%covLowChoUpp(:, 2 : ndimp1, 0), transHerm) ! The `:` instead of `1 : ndim` avoids temporary array creation.
276 13 : err%occurred = info /= 0_IK
277 13 : if (err%occurred) then
278 0 : err%msg = "Internal erorr: Cholesky factorization of proposalCovMat failed at diagonal #"//getStr(info)//SKC_": "//getStr(transpose(proposal%covLowChoUpp(:, 1 : ndim, 0)), proposal%format)
279 0 : return
280 : end if
281 13 : proposal%logSqrtDetOld = .5_RKC * getMatMulTraceLog(proposal%covLowChoUpp(:, 2 : ndimp1, 0))
282 : ! Scale the higher-stage delayed-rejection Cholesky Upper matrices.
283 13 : if (proposal%delRejEnabled) call setProposalDelRejChol(spec, proposal)
284 0 : SET_ParaDISE(proposal%negLogVolUnitBall = -getLogVolUnitBall(real(ndim, RKC))) ! only for uniform proposal?
285 0 : SET_ParaDISE(call setProposalInvCov(spec, proposal))
286 : ! read/write the first entry of the restart file
287 13 : if (spec%image%is%leader) then
288 : block
289 : real(RKC) :: meanAccRateSinceStart
290 13 : if (spec%run%is%new) then
291 : !if (.not. spec%outputRestartFileFormat%isBinary) then
292 : ! write(spec%restartFile%unit) "ndim"
293 : ! write(spec%restartFile%unit) ndim
294 : !end if
295 13 : call writeRestart(spec, proposal, meanAccRateSinceStart = 1._RKC)
296 : else
297 : !if (.not. spec%outputRestartFileFormat%isBinary) then
298 : ! read(spec%restartFile%unit)
299 : ! read(spec%restartFile%unit)
300 : !end if
301 0 : call readRestart(spec, meanAccRateSinceStart)
302 : end if
303 : end block
304 : end if
305 26 : end function constructProposal
306 :
307 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308 :
309 : !> \brief
310 : !> Update the Cholesky factorization of the higher-stage delayed rejection covariance matrices.
311 : !>
312 : !> \todo
313 : !> The performance of this update could be improved by only updating the higher-stage covariance, only when needed.<br>
314 : !> However, the gain will be likely minimal, especially in low-dimensions.<br>
315 1573 : pure subroutine setProposalDelRejChol(spec, proposal)
316 : type(proposal_type), intent(inout) :: proposal
317 : type(spec_type), intent(in) :: spec
318 : integer(IK) :: idim, ndim, istage
319 1573 : ndim = size(proposal%covLowChoUpp, 1, IK)
320 9438 : do istage = 1, spec%proposalDelayedRejectionCount%val
321 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
322 : !do concurrent(idim = 1 : ndim)
323 32148 : do idim = 1, ndim
324 78600 : proposal%covLowChoUpp(1 : idim, idim + 1, istage) = proposal%covLowChoUpp(1 : idim, idim + 1, istage - 1) * spec%proposalDelayedRejectionScaleFactor%val(istage)
325 : end do
326 : end do
327 : ! There is no need to check for positive-definiteness of the covLowChoUpp, it is already checked on the first image.
328 1573 : end subroutine setProposalDelRejChol
329 :
330 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
331 :
332 0 : subroutine readRestart(spec, meanAccRateSinceStart)
333 : real(RKC), intent(out) :: meanAccRateSinceStart
334 : type(spec_type), intent(in) :: spec
335 : integer(IK) :: idim
336 0 : if (spec%outputRestartFileFormat%isBinary) then
337 0 : read(spec%restartFile%unit) meanAccRateSinceStart
338 : else
339 0 : read(spec%restartFile%unit, *)
340 0 : read(spec%restartFile%unit, *) meanAccRateSinceStart
341 0 : do idim = 1, 8 + ndim * (ndim + 3) / 2; read(spec%restartFile%unit, *); end do
342 : end if
343 0 : end subroutine readRestart
344 :
345 60945 : subroutine writeRestart(spec, proposal, meanAccRateSinceStart)
346 : real(RKC), intent(in) :: meanAccRateSinceStart
347 : type(proposal_type), intent(in) :: proposal
348 : type(spec_type), intent(in) :: spec
349 : integer(IK) :: idim, jdim
350 60945 : if (spec%outputRestartFileFormat%isBinary) then
351 59039 : write(spec%restartFile%unit) meanAccRateSinceStart
352 : else
353 : ! The extra components are all informational for the end user, otherwise not needed for a restart.
354 : write(spec%restartFile%unit, spec%restartFile%format) "meanAcceptanceRateSinceStart" & ! LCOV_EXCL_LINE
355 : , meanAccRateSinceStart & ! LCOV_EXCL_LINE
356 : , "numFuncCall" & ! LCOV_EXCL_LINE
357 : , proposal%sampleSizeOld & ! LCOV_EXCL_LINE
358 : , "proposalAdaptiveScaleFactorSquared" & ! LCOV_EXCL_LINE
359 : , proposal%scaleFactorSq%running * proposal%scaleFactorSq%default & ! LCOV_EXCL_LINE
360 : , "proposalLogVolume" & ! LCOV_EXCL_LINE
361 : , proposal%logSqrtDetOld & ! LCOV_EXCL_LINE
362 : , "proposalMean" & ! LCOV_EXCL_LINE
363 : , proposal%meanOld & ! LCOV_EXCL_LINE
364 : , "proposalCov" & ! LCOV_EXCL_LINE
365 19633 : , ((proposal%covLowChoUpp(idim, jdim, 0), idim = jdim, ndim), jdim = 1, ndim)
366 : ! Only the lower-triangle of the covariance matrix with no delayed rejection.
367 : end if
368 60945 : flush(spec%restartFile%unit)
369 60945 : end subroutine writeRestart
370 :
371 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
372 :
373 : ! Return a newly sampled state to the kernel routine for acceptance check.
374 : ! ATTN: The colon index in place of 1 : ndim below avoids the temporary array creation.
375 : ! The purity breaks only because of the warn() call.
376 1079152 : subroutine setProposalStateNew(spec, proposal, delayedRejectionStage, stateOld, stateNew, rng, err)
377 : use pm_ellipsoid, only: isMemberEll
378 : use pm_distUnifEll, only: setUnifEllRand
379 : use pm_distMultiNorm, only: setMultiNormRand, uppDia
380 : use pm_distUnif, only: xoshiro256ssw_type
381 : type(proposal_type), intent(inout) :: proposal
382 : type(xoshiro256ssw_type), intent(inout) :: rng
383 : real(RKC), intent(in), contiguous :: stateOld(:)
384 : real(RKC), intent(out), contiguous :: stateNew(:)
385 : integer(IK), intent(in) :: delayedRejectionStage
386 : type(spec_type), intent(inout) :: spec
387 : type(err_type), intent(inout) :: err
388 : character(*, SK), parameter :: PROCEDURE_NAME = MODULE_NAME//"@setProposalStateNew()"
389 : character(*, SK), parameter :: WARNING_MSG = PROCEDURE_NAME//SKC_"Number of proposed states out of the objective function domain without any acceptance: "
390 : character(*, SK), parameter :: FATAL_PREFIX = PROCEDURE_NAME//SKC_"Number of proposed states out of the objective function domain without any acceptance has reach the maximum specified value (domainErrCountMax = "
391 : character(*, SK), parameter :: FATAL_SUFFIX = SKC_"). This can occur if the domain of the density function is incorrectly specified or if the density function definition implementation or definition is flawed."
392 : integer(IK) :: domainCheckCounter
393 1079152 : domainCheckCounter = 1_IK
394 12828 : loopBoundaryCheck: do ! Check for the support Region consistency:
395 : ! \todo
396 : ! There is an extremely low but non-zero likelihood that the resulting `stateNew` would be the same as `stateOld`.
397 : ! This could later become problematic because the same state is accepted as new state, potentially corrupting
398 : ! the proposal covariance computation based on the newly-collected samples.
399 1091980 : if (spec%proposal%is%normal) then
400 1091980 : call setMultiNormRand(rng, stateNew, stateOld, proposal%covLowChoUpp(:, 2 : ndimp1, delayedRejectionStage), uppDia)
401 0 : elseif (spec%proposal%is%uniform) then
402 0 : call setUnifEllRand(rng, stateNew, stateOld, proposal%covLowChoUpp(:, 2 : ndimp1, delayedRejectionStage), uppDia)
403 : else
404 0 : error stop "Internal error occurred: proposal distribution unrecognized."
405 : end if
406 : ! check if proposal is in the domain.
407 1091980 : if (spec%domain%isFinite) then
408 387056 : if (spec%domain%isCube) then
409 1984878 : if (all(spec%domainCubeLimitLower%val <= stateNew .and. stateNew <= spec%domainCubeLimitUpper%val)) return
410 0 : elseif (spec%domain%isBall) then
411 0 : if (isMemberEll(spec%domainBallCovMat%inv, spec%domainBallCenter%val, stateNew)) return
412 : end if
413 : else
414 : return
415 : end if
416 12828 : if (domainCheckCounter == spec%domainErrCount%val) call spec%disp%warn%show(WARNING_MSG//getStr(domainCheckCounter))
417 12828 : if (domainCheckCounter == spec%domainErrCountMax%val) exit loopBoundaryCheck
418 12828 : domainCheckCounter = domainCheckCounter + 1_IK
419 : end do loopBoundaryCheck
420 0 : err%msg = FATAL_PREFIX//spec%domainErrCountMax%str//FATAL_SUFFIX
421 0 : err%occurred = .true._LK
422 : end subroutine setProposalStateNew
423 :
424 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
425 :
426 : !> \brief
427 : !> Adapt the proposal distribution based on the specified weighted sampleState.
428 : !>
429 : !> \param[in] sampleState : The two dimensional `(ndim, nsam)` array of accepted states.
430 : !> \param[in] sampleWeight : The one dimensional integer array of the weights of the sampled states in the sampleState.
431 : !> \param[in] adaptationIsGreedy : The logical flag indicating whether the adaptation is in greedy mode..
432 : !> \param[inout] meanAccRateSinceStart : The mean acceptance rate since the start of the sampling.
433 : !> \param[out] adaptationSucceeded : The output logical flag indicating whether the sampling succeeded.
434 : !> \param[out] adaptationMeasure : The output real number in the range `[0,1]` indicating the amount of adaptation,
435 : !> with zero indicating no adaptation and one indicating extreme adaptation to the extent
436 : !> that the new adapted proposal distribution is completely different from the previous proposal.
437 : !> \warning
438 : !> This routine must be exclusively called by the leader images.
439 : !>
440 : !> \note
441 : !> Other than the call to `warn%show()`, this procedure is `pure`.<br>
442 60932 : subroutine setProposalAdapted(spec, proposal, sampleState, sampleWeight, adaptationIsGreedy, meanAccRateSinceStart, adaptationSucceeded, adaptationMeasure, err)
443 : use pm_sampleMean, only: setMean
444 : use pm_matrixTrace, only: getMatMulTraceLog
445 : use pm_sampleCov, only: setCov, setCovMeanMerged
446 : use pm_matrixChol, only: setMatChol, lowDia, transHerm
447 : character(*, SK), parameter :: PROCEDURE_NAME = MODULE_NAME//"@setProposalAdapted()"
448 : type(spec_type), intent(inout) :: spec
449 : type(proposal_type), intent(inout) :: proposal
450 : real(RKC), intent(in), contiguous :: sampleState(:,:)
451 : integer(IK), intent(in), contiguous :: sampleWeight(:)
452 : logical(LK), intent(in) :: adaptationIsGreedy
453 : real(RKC), intent(inout) :: meanAccRateSinceStart ! is intent(out) in restart mode, intent(in) in fresh mode.
454 : logical(LK), intent(out) :: adaptationSucceeded
455 : real(RKC), intent(out) :: adaptationMeasure
456 : type(err_type), intent(inout) :: err
457 121864 : real(RKC) :: sampleStateMeanNew(size(sampleState, 1, IK))
458 121864 : real(RKC) :: sampleStateMeanCurrent(size(sampleState, 1, IK))
459 121864 : real(RKC) :: sampleStateCovCholOld(size(sampleState, 1, IK), size(sampleState, 1, IK) + 1)
460 60932 : real(RKC) :: sampleStateCovLowCurrent(size(sampleState, 1, IK), size(sampleState, 1, IK) + 1)
461 : real(RKC) :: logSqrtDetNew, logSqrtDetSum, adaptiveScaleFactor, fracA
462 : logical(LK) :: adaptationMeasureComputationNeeded ! only used to avoid redundant affinity computation, if no update occurs.
463 : logical(LK) :: scalingNeeded!, singularityOccurred
464 : integer(IK) :: sampleSizeOld, sampleSizeCurrent
465 : integer(IK):: idim, jdim, ndim, nsam, info
466 : integer(IK), parameter :: dim = 2_IK
467 : scalingNeeded = .false._LK
468 : ndim = size(sampleState, 1, IK)
469 60932 : nsam = size(sampleState, 2, IK)
470 60932 : sampleSizeOld = proposal%sampleSizeOld ! this is kept only for restoration of proposal%sampleSizeOld, if needed.
471 60932 : if (spec%image%is%leader .and. spec%run%is%dry) call readRestart(spec, meanAccRateSinceStart)
472 : ! First if there are less than ndimp1 points for new covariance computation, then just scale the covariance and return.
473 60932 : blockSufficientSampleSizeCheck: if (ndim < nsam) then
474 : ! get the upper covariance matrix and mean of the new sample
475 50825 : if (adaptationIsGreedy) then
476 0 : sampleSizeCurrent = nsam
477 0 : call setMean(sampleStateMeanCurrent, sampleState, dim)
478 0 : call setCov(sampleStateCovLowCurrent(:, 1 : ndim), lowDia, sampleStateMeanCurrent, sampleState, dim)
479 : else
480 50825 : call setMean(sampleStateMeanCurrent, sampleState, dim, sampleWeight, sampleSizeCurrent)
481 50825 : call setCov(sampleStateCovLowCurrent(:, 1 : ndim), lowDia, sampleStateMeanCurrent, sampleState, dim, sampleWeight, sampleSizeCurrent)
482 : end if
483 : ! Combine old and new covariance matrices if both exist.
484 1019148 : sampleStateCovCholOld = proposal%covLowChoUpp(:, :, 0) ! This will be used to recover the old covariance in case of update failure, and to compute the adaptation measure.
485 50825 : blockMergeCovMat: if (proposal%sampleSizeOld == 1_IK) then
486 : ! There is no prior old Covariance matrix to combine with the new one from the new sampleState
487 40 : proposal%meanOld(1 : ndim) = sampleStateMeanCurrent
488 13 : proposal%sampleSizeOld = sampleSizeCurrent
489 : ! Copy and then scale the new covariance matrix by the default scale factor, which will be then used to get the Cholesky Factor.
490 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
491 : !do concurrent(jdim = 1 : ndim)
492 40 : do jdim = 1, ndim
493 93 : proposal%covLowChoUpp(jdim : ndim, jdim, 0) = sampleStateCovLowCurrent(jdim : ndim, jdim) * proposal%scaleFactorSq%default
494 : end do
495 : else blockMergeCovMat
496 : ! First scale the new covariance matrix by the default scale factor, which will be then used to get the Cholesky Factor.
497 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
498 : !do concurrent(jdim = 1 : ndim)
499 211153 : do jdim = 1, ndim
500 589665 : sampleStateCovLowCurrent(jdim : ndim, jdim) = sampleStateCovLowCurrent(jdim : ndim, jdim) * proposal%scaleFactorSq%default
501 : end do
502 : ! Now combine it with the old covariance matrix.
503 : ! Do not set the full boundaries' range `(1 : ndim)` for `proposal%covLowChoUpp` in the following subroutine call.
504 : ! Setting the boundaries forces the compiler to generate a temporary array.
505 50812 : fracA = real(proposal%sampleSizeOld, RKC) / real(proposal%sampleSizeOld + sampleSizeCurrent, RKC)
506 50812 : call setCovMeanMerged(proposal%covLowChoUpp(:, 1 : ndim, 0), sampleStateMeanNew, sampleStateCovLowCurrent(:, 1 : ndim), sampleStateMeanCurrent, sampleStateCovCholOld(:, 1 : ndim), proposal%meanOld, fracA, lowDia)
507 211153 : proposal%meanOld(1 : ndim) = sampleStateMeanNew
508 : end if blockMergeCovMat
509 : ! Now get the Cholesky factorization.
510 : ! WARNING: Do not set the full boundaries' range `(1 : ndim)` for the first index of `proposal%covLowChoUpp` in the following subroutine call.
511 : ! WARNING: Setting the boundaries forces the compiler to generate a temporary array.
512 50825 : call setMatChol(proposal%covLowChoUpp(:, 1 : ndim, 0), lowDia, info, proposal%covLowChoUpp(:, 2 : ndimp1, 0), transHerm) ! The `:` instead of `1 : ndim` avoids temporary array creation.
513 50825 : blockPosDefCheck: if (info == 0_IK) then
514 : !singularityOccurred = .false._LK
515 50825 : adaptationSucceeded = .true._LK
516 : adaptationMeasureComputationNeeded = .true._LK
517 50825 : proposal%sampleSizeOld = proposal%sampleSizeOld + sampleSizeCurrent
518 : else blockPosDefCheck
519 0 : adaptationMeasure = 0._RKC
520 : !singularityOccurred = .true._LK
521 0 : adaptationSucceeded = .false._LK
522 : adaptationMeasureComputationNeeded = .false._LK
523 0 : proposal%sampleSizeOld = sampleSizeOld
524 : ! It may be a good idea to add a warning message printed out here for the singularity occurrence.
525 0 : call spec%disp%warn%show(SKC_"Singularity occurred while adapting the proposal distribution covariance matrix.", tmsize = 0_IK, bmsize = 0_IK)
526 : ! Recover the old covariance matrix + Cholesky factor.
527 0 : proposal%covLowChoUpp(:, :, 0) = sampleStateCovCholOld
528 : end if blockPosDefCheck
529 : !! perform global adaptive scaling is requested
530 : !if (spec%targetAcceptanceRate%enabled) then
531 : ! if (meanAccRateSinceStart < spec%targetAcceptanceRate%aim) then
532 : ! proposal%scaleFactorSq%running = mc_maxScaleFactorSq**((spec%targetAcceptanceRate%aim-meanAccRateSinceStart)/spec%targetAcceptanceRate%aim)
533 : ! else
534 : ! proposal%scaleFactorSq%running = mc_maxScaleFactorSq**((spec%targetAcceptanceRate%aim-meanAccRateSinceStart)/(1._RKC-spec%targetAcceptanceRate%aim))
535 : ! end if
536 : ! proposal%scaleFactorSq%running = proposal%scaleFactorSq%running * (meanAccRateSinceStart/spec%targetAcceptanceRate%aim)**spec%ndim%inv
537 : !end if
538 50825 : if (spec%targetAcceptanceRate%enabled) scalingNeeded = .true._LK
539 : else blockSufficientSampleSizeCheck
540 : ! Singularity has occurred. If the first covariance merging has not occurred yet, set the scaling factor appropriately to shrink the covariance matrix.
541 10107 : adaptationSucceeded = .false._LK
542 10107 : if (proposal%sampleSizeOld == 1_IK .or. spec%targetAcceptanceRate%enabled) then
543 : scalingNeeded = .true._LK
544 : adaptationMeasureComputationNeeded = .true._LK
545 : ! Save the old covariance matrix for the computation of the adaptation measure.
546 : ! \todo Could these copies be somehow removed and optimized away? perhaps unimportant if sample size is sufficient most of the times.
547 159 : do jdim = 1, ndim
548 318 : do idim = jdim, ndim
549 265 : sampleStateCovCholOld(idim, jdim) = proposal%covLowChoUpp(idim, jdim, 0)
550 : end do
551 : end do
552 : else
553 10054 : adaptationMeasure = 0._RKC
554 : adaptationMeasureComputationNeeded = .false._LK
555 : end if
556 : end if blockSufficientSampleSizeCheck
557 : ! Adjust the scale of the covariance matrix and the Cholesky factor, if needed.
558 : !scalingNeeded = .true._LK
559 : if (scalingNeeded) then
560 53 : if (meanAccRateSinceStart < spec%targetAcceptanceRate%val(1) .or. spec%targetAcceptanceRate%val(2) < meanAccRateSinceStart) then
561 0 : proposal%scaleFactorSq%running = (meanAccRateSinceStart / spec%targetAcceptanceRate%aim)**spec%ndim%inv
562 : !block
563 : ! use pm_distUnif, only: getUnifRand
564 : ! integer, save :: counter = 0_IK
565 : ! counter = counter - 1
566 : ! proposal%scaleFactorSq%running = proposal%scaleFactorSq%running * exp(-counter*getUnifRand(-1._RKC, 1._RKC)/1.e4_RK)
567 : ! !use pm_statistics, only: getUnifRand
568 : ! !proposal%scaleFactorSq%running = proposal%scaleFactorSq%running * exp(real(getUnifRand(-1_IK, 1_IK), RK))
569 : ! !write(*,*) counter, proposal%scaleFactorSq%running
570 : !end block
571 0 : adaptiveScaleFactor = sqrt(proposal%scaleFactorSq%running)
572 0 : proposal%covLowChoUpp(1 : ndim, 1, 0) = proposal%covLowChoUpp(1 : ndim, 1, 0) * proposal%scaleFactorSq%running ! Update first column of covariance matrix.
573 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
574 : !do concurrent(jdim = 2 : ndim)
575 0 : do jdim = 2, ndim
576 0 : proposal%covLowChoUpp(1 : jdim - 1, jdim, 0) = proposal%covLowChoUpp(1 : jdim, jdim + 1, 0) * adaptiveScaleFactor ! Update the Cholesky factor.
577 0 : proposal%covLowChoUpp(jdim : ndim, jdim, 0) = proposal%covLowChoUpp(jdim : ndim, jdim, 0) * proposal%scaleFactorSq%running ! Update covariance matrix.
578 : end do
579 0 : proposal%covLowChoUpp(1 : ndim, ndimp1, 0) = proposal%covLowChoUpp(1 : ndim, ndimp1, 0) * adaptiveScaleFactor ! Update last column of Cholesky factor.
580 : end if
581 : end if
582 : ! Compute the adaptivity only if any updates has occurred.
583 60932 : blockAdaptationMeasureComputation: if (adaptationMeasureComputationNeeded) then
584 50878 : logSqrtDetNew = getMatMulTraceLog(proposal%covLowChoUpp(:, 2 : ndimp1, 0))
585 : ! Use the universal upper bound.
586 : !block
587 : !use pm_io
588 : !call disp%show(sampleStateCovLowCurrent)
589 : !end block
590 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
591 : !do concurrent(jdim = 1 : ndim)
592 211352 : do jdim = 1, ndim
593 590076 : sampleStateCovLowCurrent(jdim : ndim, jdim) = 0.5_RKC * (proposal%covLowChoUpp(jdim : ndim, jdim, 0) + sampleStateCovCholOld(jdim : ndim, jdim)) ! dummy
594 : end do
595 50878 : call setMatChol(sampleStateCovLowCurrent(:, 1 : ndim), lowDia, info, sampleStateCovLowCurrent(:, 2 : ndimp1), transHerm)
596 50878 : if (info == 0_IK) then !singularityOccurred = info /= 0_IK
597 50878 : logSqrtDetSum = getMatMulTraceLog(sampleStateCovLowCurrent(:, 2 : ndimp1))
598 : else ! singularityOccurred ! LCOV_EXCL_START
599 : err%occurred = .true._LK
600 : err%msg = NL1//PROCEDURE_NAME//SKC_"@line::"//getStr(__LINE__)//SKC_": "//&
601 : SKC_"Singular lower covariance matrix detected at column ("//getStr(info)//SKC_") while computing the adaptationMeasure measure:"//NL2
602 : do info = 1, size(sampleStateCovLowCurrent, 1, IK)
603 : err%msg = err%msg//getStr(sampleStateCovLowCurrent(info,:))//NL1
604 : end do
605 : err%msg = err%msg//NL1//&
606 : SKC_": Error occurred while computing the Cholesky factorization of a matrix needed for the computation of the adaptation measure. &
607 : &Such error is highly unusual, and requires an in depth investigation of the case. &
608 : &It may also be due to a runtime computational glitch, in particular, for high-dimensional simulations. &
609 : &In such case, consider increasing the value of the input variable proposalAdaptationPeriod. &
610 : &It may also be that your input objective function has been incorrectly implemented. &
611 : &Also, ensure a correct value of `ndim` to the ParaMonte sampler routine, representing the domain size. &
612 : &Otherwise, restarting the simulation might resolve the error. If the error persists, please inform the developers at: "//NL2//PARAMONTE_WEB_ISSUES
613 : return
614 : end if ! LCOV_EXCL_STOP
615 50878 : adaptationMeasure = 1._RKC - exp(0.5_RKC * (proposal%logSqrtDetOld + logSqrtDetNew) - logSqrtDetSum) ! totalVariationUpperBound
616 50878 : if (0._RKC < adaptationMeasure) then
617 50813 : adaptationMeasure = sqrt(adaptationMeasure) ! totalVariationUpperBound
618 65 : elseif (adaptationMeasure < 0._RKC) then
619 : spec%msg = PROCEDURE_NAME//SKC_": Non-positive adaptation measure detected" ! LCOV_EXCL_LINE
620 : if (abs(adaptationMeasure) < sqrt(epsilon(0._RKC))) spec%msg = spec%msg//SKC_" (possibly due to round-off error)" ! LCOV_EXCL_LINE
621 : call spec%disp%warn%show(spec%msg//SKC_": "//getStr(adaptationMeasure)) ! LCOV_EXCL_LINE
622 : adaptationMeasure = 0._RKC ! LCOV_EXCL_LINE
623 : end if
624 50878 : proposal%logSqrtDetOld = logSqrtDetNew
625 : !block
626 : !integer, save :: counter = 0
627 : !counter = counter + 1
628 : !!if (counter==1) then
629 : !if (adaptationMeasure>1._RKC) then
630 : !write(*,*)
631 : !write(*,*) proposal%logSqrtDetOld
632 : !write(*,*) logSqrtDetNew
633 : !write(*,*) logSqrtDetSum
634 : !write(*,*) proposal%logSqrtDetOld + logSqrtDetNew - 2_IK * logSqrtDetSum
635 : !write(*,*) exp( proposal%logSqrtDetOld + logSqrtDetNew - 2_IK * logSqrtDetSum )
636 : !write(*,*)
637 : !end if
638 : !end block
639 : ! update the higher-stage delayed-rejection Cholesky Upper matrices
640 50878 : if (proposal%delRejEnabled) call setProposalDelRejChol(spec, proposal)
641 0 : SET_ParaDISE(call setProposalInvCov(spec, proposal))
642 : end if blockAdaptationMeasureComputation
643 60932 : if (spec%image%is%leader .and. spec%run%is%new) call writeRestart(spec, proposal, meanAccRateSinceStart)
644 : end subroutine setProposalAdapted
645 :
646 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
647 :
648 : #undef RETURN_IF_FAILED
649 : #undef SET_ParaDISE
650 :
651 : #else
652 : !%%%%%%%%%%%%%%%%%%%%%%%%
653 : #error "Unrecognized interface."
654 : !%%%%%%%%%%%%%%%%%%%%%%%%
655 : #endif
656 : #undef SET_CAF_MPI
657 : #undef SET_CAF
658 : #undef SHARED
|
