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 CHECK_ENABLED
18 : use pm_err, only: setAsserted
19 : #define CHECK_ASSERTION(LINE,ASSERTION,MSG) \
20 : call setAsserted(ASSERTION,getFine(__FILE__,LINE)//MODULE_NAME//MSG);
21 : #else
22 : #define CHECK_ASSERTION(LINE,ASSERTION,MSG) continue;
23 : #endif
24 : use pm_kind, only: SKC => SK, SK, IK, LK
25 : use pm_matrixClass, only: isMatClass, posdefmat
26 : use pm_sampling_base, only: specbase_type, astatbase_type, sfcbase_type, NL2, NL1
27 : use pm_arrayRemove, only: getRemoved
28 : use pm_arrayResize, only: setResized
29 : use pm_strASCII, only: getStrLower
30 : use pm_val2str, only: getStr
31 : use pm_except, only: setNAN
32 : use pm_except, only: isNAN
33 : use pm_strASCII, only: SUB
34 : use pm_err, only: err_type
35 : use pm_err, only: getFine
36 :
37 : implicit none
38 :
39 : character(*,SKC) , parameter :: MODULE_NAME = SK_"@pm_sampling_mcmc"
40 :
41 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 : ! simulation declarations.
43 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44 :
45 : type :: burninLoc_type
46 : integer(IK) :: compact
47 : integer(IK) :: verbose
48 : end type
49 :
50 : ! type for chain refinements
51 : type, extends(sfcbase_type) :: sfcmcmc_type
52 : real(RKC) , allocatable :: act(:,:) !< \public The array of shape `(1:ndim, 0:nref)` containing he Integrated AutoCorrelation Time at each refinement stage.
53 : contains
54 : procedure , pass :: getErrRefinement
55 : end type
56 :
57 : type, abstract, extends(astatbase_type) :: astatmcmc_type
58 : type(burninLoc_type) :: burninLocMCMC
59 : type(sfcmcmc_type) :: sfc
60 : end type
61 :
62 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 : ! specification declarations.
64 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65 :
66 : !integer(IK) :: outputChainSize ! namelist input
67 : type :: outputChainSize_type
68 : integer(IK) :: val
69 : integer(IK) :: def
70 : integer(IK) :: null
71 : character(:,SKC) , allocatable :: desc
72 : end type
73 :
74 : !integer(IK) :: outputSampleRefinementCount ! namelist input
75 : type :: outputSampleRefinementCount_type
76 : character(:,SKC) , allocatable :: str
77 : integer(IK) :: val
78 : integer(IK) :: def
79 : integer(IK) :: null
80 : character(:,SKC) , allocatable :: desc
81 : end type
82 :
83 : !character(63,SKC) :: outputSampleRefinementMethod ! namelist input
84 : type :: outputSampleRefinementMethod_type
85 : logical(LK) :: isBatchMeansCompact = .false._LK
86 : logical(LK) :: isBatchMeansVerbose = .false._LK
87 : logical(LK) :: isMaxCumSumACF = .false._LK
88 : logical(LK) :: isBatchMeans = .false._LK
89 : character(12,SKC) :: maxCumSumACF = "maxCumSumACF"
90 : character(10,SKC) :: batchMeans = "BatchMeans"
91 : character(9,SKC) :: cutoffACF = "cutoffACF"
92 : character(:,SKC) , allocatable :: def
93 : character(:,SKC) , allocatable :: val
94 : character(:,SKC) , allocatable :: null
95 : character(:,SKC) , allocatable :: desc
96 : end type
97 :
98 : !character(63,SKC) :: proposal ! namelist input
99 : type :: proposalis_type
100 : logical(LK) :: normal = .false._LK
101 : logical(LK) :: uniform = .false._LK
102 : end type
103 : type :: proposal_type
104 : type(proposalis_type) :: is
105 : character(7,SKC) :: uniform
106 : character(6,SKC) :: normal
107 : character(:,SKC) , allocatable :: val
108 : character(:,SKC) , allocatable :: def
109 : character(:,SKC) , allocatable :: null
110 : character(:,SKC) , allocatable :: desc
111 : end type
112 :
113 : !real(RKC) , allocatable :: proposalCorMat(:,:) ! namelist input
114 : type :: proposalCorMat_type
115 : real(RKC) , allocatable :: val(:,:)
116 : real(RKC) , allocatable :: def(:,:)
117 : character(:,SKC) , allocatable :: desc
118 : end type
119 :
120 : !real(RKC) , allocatable :: proposalCovMat(:,:) ! namelist input
121 : type :: proposalCovMat_type
122 : logical(LK) :: isUserSet
123 : real(RKC) , allocatable :: def(:,:)
124 : real(RKC) , allocatable :: val(:,:)
125 : character(:,SKC) , allocatable :: desc
126 : end type
127 :
128 : !character(127,SKC) :: proposalScaleFactor
129 : type :: proposalScaleFactor_type
130 : real(RKC) :: val, valdef
131 : character(:,SKC) , allocatable :: str, strdef, null, desc
132 : end type
133 :
134 : !real(RKC) , allocatable :: proposalStart(:) ! namelist input
135 : type :: proposalStart_type
136 : real(RKC) , allocatable :: val(:)
137 : character(:,SKC) , allocatable :: desc
138 : end type
139 :
140 : !real(RKC) , allocatable :: proposalStartDomainCubeLimitLower(:) ! namelist input
141 : type :: proposalStartDomainCubeLimitLower_type
142 : real(RKC) :: null
143 : real(RKC) :: def
144 : real(RKC) , allocatable :: val(:)
145 : character(:,SKC) , allocatable :: desc
146 : end type
147 :
148 : !real(RKC) , allocatable :: proposalStartDomainCubeLimitUpper(:) ! namelist input
149 : type :: proposalStartDomainCubeLimitUpper_type
150 : real(RKC) :: def
151 : real(RKC) , allocatable :: val(:)
152 : character(:,SKC) , allocatable :: desc
153 : end type
154 :
155 : !logical(LK) :: proposalStartRandomized ! namelist input
156 : type :: proposalStartRandomized_type
157 : logical(LK) :: val
158 : logical(LK) :: def
159 : character(:,SKC) , allocatable :: desc
160 : end type
161 :
162 : !real(RKC) , allocatable :: proposalStdVec(:) ! namelist input
163 : type :: proposalStdVec_type
164 : real(RKC) , allocatable :: val(:)
165 : real(RKC) , allocatable :: def(:)
166 : character(:,SKC) , allocatable :: desc
167 : end type
168 :
169 : type, extends(specbase_type) :: specmcmc_type
170 : type(outputChainSize_type) :: outputChainSize
171 : type(proposal_type) :: proposal
172 : type(ProposalScaleFactor_type) :: ProposalScaleFactor
173 : type(proposalStart_type) :: proposalStart
174 : type(proposalStdVec_type) :: proposalStdVec
175 : type(proposalCorMat_type) :: proposalCorMat
176 : type(proposalCovMat_type) :: proposalCovMat
177 : type(outputSampleRefinementCount_type) :: outputSampleRefinementCount
178 : type(outputSampleRefinementMethod_type) :: outputSampleRefinementMethod
179 : type(proposalStartRandomized_type) :: proposalStartRandomized
180 : type(proposalStartDomainCubeLimitLower_type) :: proposalStartDomainCubeLimitLower
181 : type(proposalStartDomainCubeLimitUpper_type) :: proposalStartDomainCubeLimitUpper
182 : contains
183 : procedure, pass, private :: sanitize
184 : procedure, pass, private :: report
185 : procedure, pass, public :: set
186 : end type
187 :
188 : interface specmcmc_type
189 : module procedure :: construct
190 : end interface
191 :
192 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
193 :
194 : type, private :: method_type
195 : logical(LK) :: isMax = .false._LK
196 : logical(LK) :: isMed = .false._LK
197 : logical(LK) :: isMin = .false._LK
198 : logical(LK) :: isAvg = .false._LK
199 : logical(LK) :: isCumSum = .false._LK
200 : logical(LK) :: isCumSumMax = .false._LK
201 : logical(LK) :: isBatchMeans = .false._LK
202 : logical(LK) :: isBatchMeansMax = .false._LK
203 : logical(LK) :: isViaCompactChain = .true._LK
204 : logical(LK) :: isViaVerboseChain = .true._LK
205 : end type
206 :
207 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
208 :
209 : contains
210 :
211 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
212 :
213 0 : subroutine killMeAlreadyCMake1_RK5(); use pm_sampling_scio_RK5, only: RKC; end subroutine
214 0 : subroutine killMeAlreadyCMake1_RK4(); use pm_sampling_scio_RK4, only: RKC; end subroutine
215 0 : subroutine killMeAlreadyCMake1_RK3(); use pm_sampling_scio_RK3, only: RKC; end subroutine
216 0 : subroutine killMeAlreadyCMake1_RK2(); use pm_sampling_scio_RK2, only: RKC; end subroutine
217 0 : subroutine killMeAlreadyCMake1_RK1(); use pm_sampling_scio_RK1, only: RKC; end subroutine
218 :
219 0 : subroutine killMeAlreadyCMake2_RK5(); use pm_sampling_base_RK5, only: RKC; end subroutine
220 0 : subroutine killMeAlreadyCMake2_RK4(); use pm_sampling_base_RK4, only: RKC; end subroutine
221 0 : subroutine killMeAlreadyCMake2_RK3(); use pm_sampling_base_RK3, only: RKC; end subroutine
222 0 : subroutine killMeAlreadyCMake2_RK2(); use pm_sampling_base_RK2, only: RKC; end subroutine
223 0 : subroutine killMeAlreadyCMake2_RK1(); use pm_sampling_base_RK1, only: RKC; end subroutine
224 :
225 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
226 :
227 14 : function construct(modelr, method, ndim) result(spec)
228 : #if __INTEL_COMPILER && DLL_ENABLED && (_WIN32 || _WIN64)
229 : !DEC$ ATTRIBUTES DLLEXPORT :: construct
230 : #endif
231 : use pm_kind, only: modelr_type
232 : type(modelr_type), intent(in) :: modelr
233 : character(*,SKC), intent(in) :: method
234 : integer(IK), intent(in) :: ndim
235 : type(specmcmc_type) :: spec
236 :
237 14 : spec%specbase_type = specbase_type(modelr, method, ndim)
238 :
239 : outputChainSize_block: block
240 : use pm_sampling_scio, only: outputChainSize
241 14 : spec%outputChainSize%null = -huge(0_IK)
242 14 : spec%outputChainSize%def = 100000_IK
243 : spec%outputChainSize%desc = &
244 : SKC_"The simulation specification `outputChainSize` is a positive scalar of type `integer` whose value determines the number &
245 : &of non-refined, potentially auto-correlated, but unique samples drawn by the MCMC sampler before stopping the sampling. &
246 : &For example, if `outputChainSize = 10000`, then `10000` unique sample points (with no duplicates) will be drawn from the &
247 : &target objective function that the user has provided. &
248 : &The input value for `outputChainSize` must be a positive integer of a minimum value `ndim + 1` or larger, &
249 : &where `ndim` is the number of dimensions of the domain of the objective function to be sampled. &
250 : &Note that `outputChainSize` is different from and always smaller than the length of the constructed MCMC chain. &
251 14 : &The default value for `outputChainSize` is `"//getStr(spec%outputChainSize%def)//SKC_"`."
252 : !!$omp master
253 14 : outputChainSize = spec%outputChainSize%null
254 : !!$omp end master
255 : end block outputChainSize_block
256 :
257 : outputSampleRefinementCount_block: block
258 : use pm_sampling_scio, only: outputSampleRefinementCount
259 14 : spec%outputSampleRefinementCount%null = -huge(0_IK)
260 14 : spec%outputSampleRefinementCount%def = huge(0_IK)
261 : spec%outputSampleRefinementCount%desc = &
262 : SKC_"The simulation specification `outputSampleRefinementCount` is a positive-valued scalar of type `integer`. &
263 : &When `outputSampleSize < 0`, the value of `outputSampleRefinementCount` dictates the maximum number of times &
264 : &the MCMC chain will be refined to remove the autocorrelation within the output MCMC sample. For example,"//NL2//&
265 : SKC_"+ if `outputSampleRefinementCount = 0`,"//NL2//&
266 : SKC_" no refinement of the output MCMC chain will be performed. The resulting MCMC sample will correspond &
267 : &to the full MCMC chain in verbose format (i.e., each sampled state has a weight of one)."//NL2//&
268 : SKC_"+ if `outputSampleRefinementCount = 1`,"//NL2//&
269 : SKC_" the refinement of the output MCMC chain will be done only once if needed, and no more, &
270 : &even though some residual autocorrelation in the output MCMC sample may still exist. &
271 : &In practice, only one refinement of the final output MCMC chain should be enough to remove &
272 : &the existing autocorrelations in the final output sample. Exceptions occur when the integrated &
273 : &Autocorrelation (ACT) of the output MCMC chain is comparable to or larger than the length of the chain. &
274 : &In such cases, neither the BatchMeans method nor any other method of ACT computation will be able to &
275 : &accurately compute the ACT. Consequently, the samples generated based on the computed ACT values will &
276 : &likely not be i.i.d. and will still be significantly autocorrelated. In such scenarios, more than &
277 : &one refinement of the MCMC chain will be necessary. Very small sample size resulting from multiple &
278 : &refinements of the sample could be a strong indication of the bad mixing of the MCMC chain and &
279 : &the lack of convergence to the target objective function."//NL2//&
280 : SKC_"+ if `outputSampleRefinementCount > 1`,"//NL2//&
281 : SKC_" the refinement of the output MCMC chain will be done for a maximum `outputSampleRefinementCount` number of &
282 : ×, even though some residual autocorrelation in the final output MCMC sample may still exist."//NL2//&
283 : SKC_"+ if `outputSampleRefinementCount >> 1` (e.g., comparable to or larger than the length of the MCMC chain),"//NL2//&
284 : SKC_" the refinement of the output MCMC chain will continue until the integrated autocorrelation of the resulting &
285 : &final sample is less than 2, virtually implying that an independent identically-distributed (i.i.d.) sample &
286 : &from the target objective function has finally been obtained."//NL2//&
287 : SKC_"Note that to obtain i.i.d. samples from a multidimensional chain, the sampler will, by default, use the maximum of &
288 : &integrated Autocorrelation (ACT) among all chain dimensions to refine the chain. &
289 : &Note that the value specified for `outputSampleRefinementCount` is used only when the variable `outputSampleSize < 0`, &
290 14 : &otherwise, it will be ignored. The default value for `outputSampleRefinementCount` is `"//getStr(spec%outputSampleRefinementCount%def)//SKC_"`."
291 : !!$omp master
292 14 : outputSampleRefinementCount = spec%outputSampleRefinementCount%null
293 : !!$omp end master
294 : end block outputSampleRefinementCount_block
295 :
296 : outputSampleRefinementMethod_block: block
297 : use pm_sampling_scio, only: outputSampleRefinementMethod
298 14 : spec%outputSampleRefinementMethod%def = spec%outputSampleRefinementMethod%batchMeans
299 14 : spec%outputSampleRefinementMethod%null = repeat(SUB, len(outputSampleRefinementMethod, IK))
300 : spec%outputSampleRefinementMethod%desc = &
301 : SKC_"The simulation specification `outputSampleRefinementMethod` is a scalar string of maximum length `"//getStr(len(outputSampleRefinementMethod, IK))//&
302 : SKC_"` representing the method of computing the (integrated) AutoCorrelation Time (ACT) to be used in the simulation for refining &
303 : &the final output MCMC chain and sample. If specified within an external input file, it must be either singly or doubly quoted. &
304 : &Methods that are currently supported include:"//NL2//&
305 : SKC_"+ `outputSampleRefinementMethod = '"//spec%outputSampleRefinementMethod%batchMeans//SKC_"'`"//NL2//&
306 : SKC_" This method of computing the integrated Autocorrelation Time is based on the approach described in"//NL2//&
307 : SKC_" SCHMEISER, B., 1982, Batch size effects in the analysis of simulation output, Oper. Res. 30 556-568."//NL2//&
308 : SKC_" The batch sizes in the BatchMeans method are chosen to be `int(N^(2/3))`, where `N` is the length of the MCMC chain. &
309 : &As long as the batch size is larger than the ACT of the chain and there are significantly more than `10` batches, &
310 : &the BatchMeans method will provide reliable estimates of the ACT. &
311 : &Note that the refinement strategy involves two separate phases of sample decorrelation. At the first stage, &
312 : &the Markov chain is decorrelated recursively (for as long as needed) based on the ACT of its compact format, &
313 : &where only the uniquely visited states are kept in the (compact) chain. Once the Markov chain is refined &
314 : &such that its compact format is fully decorrelated, the second phase of the decorrelation begins, during which &
315 : &the Markov chain is decorrelated based on the ACT of the chain in its verbose (Markov) format. This process &
316 : &is repeated recursively for as long as residual autocorrelation exists in the refined sample."//NL2//&
317 : SKC_"+ `outputSampleRefinementMethod = '"//spec%outputSampleRefinementMethod%batchMeans//SKC_"-compact'`"//NL2//&
318 : SKC_" This is the same as the first case in the above, except that only the first phase of the sample refinement &
319 : &described in the above will be performed; that is, the (verbose) Markov chain is refined only based on the &
320 : &ACT computed from the compact format of the Markov chain. This will lead to a larger, final, refined sample. &
321 : &However, the final sample will likely not be fully decorrelated."//NL2//&
322 : SKC_"+ `outputSampleRefinementMethod = '"//spec%outputSampleRefinementMethod%batchMeans//SKC_"-verbose'`"//NL2//&
323 : SKC_" This is the same as the first case in the above, except that only the second phase of the sample refinement &
324 : &described in the above will be performed; that is, the (verbose) Markov chain is refined only based on the ACT &
325 : &computed from the verbose format of the Markov chain. While the resulting refined sample will be fully decorrelated, &
326 : &the size of the refined sample may be smaller than the default choice in the first case in the above."//NL2//&
327 : SKC_"Note that to obtain i.i.d. samples from a multidimensional chain, the sampler will use the average of the &
328 : &ACT among all dimensions of the chain to refine the chain. If the maximum, minimum, or the median of IACs is preferred &
329 : &add `'-max'` (or `'-maximum'`), `'-min'` (or `'-minimum'`), `'-med'` (or `'-median'`), respectively, to the value of &
330 : &`outputSampleRefinementMethod`. For example, "//NL2//&
331 : SKC_"+ `outputSampleRefinementMethod = '"//spec%outputSampleRefinementMethod%batchMeans//SKC_"-max'`"//NL2//&
332 : SKC_"or, "//NL2//&
333 : SKC_"+ `outputSampleRefinementMethod = '"//spec%outputSampleRefinementMethod%batchMeans//SKC_"-compact-max'`"//NL2//&
334 : SKC_"or, "//NL2//&
335 : SKC_"+ `outputSampleRefinementMethod = '"//spec%outputSampleRefinementMethod%batchMeans//SKC_"-max-compact'`"//NL2//&
336 : SKC_"Note that the specified `outputSampleRefinementCount` is used only when the condition `outputSampleSize < 0` holds. &
337 : &Otherwise, it is ignored. The default value for `outputSampleRefinementMethod` is `'"//spec%outputSampleRefinementMethod%def//SKC_"'`. &
338 14 : &Note that the input values are case-INsensitive and white-space characters are ignored."
339 : !!$omp master
340 14 : outputSampleRefinementMethod = spec%outputSampleRefinementMethod%null
341 : !!$omp end master
342 : end block outputSampleRefinementMethod_block
343 :
344 : proposal_block: block
345 : use pm_sampling_scio, only: proposal
346 14 : spec%proposal%is%uniform = .false._LK
347 14 : spec%proposal%is%normal = .false._LK
348 14 : spec%proposal%uniform = "uniform"
349 14 : spec%proposal%normal = "normal"
350 14 : spec%proposal%def = spec%proposal%normal
351 14 : spec%proposal%null = repeat(SUB, len(proposal, IK))
352 : spec%proposal%desc = &
353 : SKC_"The simulation specification `proposal` is a scalar string of maximum length `"//getStr(len(proposal, IK))//SKC_"` containing &
354 : &the name of the proposal distribution for the MCMC sampler. When specified from within an external input file, it must &
355 : &be singly or doubly quoted. Options that are currently supported include:"//NL2//&
356 : SKC_"+ `proposal = '"//spec%proposal%normal//SKC_"'`"//NL2//&
357 : SKC_" This is equivalent to the multivariate normal distribution, &
358 : &which is the most widely-used proposal model along with MCMC samplers."//NL2//&
359 : SKC_"+ `proposal = '"//spec%proposal%uniform//SKC_"'`"//NL2//&
360 : SKC_" The proposals will be drawn uniformly from within a ndim-dimensional ellipsoid whose covariance matrix &
361 : &and scale are initialized by the user and optionally adaptively updated throughout the simulation."//NL2//&
362 14 : SKC_"The default value for `proposal` is `'"//spec%proposal%def//SKC_"'`."
363 : !!$omp master
364 14 : proposal = spec%proposal%null
365 : !!$omp end master
366 : end block proposal_block
367 :
368 : proposalCorMat_block: block
369 : use pm_except, only: setNAN
370 : use pm_sampling_scio, only: proposalCorMat
371 : use pm_matrixInit, only: getMatInit, uppLowDia
372 272 : spec%proposalCorMat%def = getMatInit([ndim, ndim], uppLowDia, 0._RKC, 0._RKC, 1._RKC)
373 : spec%proposalCorMat%desc = &
374 : SKC_"The simulation specification `proposalCorMat` is a positive-definite square matrix of type `real` of the highest precision &
375 : &available within the ParaMonte library matrix of size `(ndim, ndim)`, where `ndim` is the dimension of the sampling space. &
376 : &It serves as the best-guess starting correlation matrix of the proposal distribution used by the sampler. &
377 : &It is used (along with the input vector `proposalStdVec`) to construct the covariance matrix of the proposal &
378 : &distribution when the input covariance matrix is missing in the input list of variables. &
379 : &If the covariance matrix is specified as input to the sampler, any input values for `proposalCorMat`, &
380 : &and `proposalStdVec` will be automatically ignored. Specifying `proposalCorMat` along with `proposalStdVec` is especially &
381 : &useful when obtaining the best-guess covariance matrix is not trivial. &
382 14 : &The default value for `proposalCorMat` is a square Identity matrix of rank `ndim`."
383 : !!$omp master
384 14 : if (allocated(proposalCorMat)) deallocate(proposalCorMat)
385 27 : allocate(proposalCorMat(ndim, ndim))
386 122 : call setNAN(proposalCorMat)
387 : !!$omp end master
388 : end block proposalCorMat_block
389 :
390 : proposalCovMat_block: block
391 : use pm_except, only: setNAN
392 : use pm_sampling_scio, only: proposalCovMat
393 : use pm_matrixInit, only: getMatInit, uppLowDia
394 272 : spec%proposalCovMat%def = getMatInit([ndim, ndim], uppLowDia, 0._RKC, 0._RKC, 1._RKC)
395 : ! This must be set here. It is important for the proper setting from inputFile and inputArg.
396 14 : spec%proposalCovMat%isUserSet = .false._LK
397 : spec%proposalCovMat%desc = &
398 : SKC_"The simulation specification `proposalCovMat` is a square positive-definite matrix of type `real` of the highest precision &
399 : &available within the ParaMonte library, of shape `(ndim, ndim)`, where `ndim` is the number of dimensions of the sampling space. &
400 : &It serves as the best-guess starting covariance matrix of the proposal distribution. &
401 : &To bring the sampling efficiency of the sampler to within the desired requested range, the covariance matrix will &
402 : &be adaptively updated throughout the simulation, according to the user-specified schedule. If `proposalCovMat` &
403 : &is not provided by the user or it is completely missing from the input file, its value will be automatically &
404 : &computed via the input variables `proposalCorMat` and `proposalStdVec` (or via their default values, if not provided). &
405 : &If the simulation specification `outputStatus` is set to ""extend"" and a successful prior simulation run exists, &
406 : &then `proposalCovMat` will be set to the covariance matrix of the output sample from the most recent simulation run. &
407 : &In this case, the computed `proposalCovMat` will override any user-specified value. &
408 14 : &Otherwise, the default value for `proposalCovMat` is a square Identity matrix of rank `ndim`."
409 : !!$omp master
410 14 : if (allocated(proposalCovMat)) deallocate(proposalCovMat)
411 27 : allocate(proposalCovMat(ndim, ndim))
412 122 : call setNAN(proposalCovMat)
413 : !!$omp end master
414 : end block proposalCovMat_block
415 :
416 : proposalScaleFactor_block: block
417 : use pm_sampling_scio, only: proposalScaleFactor
418 14 : spec%proposalScaleFactor%strdef = SKC_"gelman"
419 14 : spec%proposalScaleFactor%valdef = 2.38_RKC / sqrt(real(ndim, RKC)) ! Gelman, Roberts, Gilks (1996): Efficient Metropolis Jumping Rules.
420 14 : spec%proposalScaleFactor%null = repeat(SUB, len(proposalScaleFactor, IK))
421 : spec%proposalScaleFactor%desc = &
422 : SKC_"The simulation specification `proposalScaleFactor` is a scalar string of maximum length `"//getStr(len(proposalScaleFactor, IK))//SKC_"` &
423 : &containing a positive real-valued number whose square will be multiplied with the covariance matrix of the proposal distribution &
424 : &of the MCMC sampler to shrink or enlarge it. In other words, the proposal distribution will be scaled in every direction &
425 : &by the specified numeric value of `proposalScaleFactor`. It can also be given in units of the string keyword `'gelman'` &
426 : &(which is case-INsensitive) after the paper:"//NL2//&
427 : SKC_" Gelman, Roberts, and Gilks (1996), Efficient Metropolis Jumping Rules."//NL2//&
428 : SKC_"The paper finds that the optimal scaling factor for a Multivariate Gaussian proposal distribution for the Metropolis-Hastings &
429 : &Markov Chain Monte Carlo sampling of a target Multivariate Normal Distribution of dimension `ndim` is given by:"//NL2//&
430 : SKC_" proposalScaleFactor = 2.38 / sqrt(ndim) , in the limit of ndim -> Infinity."//NL2//&
431 : SKC_"Multiples of the Gelman scale factors are also acceptable as input and can be specified like the following examples:"//NL2//&
432 : SKC_"+ `proposalScaleFactor = '1'`"//NL2//&
433 : SKC_" multiplies the ndim-dimensional proposal covariance matrix by 1, &
434 : &essentially no change occurs to the covariance matrix."//NL2//&
435 : SKC_"+ `proposalScaleFactor = ""1""`"//NL2//&
436 : SKC_" same as the previous example. The double-quotation marks act the same way as single-quotation marks."//NL2//&
437 : SKC_"+ `proposalScaleFactor = '2.5'`"//NL2//&
438 : SKC_" multiplies the ndim-dimensional proposal covariance matrix by 2.5."//NL2//&
439 : SKC_"+ `proposalScaleFactor = '2.5*Gelman'`"//NL2//&
440 : SKC_" multiplies the `ndim`-dimensional proposal covariance matrix by 2.5 * 2.38/sqrt(ndim)."//NL2//&
441 : SKC_"+ `proposalScaleFactor = ""2.5 * gelman""`"//NL2//&
442 : SKC_" same as the previous example but with double-quotation marks. space characters are ignored."//NL2//&
443 : SKC_"+ `proposalScaleFactor = ""2.5 * gelman*gelman*2""`"//NL2//&
444 : SKC_" equivalent to gelmanFactor-squared multiplied by `5`."//NL2//&
445 : SKC_"Note, however, that the result of Gelman et al. paper applies only to multivariate normal proposal distributions, in the limit &
446 : &of infinite dimensions. Therefore, care must be taken when using Gelman's scaling factor with non-Gaussian proposals and target &
447 : &objective functions. Only the product symbol `*` can be parsed in the string value of `proposalScaleFactor`. &
448 : &The presence of other mathematical symbols or multiple appearances of the product symbol will lead to a simulation crash. &
449 : &Also, note that the prescription of an acceptance range specified by the input variable `targetAcceptanceRate` will lead &
450 : &to dynamic modification of the initial input value of `proposalScaleFactor` throughout sampling for `proposalAdaptationCount` times. &
451 : &The default string value for `proposalScaleFactor` is `""gelman""` (for all proposal distributions), &
452 14 : &which is subsequently converted to `2.38 / sqrt(ndim)`."
453 : !!$omp master
454 14 : proposalScaleFactor = spec%proposalScaleFactor%null
455 : !!$omp end master
456 : end block proposalScaleFactor_block
457 :
458 : proposalStart_block: block
459 : use pm_sampling_scio, only: proposalStart
460 : spec%proposalStart%desc = &
461 : SKC_"The simulation specification `proposalStart` is a vector type of `real` of the highest precision available &
462 : &within the ParaMonte library of length `ndim` where `ndim` is the dimension of the domain of the objective function. &
463 : &For every element of `proposalStart` that is not provided as input, the default value will be the center of the sampling domain &
464 : &as determined by `domainCubeLimitLower` and `domainCubeLimitUpper` input specifications. &
465 : &If the condition `proposalStartRandomized` is set to the logical/Boolean true value, then the missing &
466 : &elements of `proposalStart` will be initialized to values drawn randomly from within the corresponding &
467 : &ranges specified by the input variables `proposalStartDomainCubeLimitLower` and `proposalStartDomainCubeLimitUpper`. &
468 : &If the simulation specification `outputStatus` is set to ""extend"" and a prior successful simulation run exists, &
469 : &then `proposalStart` will be set to the average of the sampled states of the most recent successful run. &
470 14 : &In this case, any user-specified value will be overridden by the computed `proposalStart`."
471 : !!$omp master
472 14 : if (allocated(proposalStart)) deallocate(proposalStart)
473 14 : allocate(proposalStart(ndim))
474 42 : call setNAN(proposalStart)
475 : !!$omp end master
476 : end block proposalStart_block
477 :
478 : proposalStartDomainCubeLimitLower_block: block
479 : use pm_sampling_scio, only: proposalStartDomainCubeLimitLower
480 14 : spec%proposalStartDomainCubeLimitLower%def = spec%domainCubeLimitLower%def
481 : spec%proposalStartDomainCubeLimitLower%desc = &
482 : SKC_"The simulation specification `proposalStartDomainCubeLimitLower` is a vector of type `real` of the highest precision available &
483 : &in the ParaMonet library of size `ndim` is the number of dimensions of the domain of the objective function. It contains the &
484 : &lower boundaries of the cubical domain from which the starting point(s) of the MCMC chain(s) will be initialized randomly &
485 : &(only if requested via the input variable `proposalStartRandomized`). &
486 : &This happens only when some or all of the input specification `proposalStart` elements are missing. &
487 : &In such cases, every missing value of the input `proposalStart` will be set to the center point between &
488 : &`proposalStartDomainCubeLimitLower` and `proposalStartDomainCubeLimitUpper` in the corresponding dimension. &
489 : &If `proposalStartRandomized` is set to the logical/Boolean true value, then the missing elements of &
490 : &`proposalStart` will be initialized to values drawn randomly from within the corresponding ranges &
491 : &whose lower limits are specified by the input `proposalStartDomainCubeLimitLower`. &
492 : &When specified from within an external input file to the sampler, it is also possible to assign only select &
493 : &values of `proposalStartDomainCubeLimitLower` and leave the rest of the components to be assigned the default value. &
494 : &For example, having the following inside the input file, "//NL2//&
495 : SKC_"+ `proposalStartDomainCubeLimitLower(3:5) = -100`"//NL2//&
496 : SKC_" will only set the lower limits of the third, fourth, and the fifth dimensions to -100, or,"//NL2//&
497 : SKC_"+ `proposalStartDomainCubeLimitLower(1) = -100, proposalStartDomainCubeLimitLower(2) = -1.e6`"//NL2//&
498 : SKC_" will set the lower limit on the first dimension to -100, and 1.e6 on the second dimension, or,"//NL2//&
499 : SKC_"+ `proposalStartDomainCubeLimitLower = 3*-2.5e100`"//NL2//&
500 : SKC_" will only set the lower limits on the first, second, and the third dimensions to `-2.5*10^100`, &
501 : &while the rest of the lower limits for the missing dimensions will be automatically set to the default value."//NL2//&
502 14 : SKC_"The default for all `proposalStartDomainCubeLimitLower` elements are taken from the corresponding elements of `domainCubeLimitLower`."
503 : !!$omp master
504 14 : call setResized(proposalStartDomainCubeLimitLower, ndim)
505 42 : call setNAN(proposalStartDomainCubeLimitLower)
506 : !!$omp end master
507 : end block proposalStartDomainCubeLimitLower_block
508 :
509 : proposalStartDomainCubeLimitUpper_block: block
510 : use pm_sampling_scio, only: proposalStartDomainCubeLimitUpper
511 14 : spec%proposalStartDomainCubeLimitUpper%def = spec%domainCubeLimitUpper%def
512 : spec%proposalStartDomainCubeLimitUpper%desc = &
513 : SKC_"The simulation specification `proposalStartDomainCubeLimitUpper` is a vector of type `real` of the highest precision available &
514 : &in the ParaMonet library of size `ndim` is the number of dimensions of the domain of the objective function. It contains the &
515 : &upper boundaries of the cubical domain from which the starting point(s) of the MCMC chain(s) will be initialized randomly &
516 : &(only if requested via the input variable `proposalStartRandomized`). &
517 : &This happens only when some or all of the input specification `proposalStart` elements are missing. &
518 : &In such cases, every missing value of the input `proposalStart` will be set to the center point between &
519 : &`proposalStartDomainCubeLimitLower` and `proposalStartDomainCubeLimitUpper` in the corresponding dimension. &
520 : &If `proposalStartRandomized` is set to the logical/Boolean true value, then the missing elements of &
521 : &`proposalStart` will be initialized to values drawn randomly from within the corresponding ranges &
522 : &whose upper limits are specified by the input `proposalStartDomainCubeLimitUpper`. &
523 : &When specified from within an external input file to the sampler, it is also possible to assign only select &
524 : &values of `proposalStartDomainCubeLimitUpper` and leave the rest of the components to be assigned the default value. &
525 : &For example, having the following inside the input file, "//NL2//&
526 : SKC_"+ `proposalStartDomainCubeLimitUpper(3:5) = -100`"//NL2//&
527 : SKC_" will only set the upper limits of the third, fourth, and the fifth dimensions to `-100`, or,"//NL2//&
528 : SKC_"+ `proposalStartDomainCubeLimitUpper(1) = -100, proposalStartDomainCubeLimitUpper(2) = -1.e6`"//NL2//&
529 : SKC_" will set the upper limit on the first dimension to -100, and 1.e6 on the second dimension, or,"//NL2//&
530 : SKC_"+ `proposalStartDomainCubeLimitUpper = 3*-2.5e100`"//NL2//&
531 : SKC_" will only set the upper limits on the first, second, and the third dimensions to `-2.5*10**100`, while &
532 : &the rest of the upper limits for the missing dimensions will be automatically set to the default value."//NL2//&
533 : SKC_"The default values for all elements of `proposalStartDomainCubeLimitUpper` are &
534 14 : &taken from the corresponding values in the input variable `domainCubeLimitUpper`."
535 : !!$omp master
536 14 : call setResized(proposalStartDomainCubeLimitUpper, ndim)
537 42 : call setNAN(proposalStartDomainCubeLimitUpper)
538 : !!$omp end master
539 : end block proposalStartDomainCubeLimitUpper_block
540 :
541 : proposalStartRandomized_block: block
542 : use pm_sampling_scio, only: proposalStartRandomized
543 14 : spec%proposalStartRandomized%def = .false._LK
544 : spec%proposalStartRandomized%desc = &
545 : SKC_"The simulation specification `proposalStartRandomized` is scalar of type `logical` (Boolean). &
546 : &If `true` (or `.true.` or `TRUE` or `.t.` from within an external input file), then the variable `proposalStart` &
547 : &will be initialized randomly for each MCMC chain that is to be generated by the sampler. The random values will be &
548 : &drawn from the specified or the default domain of `proposalStart`, given by `proposalStartDomainCubeLimitLower` and &
549 : &`proposalStartDomainCubeLimitUpper` variable. Note that the value of `proposalStart`, if provided, has precedence over &
550 : &random initialization. In other words, only uninitialized elements of `proposalStart` will be randomly initialized only &
551 : &if `proposalStartRandomized` is set to the logical true value. Note that even if `proposalStart` is randomly initialized, &
552 : &its random value will be deterministic between different independent simulation runs if the input variable `randomSeed` &
553 14 : &is specified by the user. The default value is `"//getStr(spec%proposalStartRandomized%def)//SKC_"`."
554 : !!$omp master
555 14 : proposalStartRandomized = spec%proposalStartRandomized%def
556 : !!$omp end master
557 : end block proposalStartRandomized_block
558 :
559 : proposalStdVec_block: block
560 : use pm_arrayFill, only: getFilled
561 : use pm_sampling_scio, only: proposalStdVec
562 56 : spec%proposalStdVec%def = getFilled(1._RKC, ndim)
563 : spec%proposalStdVec%desc = &
564 : SKC_"The simulation specification `proposalStdVec` is a positive-valued vector of type `real` of the highest precision available &
565 : &within the ParaMonte library, of size `ndim`, where `ndim` is the dimension of the domain of the objective function. &
566 : &It serves as the best-guess starting standard deviation for each component of the proposal distribution. &
567 : &If the initial covariance matrix (`proposalCovMat`) is missing as an input specification to the sampler, &
568 : &then `proposalStdVec` (along with the specified `proposalCorMat`) will be used to construct &
569 : &the initial covariance matrix of the proposal distribution of the MCMC sampler. &
570 : &However, if `proposalCovMat` is specified for the sampler, then the input `proposalStdVec` and `proposalCorMat` will &
571 : &be completely ignored, and the input value for `proposalCovMat` will be used to construct the initial covariance matrix &
572 14 : &of the proposal distribution. The default value of `proposalStdVec` is a vector of size `ndim` of unit values (i.e., ones)."
573 : !!$omp master
574 14 : call setResized(proposalStdVec, ndim)
575 42 : call setNAN(proposalStdVec)
576 : !!$omp end master
577 : end block proposalStdVec_block
578 :
579 : !!$omp barrier
580 :
581 70 : end function
582 :
583 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
584 :
585 13 : function set(spec, sampler) result(err)
586 : #if __INTEL_COMPILER && DLL_ENABLED && (_WIN32 || _WIN64)
587 : !DEC$ ATTRIBUTES DLLEXPORT :: set
588 : #endif
589 : use pm_err, only: err_type
590 : use pm_sampling, only: paramcmc_type
591 : use pm_sampling, only: sampler_type
592 : class(specmcmc_type), intent(inout) :: spec
593 : class(sampler_type), intent(in), optional :: sampler
594 : type(err_type) :: err
595 :
596 : character(*, SK), parameter :: PROCEDURE_NAME = MODULE_NAME//"@set()"
597 :
598 : select type(sampler)
599 : type is (paramcmc_type)
600 :
601 13 : err = spec%specbase_type%set(sampler%sampler_type)
602 13 : if (err%occurred) return
603 :
604 : outputChainSize_block: block
605 : use pm_sampling_scio, only: outputChainSize
606 13 : if (spec%overridable .and. allocated(sampler%outputChainSize)) then
607 9 : spec%outputChainSize%val = sampler%outputChainSize
608 : else
609 4 : spec%outputChainSize%val = outputChainSize
610 : end if
611 13 : if (spec%outputChainSize%val == spec%outputChainSize%null) spec%outputChainSize%val = spec%outputChainSize%def
612 : end block outputChainSize_block
613 :
614 : outputSampleRefinementCount_block: block
615 : use pm_sampling_scio, only: outputSampleRefinementCount
616 13 : if (spec%overridable .and. allocated(sampler%outputSampleRefinementCount)) then
617 0 : spec%outputSampleRefinementCount%val = sampler%outputSampleRefinementCount
618 : else
619 13 : spec%outputSampleRefinementCount%val = outputSampleRefinementCount
620 : end if
621 13 : if (spec%outputSampleRefinementCount%val == spec%outputSampleRefinementCount%null) spec%outputSampleRefinementCount%val = spec%outputSampleRefinementCount%def
622 13 : spec%outputSampleRefinementCount%str = getStr(spec%outputSampleRefinementCount%val)
623 : end block outputSampleRefinementCount_block
624 :
625 13 : outputSampleRefinementMethod_block: block
626 : use pm_sampling_scio, only: outputSampleRefinementMethod
627 : character(:,SKC), allocatable :: lowval
628 13 : if (spec%overridable .and. allocated(sampler%outputSampleRefinementMethod)) then
629 0 : spec%outputSampleRefinementMethod%val = trim(adjustl(getRemoved(sampler%outputSampleRefinementMethod, SKC_" ")))
630 : else
631 13 : spec%outputSampleRefinementMethod%val = trim(adjustl(getRemoved(outputSampleRefinementMethod, SKC_" ")))
632 : end if
633 13 : if (spec%outputSampleRefinementMethod%val == trim(adjustl(spec%outputSampleRefinementMethod%null))) spec%outputSampleRefinementMethod%val = spec%outputSampleRefinementMethod%def
634 13 : lowval = getStrLower(spec%outputSampleRefinementMethod%val)
635 13 : spec%outputSampleRefinementMethod%isMaxCumSumACF = logical(lowval == SKC_"maxcumsumacf", LK)
636 13 : spec%outputSampleRefinementMethod%isBatchMeans = logical(0 < index(lowval, SKC_"batchmeans"), LK)
637 13 : spec%outputSampleRefinementMethod%isBatchMeansCompact = spec%outputSampleRefinementMethod%isBatchMeans .and. logical(0 < index(lowval, SKC_"compact"), LK)
638 39 : spec%outputSampleRefinementMethod%isBatchMeansVerbose = spec%outputSampleRefinementMethod%isBatchMeans .and. logical(0 < index(lowval, SKC_"verbose"), LK)
639 : end block outputSampleRefinementMethod_block
640 :
641 : proposal_block: block
642 : use pm_sampling_scio, only: proposal
643 13 : if (spec%overridable .and. allocated(sampler%proposal)) then
644 0 : spec%proposal%val = getStrLower(trim(adjustl(sampler%proposal)))
645 : else
646 13 : spec%proposal%val = getStrLower(trim(adjustl(proposal)))
647 : end if
648 13 : if (spec%proposal%val == spec%proposal%null) spec%proposal%val = spec%proposal%def
649 13 : spec%proposal%is%uniform = spec%proposal%val == spec%proposal%uniform
650 13 : spec%proposal%is%normal = spec%proposal%val == spec%proposal%normal
651 : end block proposal_block
652 :
653 : proposalCorMat_block: block
654 : use pm_sampling_scio, only: proposalCorMat
655 13 : if (spec%overridable .and. allocated(sampler%proposalCorMat)) then
656 0 : spec%proposalCorMat%val = real(sampler%proposalCorMat, RKC)
657 : else
658 158 : spec%proposalCorMat%val = proposalCorMat
659 : end if
660 119 : where (isNAN(spec%proposalCorMat%val))
661 : spec%proposalCorMat%val = spec%proposalCorMat%def
662 : end where
663 : end block proposalCorMat_block
664 :
665 : proposalCovMat_block: block
666 : use pm_sampling_scio, only: proposalCovMat
667 : integer(IK) :: idim, jdim
668 13 : if (spec%overridable .and. allocated(sampler%proposalCovMat)) then
669 0 : spec%proposalCovMat%val = real(sampler%proposalCovMat, RKC)
670 : else
671 158 : spec%proposalCovMat%val = proposalCovMat
672 : end if
673 40 : do jdim = 1, spec%ndim%val
674 119 : do idim = 1, spec%ndim%val
675 106 : if (isNAN(spec%proposalCovMat%val(idim, jdim))) then
676 59 : spec%proposalCovMat%val(idim, jdim) = spec%proposalCovMat%def(idim, jdim)
677 : else
678 20 : spec%proposalCovMat%isUserSet = .true._LK
679 : end if
680 : end do
681 : end do
682 : end block proposalCovMat_block
683 :
684 : proposalScaleFactor_block: block
685 : use pm_sampling_scio, only: proposalScaleFactor
686 13 : if (spec%overridable .and. allocated(sampler%proposalScaleFactor)) then
687 1 : spec%proposalScaleFactor%str = trim(adjustl(sampler%proposalScaleFactor))
688 : else
689 12 : spec%proposalScaleFactor%str = trim(adjustl(proposalScaleFactor))
690 : end if
691 13 : if (spec%proposalScaleFactor%str == spec%proposalScaleFactor%null) then
692 10 : spec%proposalScaleFactor%val = spec%proposalScaleFactor%valdef
693 10 : spec%proposalScaleFactor%str = spec%proposalScaleFactor%strdef
694 : end if
695 : end block proposalScaleFactor_block
696 :
697 : proposalStart_block: block
698 : use pm_sampling_scio, only: proposalStart
699 13 : if (spec%overridable .and. allocated(sampler%proposalStart)) then
700 21 : spec%proposalStart%val = real(sampler%proposalStart, RKC)
701 : else
702 38 : spec%proposalStart%val = proposalStart
703 : end if
704 : end block proposalStart_block
705 :
706 : proposalStartDomainCubeLimitLower_block: block
707 : use pm_sampling_scio, only: proposalStartDomainCubeLimitLower
708 13 : if (spec%overridable .and. allocated(sampler%proposalStartDomainCubeLimitLower)) then
709 0 : spec%proposalStartDomainCubeLimitLower%val = real(sampler%proposalStartDomainCubeLimitLower, RKC)
710 : else
711 66 : spec%proposalStartDomainCubeLimitLower%val = proposalStartDomainCubeLimitLower
712 : end if
713 40 : where(isNAN(spec%proposalStartDomainCubeLimitLower%val))
714 : spec%proposalStartDomainCubeLimitLower%val = spec%domainCubeLimitLower%val
715 : end where
716 : end block proposalStartDomainCubeLimitLower_block
717 :
718 : proposalStartDomainCubeLimitUpper_block: block
719 : use pm_sampling_scio, only: proposalStartDomainCubeLimitUpper
720 13 : if (spec%overridable .and. allocated(sampler%proposalStartDomainCubeLimitUpper)) then
721 0 : spec%proposalStartDomainCubeLimitUpper%val = real(sampler%proposalStartDomainCubeLimitUpper, RKC)
722 : else
723 66 : spec%proposalStartDomainCubeLimitUpper%val = proposalStartDomainCubeLimitUpper
724 : end if
725 40 : where(isNAN(spec%proposalStartDomainCubeLimitUpper%val))
726 : spec%proposalStartDomainCubeLimitUpper%val = spec%domainCubeLimitUpper%val
727 : end where
728 : end block proposalStartDomainCubeLimitUpper_block
729 :
730 : proposalStartRandomized_block: block
731 : use pm_sampling_scio, only: proposalStartRandomized
732 13 : if (spec%overridable .and. allocated(sampler%proposalStartRandomized)) then
733 0 : spec%proposalStartRandomized%val = sampler%proposalStartRandomized
734 : else
735 13 : spec%proposalStartRandomized%val = proposalStartRandomized
736 : end if
737 : end block proposalStartRandomized_block
738 :
739 : proposalStdVec_block: block
740 : use pm_sampling_scio, only: proposalStdVec
741 13 : if (spec%overridable .and. allocated(sampler%proposalStdVec)) then
742 0 : spec%proposalStdVec%val = real(sampler%proposalStdVec, RKC)
743 : else
744 66 : spec%proposalStdVec%val = proposalStdVec
745 : end if
746 40 : where (isNAN(spec%proposalStdVec%val))
747 : spec%proposalStdVec%val = spec%proposalStdVec%def
748 : end where
749 : end block proposalStdVec_block
750 :
751 : ! Resolve the conflicting cases.
752 :
753 13 : if (spec%outputStatus%is%extend .and. 1 < spec%run%id) then
754 1 : block
755 : use pm_sampleMean, only: setMean
756 : use pm_sampleCov, only: setCov, uppDia
757 : use pm_io, only: getErrTableRead, LEN_IOMSG
758 : use pm_matrixCopy, only: setMatCopy, rdpack, upp, transHerm
759 : real(RKC), allocatable :: logFuncState(:,:)
760 1 : call setResized(err%msg, LEN_IOMSG)
761 1 : err%stat = getErrTableRead(spec%sampleFile%list(spec%run%id - 1)%val, logFuncState, sep = spec%outputSeparator%val, roff = 1_IK, iomsg = err%msg)
762 1 : err%occurred = err%stat /= 0
763 1 : if (err%occurred) then
764 : err%msg = PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SK_": Failed to read the sample file """//spec%sampleFile%list(spec%run%id - 1)%val//& ! LCOV_EXCL_LINE
765 : SK_""" contents from the previous simulation run. "//trim(err%msg) ! LCOV_EXCL_LINE
766 : return ! LCOV_EXCL_LINE
767 : end if
768 1 : call setMean(spec%proposalStart%val, logFuncState(:, 2 : spec%ndim%val + 1), dim = 1_IK)
769 1 : call setCov(spec%proposalCovMat%val, uppDia, logFuncState(:, 2 : spec%ndim%val + 1), dim = 1_IK)
770 1 : call setMatCopy(spec%proposalCovMat%val, rdpack, spec%proposalCovMat%val, rdpack, upp, transHerm)
771 : end block
772 : else
773 : block
774 : use pm_sampleCov, only: getCov, uppDia
775 : use pm_distUnif, only: setUnifRand
776 : integer(IK) :: idim
777 183 : if (.not. spec%proposalCovMat%isUserSet) spec%proposalCovMat%val = getCov(spec%proposalCorMat%val, uppDia, spec%proposalStdVec%val)
778 37 : do idim = 1, size(spec%proposalStart%val, 1, IK)
779 37 : if (isNAN(spec%proposalStart%val(idim))) then
780 18 : if (spec%proposalStartRandomized%val) then
781 0 : call setUnifRand(spec%rng, spec%proposalStart%val(idim), spec%proposalStartDomainCubeLimitLower%val(idim), spec%proposalStartDomainCubeLimitUpper%val(idim))
782 : else
783 18 : spec%proposalStart%val(idim) = 0.5_RKC * spec%proposalStartDomainCubeLimitLower%val(idim) + 0.5_RKC * spec%proposalStartDomainCubeLimitUpper%val(idim)
784 : end if
785 : end if
786 : end do
787 : end block
788 : end if
789 :
790 : ! open output files, report and sanitize.
791 :
792 13 : if (spec%image%is%leader) call spec%report() ! if (spec%run%is%new)
793 26 : call spec%sanitize(err)
794 :
795 : class default
796 : error stop "The input `sampler` must be of type `paramcmc_type`." ! LCOV_EXCL_LINE
797 : end select
798 :
799 26 : end function
800 :
801 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
802 :
803 13 : subroutine report(spec)
804 : #if __INTEL_COMPILER && DLL_ENABLED && (_WIN32 || _WIN64)
805 : !DEC$ ATTRIBUTES DLLEXPORT :: report
806 : #endif
807 : class(specmcmc_type), intent(inout) :: spec
808 :
809 13 : call spec%disp%text%wrap(NL1//spec%method%val//SKC_".simulation.specifications.mcmc"//NL1)
810 :
811 : associate(ndim => spec%ndim%val, format => spec%reportFile%format%generic)
812 :
813 13 : call spec%disp%show("outputChainSize")
814 13 : call spec%disp%show(spec%outputChainSize%val, format = format)
815 13 : call spec%disp%note%show(spec%outputChainSize%desc)
816 :
817 13 : call spec%disp%show("outputSampleRefinementCount")
818 13 : call spec%disp%show(spec%outputSampleRefinementCount%val, format = format)
819 13 : call spec%disp%note%show(spec%outputSampleRefinementCount%desc)
820 :
821 13 : call spec%disp%show("outputSampleRefinementMethod")
822 13 : call spec%disp%show(spec%outputSampleRefinementMethod%val, format = format)
823 13 : call spec%disp%note%show(spec%outputSampleRefinementMethod%desc)
824 :
825 13 : call spec%disp%show("proposal")
826 13 : call spec%disp%show(spec%proposal%val, format = format)
827 13 : call spec%disp%note%show(spec%proposal%desc)
828 :
829 13 : call spec%disp%show("proposalCorMat")
830 13 : call spec%disp%show(spec%proposalCorMat%val, format = format)
831 13 : call spec%disp%note%show(spec%proposalCorMat%desc)
832 :
833 13 : call spec%disp%show("proposalCovMat")
834 13 : call spec%disp%show(spec%proposalCovMat%val, format = format)
835 13 : call spec%disp%note%show(spec%proposalCovMat%desc)
836 :
837 13 : call spec%disp%show("proposalScaleFactor")
838 13 : call spec%disp%show(spec%proposalScaleFactor%str//SKC_" (gelman(ndim = "//getStr(spec%ndim%val)//SKC_") = "//getStr(spec%proposalScaleFactor%valdef)//SKC_")", format = format)
839 13 : call spec%disp%note%show(spec%proposalScaleFactor%desc)
840 :
841 13 : call spec%disp%show("proposalStart")
842 39 : call spec%disp%show(reshape(spec%proposalStart%val, [ndim, 1_IK]), format = format)
843 13 : call spec%disp%note%show(spec%proposalStart%desc)
844 :
845 13 : call spec%disp%show("proposalStartDomainCubeLimitLower")
846 39 : call spec%disp%show(reshape(spec%proposalStartDomainCubeLimitLower%val, [ndim, 1_IK]), format = format)
847 13 : call spec%disp%note%show(spec%proposalStartDomainCubeLimitLower%desc)
848 :
849 13 : call spec%disp%show("proposalStartDomainCubeLimitUpper")
850 39 : call spec%disp%show(reshape(spec%proposalStartDomainCubeLimitUpper%val, [ndim, 1_IK]), format = format)
851 13 : call spec%disp%note%show(spec%proposalStartDomainCubeLimitUpper%desc)
852 :
853 13 : call spec%disp%show("proposalStartRandomized")
854 13 : call spec%disp%show(spec%proposalStartRandomized%val, format = format)
855 13 : call spec%disp%note%show(spec%proposalStartRandomized%desc)
856 :
857 13 : call spec%disp%show("proposalStdVec")
858 39 : call spec%disp%show(reshape(spec%proposalStdVec%val, [ndim, 1_IK]), format = format)
859 26 : call spec%disp%note%show(spec%proposalStdVec%desc)
860 :
861 : end associate
862 :
863 13 : end subroutine
864 :
865 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
866 :
867 13 : subroutine sanitize(spec, err)
868 : #if __INTEL_COMPILER && DLL_ENABLED && (_WIN32 || _WIN64)
869 : !DEC$ ATTRIBUTES DLLEXPORT :: sanitize
870 : #endif
871 : use pm_err, only: err_type
872 : type(err_type), intent(inout) :: err
873 : class(specmcmc_type), intent(inout) :: spec
874 : character(*,SKC), parameter :: PROCEDURE_NAME = MODULE_NAME//SKC_"@sanitize()"
875 :
876 : outputChainSize_block: block
877 13 : if (spec%outputChainSize%val < spec%ndim%val + 1_IK) then
878 0 : err%occurred = .true._LK
879 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. The input requested value for `outputChainSize` ("//getStr(spec%outputChainSize%val)//&
880 : SKC_") can neither be negative nor smaller than `ndim + 1`, where `ndim` represents the dimension of the domain of the objective &
881 : &function, here `ndim = "//getStr(spec%ndim%val)//SKC_". If you do not know an appropriate value for `outputChainSize`, &
882 0 : &drop it from the input list. The MCMC sampler will automatically assign an appropriate value to it."
883 : end if
884 : end block outputChainSize_block
885 :
886 : outputSampleRefinementCount_block: block
887 13 : if (spec%outputSampleRefinementCount%val < 0_IK) then
888 0 : err%occurred = .true._LK
889 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
890 : &The input value for variable `outputSampleRefinementCount` must be a non-negative integer. &
891 : &If you are unsure about the appropriate value for this variable, simply drop it from the input. &
892 0 : &The sampler will automatically assign an appropriate value to it."
893 : end if
894 : end block outputSampleRefinementCount_block
895 :
896 13 : outputSampleRefinementMethod_block: block
897 : character(:,SKC), allocatable :: lowerCaseVal
898 13 : lowerCaseVal = getStrLower(spec%outputSampleRefinementMethod%val)
899 : if (& ! LCOV_EXCL_LINE
900 : index(lowerCaseVal, getStrLower(spec%outputSampleRefinementMethod%maxCumSumACF)) == 0_IK & ! LCOV_EXCL_LINE
901 : .and. & ! LCOV_EXCL_LINE
902 : index(lowerCaseVal, getStrLower(spec%outputSampleRefinementMethod%batchMeans)) == 0_IK & ! LCOV_EXCL_LINE
903 : .and. & ! LCOV_EXCL_LINE
904 : index(lowerCaseVal, getStrLower(spec%outputSampleRefinementMethod%cutoffACF)) == 0_IK & ! LCOV_EXCL_LINE
905 26 : ) then
906 0 : err%occurred = .true._LK
907 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
908 : &The input requested method for the computation of the integrated AutoCorrelation Time ("//spec%outputSampleRefinementMethod%val//&
909 : SKC_") assigned to the variable `outputSampleRefinementMethod` cannot be anything other than "//spec%outputSampleRefinementMethod%batchMeans//SKC_". "//&
910 : ! " or "//spec%maxCumSumACF//". &
911 : SKC_"If you are unsure of the appropriate value for `outputSampleRefinementMethod`, &
912 13 : &drop it from the input list. The sampler will automatically assign an appropriate value to it."
913 : end if
914 : end block outputSampleRefinementMethod_block
915 :
916 : proposal_block: block
917 13 : if (.not. (spec%proposal%is%normal .or. spec%proposal%is%uniform)) then
918 0 : err%occurred = .true._LK
919 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
920 : &The input requested proposal model ("//spec%proposal%val//") is not supported. &
921 : &The variable proposal cannot be set to anything other than '"//&
922 : spec%proposal%normal//SKC_"', or '"//&
923 0 : spec%proposal%uniform//SKC_"'."
924 : end if
925 : end block proposal_block
926 :
927 : proposalCorMat_block: block
928 13 : if (.not. isMatClass(spec%proposalCorMat%val, posdefmat)) then
929 0 : err%occurred = .true._LK
930 0 : err%msg = err%msg//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. The specified `proposalCorMat` is not positive-definite: "//getStr(spec%proposalCorMat%val)//SKC_""
931 : end if
932 : end block proposalCorMat_block
933 :
934 : proposalCovMat_block: block
935 13 : if (.not. isMatClass(spec%proposalCovMat%val, posdefmat)) then
936 0 : err%occurred = .true._LK
937 0 : err%msg = err%msg//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. The specified `proposalCovMat` is not positive-definite: "//getStr(spec%proposalCovMat%val)//SKC_""
938 : end if
939 : end block proposalCovMat_block
940 :
941 13 : proposalScaleFactor_block: block
942 : use pm_val2real, only: setReal
943 : use pm_arraySplit, only: setSplit
944 : integer(IK), allocatable :: sindex(:,:)
945 : character(:,SKC), allocatable :: str
946 : real(RKC) :: conversion
947 : integer(IK) :: ipart
948 : ! First convert the proposalScaleFactor string to real value:
949 13 : str = getRemoved(spec%proposalScaleFactor%str, SKC_" ") ! remove the white spaces.
950 13 : if (len_trim(adjustl(str)) == 0_IK) then
951 0 : err%occurred = .true._LK
952 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
953 : &The input string value ("//spec%proposalScaleFactor%str//SKC_") for the variable `proposalScaleFactor` is empty. &
954 : &Ensure the specified string follows the syntax rules of the sampler for this variable. Otherwise drop &
955 0 : &it from the input list. The sampler will automatically assign an appropriate value to it."
956 : end if
957 : ! Now split the string by "*" to real coefficient and character (gelman) parts for further evaluations.
958 13 : call setSplit(sindex, str, sep = SKC_"*")
959 13 : spec%proposalScaleFactor%val = 1._RKC
960 31 : do ipart = 1, size(sindex, 2, IK)
961 31 : if (getStrLower(str(sindex(1, ipart) : sindex(2, ipart))) == SKC_"gelman") then
962 13 : spec%proposalScaleFactor%val = spec%proposalScaleFactor%val * spec%proposalScaleFactor%valdef
963 : else
964 5 : call setReal(conversion, str(sindex(1, ipart) : sindex(2, ipart)), iostat = err%stat)
965 23 : if (err%stat /= 0_IK) then
966 0 : err%occurred = .true._LK
967 : err%msg = err%msg//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred while reading real number. &
968 : &The input string value for the variable `proposalScaleFactor` ("//spec%proposalScaleFactor%str//SKC_") does not appear to follow &
969 : &the standard syntax rules of the sampler for this variable. The slice `"//str(sindex(1, ipart) : sindex(1, ipart))//& ! LCOV_EXCL_LINE
970 : SKC_"` cannot be parsed into any meaningful token. Please correct the input value, or drop it from the input list &
971 0 : &in which case, the sampler will automatically assign an appropriate value to it."
972 : else
973 5 : spec%proposalScaleFactor%val = spec%proposalScaleFactor%val * conversion
974 : end if
975 : end if
976 : end do
977 : ! Now check if the real value is positive
978 26 : if (spec%proposalScaleFactor%val <= 0_IK) then
979 0 : err%occurred = .true._LK
980 : err%msg = err%msg//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
981 : &The input string value ("""//spec%proposalScaleFactor%str//""") translates to a non-positive real value: "//getStr(spec%proposalScaleFactor%val)//SKC_". &
982 : &Make sure the input string follows the syntax rules of the sampler for this variable. &
983 0 : &Otherwise drop it from the input list. The sampler will automatically assign an appropriate value to it."
984 : end if
985 : end block proposalScaleFactor_block
986 :
987 : proposalStart_block: block
988 : use pm_sampling_scio, only: proposalStart
989 : !!$omp barrier
990 : !!$omp master
991 13 : if (allocated(proposalStart)) deallocate(proposalStart)
992 : !!$omp end master
993 : end block proposalStart_block
994 :
995 : proposalStartDomainCubeLimitLower_block: block
996 : use pm_sampling_scio, only: proposalStartDomainCubeLimitLower
997 : !!$omp barrier
998 : !!$omp master
999 13 : if (allocated(proposalStartDomainCubeLimitLower)) deallocate(proposalStartDomainCubeLimitLower)
1000 : !!$omp end master
1001 : end block proposalStartDomainCubeLimitLower_block
1002 :
1003 : proposalStartDomainCubeLimitUpper_block: block
1004 : use pm_sampling_scio, only: proposalStartDomainCubeLimitUpper
1005 : integer(IK) :: idim
1006 40 : do idim = 1, size(spec%proposalStartDomainCubeLimitUpper%val, kind = IK)
1007 : ! Check if the domain is set when random start point is requested.
1008 27 : if (spec%proposalStartRandomized%val .and. spec%proposalStartDomainCubeLimitUpper%val(idim) == spec%domainCubeLimitUpper%def) then
1009 0 : err%occurred = .true._LK
1010 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
1011 : &You have requested a random start point by setting `proposalStartRandomized` to TRUE while the &
1012 : &element #"//getStr(idim)//SKC_" of `proposalStartDomainCubeLimitLower` has not been preset to a finite value. &
1013 0 : &This information is essential otherwise, how could the sampler draw points randomly from within an unspecified domain?"
1014 : end if
1015 : ! The upper boundary of the domain of random-start-point must be smaller than the upper boundary of the target's domain.
1016 27 : if (spec%domainCubeLimitUpper%val(idim) < spec%proposalStartDomainCubeLimitUpper%val(idim)) then
1017 0 : err%occurred = .true._LK
1018 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. The component "//getStr(idim)//& ! LCOV_EXCL_LINE
1019 : SKC_" of the variable `proposalStartDomainCubeLimitUpper` ("//getStr(spec%proposalStartDomainCubeLimitUpper%val(idim))//SKC_") cannot be "//& ! LCOV_EXCL_LINE
1020 : SKC_"larger than the corresponding component of the variable domainCubeLimitUpper ("//& ! LCOV_EXCL_LINE
1021 : getStr(spec%domainCubeLimitUpper%val(idim))//SKC_"). If you do not know an appropriate &
1022 : &value to set for `proposalStartDomainCubeLimitUpper`, drop it from the input list. &
1023 0 : &The sampler will automatically assign an appropriate value to it."
1024 : end if
1025 : ! The upper boundary of the domain of random-start-point must be smaller than the corresponding lower boundary.
1026 40 : if (spec%proposalStartDomainCubeLimitUpper%val(idim) <= spec%proposalStartDomainCubeLimitLower%val(idim)) then
1027 0 : err%occurred = .true._LK
1028 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. The input upper limit value in the component "//getStr(idim)//&
1029 : SKC_" of the variable `proposalStartDomainCubeLimitUpper` cannot be smaller than or equal to the corresponding input &
1030 : &lower limit value in `proposalStartDomainCubeLimitLower`:"//NL1//&
1031 : SKC_" proposalStartDomainCubeLimitLower("//getStr(idim)//SKC_") = "//getStr(spec%proposalStartDomainCubeLimitLower%val(idim))//NL1//&
1032 0 : SKC_" proposalStartDomainCubeLimitUpper("//getStr(idim)//SKC_") = "//getStr(spec%proposalStartDomainCubeLimitUpper%val(idim))//SKC_""
1033 : end if
1034 :
1035 : end do
1036 : !!$omp barrier
1037 : !!$omp master
1038 13 : if (allocated(proposalStartDomainCubeLimitUpper)) deallocate(proposalStartDomainCubeLimitUpper)
1039 : !!$omp end master
1040 : end block proposalStartDomainCubeLimitUpper_block
1041 :
1042 : proposalStartRandomized_block: block
1043 : end block proposalStartRandomized_block
1044 :
1045 : proposalStdVec_block: block
1046 : integer(IK) :: idim
1047 40 : do idim = 1, spec%ndim%val
1048 40 : if (spec%proposalStdVec%val(idim) <= 0._RKC) then
1049 0 : err%occurred = .true._LK
1050 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
1051 : &The input requested value ("//getStr(spec%proposalStdVec%val(idim))//SKC_") for the component "//getStr(idim)//&
1052 0 : SKC_" of the variable proposalStdVec for the proposal distribution of the sampler must be a positive real number."
1053 : end if
1054 : end do
1055 : end block proposalStdVec_block
1056 :
1057 : ! Resolve conflicts.
1058 :
1059 : block
1060 : ! proposalStart must be within the domain.
1061 : integer(IK) :: idim
1062 40 : do idim = 1, size(spec%proposalStart%val, 1, IK)
1063 40 : if (spec%proposalStart%val(idim) < spec%domainCubeLimitLower%val(idim) .or. spec%domainCubeLimitUpper%val(idim) < spec%proposalStart%val(idim)) then
1064 0 : err%occurred = .true._LK
1065 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
1066 : &The input requested value for the component "//getStr(idim)//SKC_" of the vector `proposalStart` ("//& ! LCOV_EXCL_LINE
1067 : getStr(spec%proposalStart%val(idim))//SKC_") must be within the range of the sampling domain defined in the program: ["//& ! LCOV_EXCL_LINE
1068 : getStr([spec%domainCubeLimitLower%val(idim), spec%domainCubeLimitUpper%val(idim)])//SKC_"]. &
1069 : &If you do not know an appropriate value for proposalStart, drop it from the input list. &
1070 27 : &The sampler will automatically assign an appropriate value to it."
1071 : end if
1072 : end do
1073 : end block
1074 :
1075 : block
1076 : ! `proposalStartDomainCubeLimitLower` must be larger than `domainCubeLimitLower` when start point is randomized.
1077 : integer(IK) :: idim
1078 40 : do idim = 1, size(spec%proposalStartDomainCubeLimitLower%val, kind = IK)
1079 : ! check if the domain is set when random start point is requested.
1080 27 : if (spec%proposalStartRandomized%val .and. spec%proposalStartDomainCubeLimitLower%val(idim) == spec%domainCubeLimitLower%def) then
1081 0 : err%occurred = .true._LK
1082 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. &
1083 : &You have requested a random start point by setting `proposalStartRandomized` to the logical/Boolean true value &
1084 : &while the element #"//getStr(idim)//SKC_" of `proposalStartDomainCubeLimitLower` has not been preset to a finite value. &
1085 0 : &This information is essential otherwise, how could the sampler draw points randomly from within an unspecified domain?"
1086 : end if
1087 : ! check if the random start point domain is within the boundaries of the domain of the target.
1088 40 : if (spec%proposalStartDomainCubeLimitLower%val(idim) < spec%domainCubeLimitLower%val(idim)) then
1089 0 : err%occurred = .true._LK
1090 : err%msg = err%msg//NL2//PROCEDURE_NAME//getFine(__FILE__, __LINE__)//SKC_": Error occurred. The component "//getStr(idim)//& ! LCOV_EXCL_LINE
1091 : SKC_" of the variable `proposalStartDomainCubeLimitLower` ("//getStr(spec%proposalStartDomainCubeLimitLower%val(idim))//& ! LCOV_EXCL_LINE
1092 : SKC_") cannot be smaller than the corresponding component of the variable `domainCubeLimitLower` ("//& ! LCOV_EXCL_LINE
1093 : getStr(spec%domainCubeLimitLower%val(idim))//"). If you do not know &an appropriate value to &
1094 : &set for `proposalStartDomainCubeLimitLower`, drop it from the input list. The sampler will &
1095 0 : &automatically assign an appropriate value to it."
1096 : end if
1097 :
1098 : end do
1099 : end block
1100 :
1101 13 : end subroutine
1102 :
1103 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1104 :
1105 : !> Return the refined Markov chain, given the input Markov chain and its specifications.
1106 : !> This procedure is a method of the [sfc_type](@ref sfc_type) class.
1107 : !>
1108 : !> \param[inout] sfc : An object of class [sfc_type](@ref sfc_type).
1109 : !> \param[in] outputSampleSize : The requested refined sample size (**optional**). If the size of the refined sample is given as input,
1110 : !! then the requested sample is directly generated based on the input size.
1111 : !> \param[in] outputSampleRefinementCount : The maximum number of times the sample can be refined (**optional**, default = `Infinity`).
1112 : !! : For example, if set to 1, then only one round of refinement will be performed on the Markov chain.
1113 : !> \param[in] outputSampleRefinementMethod : The requested method of refining the sample (**optional**, default = "BatchMeans").
1114 13 : function getErrRefinement(sfc, sampleWeight, sampleLogFunc, sampleState, outputSampleRefinementCount, outputSampleRefinementMethod) result(err)
1115 : #if __INTEL_COMPILER && DLL_ENABLED && (_WIN32 || _WIN64)
1116 : !DEC$ ATTRIBUTES DLLEXPORT :: getErrRefinement
1117 : #endif
1118 : use, intrinsic :: iso_fortran_env, only: output_unit
1119 : use pm_arrayRefine, only: setRefined
1120 : use pm_arrayRebind, only: setRebound
1121 : use pm_arrayResize, only: setResized
1122 : use pm_sampleACT, only: getACT, batchMeans, batchMeansMax, cumSum, cumSumMax
1123 : use pm_sampleQuan, only: getQuan, piwilin
1124 :
1125 : type(err_type) :: err
1126 : class(sfcmcmc_type) , intent(inout) :: sfc
1127 : integer(IK) , intent(in) , contiguous :: sampleWeight(:)
1128 : real(RKC) , intent(in) , contiguous :: sampleLogFunc(:), sampleState(:,:)
1129 : integer(IK) , intent(in) , optional :: outputSampleRefinementCount
1130 : character(*, SK) , intent(in) , optional :: outputSampleRefinementMethod
1131 : real(RKC) , allocatable :: logFuncStateAxis(:) ! ensures contiguity.
1132 : character(*, SK) , parameter :: PROCEDURE_NAME = SK_"@getErrRefinement()"
1133 : integer(IK) :: idim, ndim, nsam, ndimp1, nrefAlloc, outputSampleRefinementCount_def
1134 : real(RKC) :: act, ndimp1inv
1135 : type(method_type) :: method
1136 :
1137 13 : CHECK_ASSERTION(__LINE__, 0_IK < size(sampleState, 1, IK), SK_"The condition `0 < size(sampleState, 1)` must hold. shape(sampleState) = "//getStr(shape(sampleState, IK)))
1138 78 : CHECK_ASSERTION(__LINE__, size(sampleWeight, 1, IK) == size(sampleState, 2, IK), SK_"The condition `size(sampleWeight) == size(sampleState, 2)` must hold. size(sampleWeight), shape(sampleState) = "//getStr([size(sampleWeight, 1, IK), shape(sampleState, IK)]))
1139 39 : CHECK_ASSERTION(__LINE__, size(sampleWeight, 1, IK) == size(sampleLogFunc, 1, IK), SK_"The condition `size(sampleWeight) == size(sampleLogFunc)` must hold. size(sampleWeight), size(sampleLogFunc) = "//getStr([size(sampleWeight, 1, IK), size(sampleLogFunc, 1, IK)]))
1140 13 : err = getErrRefinementMethod(method, outputSampleRefinementMethod)
1141 : ndim = size(sampleState, 1, IK)
1142 : nsam = size(sampleState, 2, IK)
1143 13 : ndimp1 = ndim + 1_IK
1144 :
1145 : ! This is to avoid memory overflow due to extremely large `outputSampleRefinementCount` requested by the user.
1146 :
1147 13 : ndimp1inv = 1._RKC / ndimp1
1148 : outputSampleRefinementCount_def = 20_IK ! this is a temporary maximum value, to be increased later if needed
1149 13 : if (present(outputSampleRefinementCount)) outputSampleRefinementCount_def = outputSampleRefinementCount
1150 13 : nrefAlloc = min(2_IK, outputSampleRefinementCount_def)
1151 :
1152 : ! Allocate components.
1153 :
1154 39 : call setRebound(sfc%act, [0_IK, 0_IK], [ndim, nrefAlloc])
1155 13 : call setRebound(sfc%nsam, 0_IK, nrefAlloc)
1156 13 : call setRebound(sfc%sumw, 0_IK, nrefAlloc)
1157 :
1158 13 : sfc%nref = 0_IK
1159 13 : sfc%nsam(sfc%nref) = size(sampleWeight, 1, IK)
1160 13 : call setResized(logFuncStateAxis, sfc%nsam(sfc%nref))
1161 39 : call setRebound(sfc%sampleLogFuncState, [0_IK, 1_IK], [ndim, nsam])
1162 1300003 : sfc%sampleLogFuncState(1 : ndim, :) = sampleState
1163 285031 : sfc%sampleLogFuncState(0, :) = sampleLogFunc
1164 285044 : sfc%sampleWeight = sampleWeight
1165 :
1166 : ! Perform chain refinement.
1167 :
1168 13 : blockSufficientChainSize: if (ndim < nsam) then
1169 :
1170 : loopRefinement: do
1171 :
1172 593466 : sfc%sumw(sfc%nref) = sum(sfc%sampleWeight(1 : sfc%nsam(sfc%nref)))
1173 :
1174 : ! Obtain the ACT for each individual axes.
1175 : block
1176 216 : do idim = 0, ndim
1177 2773783 : logFuncStateAxis(1 : sfc%nsam(sfc%nref)) = sfc%sampleLogFuncState(idim, 1 : sfc%nsam(sfc%nref))
1178 216 : if (method%isViaCompactChain) then
1179 112 : if (method%isBatchMeansMax) then
1180 0 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), batchMeansMax)
1181 112 : elseif (method%isBatchMeans) then
1182 112 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), batchMeans)
1183 0 : elseif (method%isCumSumMax) then
1184 0 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), cumSumMax)
1185 0 : elseif (method%isCumSum) then
1186 0 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), cumSum)
1187 : else
1188 : error stop PROCEDURE_NAME//SK_"Internal Error: Unrecognized ACT computation method." ! LCOV_EXCL_LINE
1189 : end if
1190 52 : elseif (method%isViaVerboseChain) then
1191 52 : if (method%isBatchMeansMax) then
1192 0 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), sfc%sampleWeight(1 : sfc%nsam(sfc%nref)), batchMeansMax)
1193 52 : elseif (method%isBatchMeans) then
1194 52 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), sfc%sampleWeight(1 : sfc%nsam(sfc%nref)), batchMeans)
1195 0 : elseif (method%isCumSumMax) then
1196 0 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), sfc%sampleWeight(1 : sfc%nsam(sfc%nref)), cumSumMax)
1197 0 : elseif (method%isCumSum) then
1198 0 : sfc%act(idim, sfc%nref) = getACT(logFuncStateAxis(1 : sfc%nsam(sfc%nref)), sfc%sampleWeight(1 : sfc%nsam(sfc%nref)), cumSum)
1199 : else
1200 : error stop PROCEDURE_NAME//SK_"Internal Error: Unrecognized ACT computation method." ! LCOV_EXCL_LINE
1201 : end if
1202 : else
1203 : error stop PROCEDURE_NAME//SK_": Internal error: Unrecognized chain format." ! LCOV_EXCL_LINE
1204 : endif
1205 : end do
1206 : end block
1207 :
1208 : ! act is guaranteed to be larger than `1` by definition.
1209 :
1210 52 : if (method%isAvg) then
1211 216 : act = sum(sfc%act(0 : ndim, sfc%nref)) * ndimp1inv
1212 0 : elseif (method%isMed) then
1213 0 : act = getQuan(piwilin, .5_RKC, sfc%act(0 : ndim, sfc%nref))
1214 0 : elseif (method%isMax) then
1215 0 : act = maxval(sfc%act(0 : ndim, sfc%nref))
1216 0 : elseif (method%isMin) then
1217 0 : act = minval(sfc%act(0 : ndim, sfc%nref))
1218 : else
1219 : error stop PROCEDURE_NAME//SK_": Internal error: Unrecognized refinement method." ! LCOV_EXCL_LINE
1220 : end if
1221 :
1222 : ! So far, we have computed the max ACT of the sample, but no refinement. Refine the sample only if needed.
1223 : ! The output sample size can be less than the requested outputSampleSize if nsam < outputSampleSize.
1224 : ! This is the right place to quit if needed, as all components for the current stage of the refinement are set.
1225 :
1226 52 : if (sfc%nref == outputSampleRefinementCount_def) exit loopRefinement
1227 52 : if (act < 2._RKC) then
1228 26 : if (method%isViaCompactChain .and. method%isViaVerboseChain) then
1229 13 : method%isViaCompactChain = .false._LK
1230 13 : cycle loopRefinement
1231 : else
1232 : exit loopRefinement
1233 : end if
1234 : end if
1235 :
1236 : ! Generate the refined sample, then put it back into sfc to start over again.
1237 :
1238 26 : sfc%nref = sfc%nref + 1_IK
1239 :
1240 : ! Reallocate to bigger array if needed.
1241 :
1242 26 : if (nrefAlloc < sfc%nref) then
1243 3 : nrefAlloc = min(nrefAlloc * 2, outputSampleRefinementCount_def)
1244 9 : call setRebound(sfc%act, [0_IK, 0_IK], [ndim, nrefAlloc])
1245 3 : call setRebound(sfc%nsam, 0_IK, nrefAlloc)
1246 3 : call setRebound(sfc%sumw, 0_IK, nrefAlloc)
1247 : end if
1248 : !block
1249 : ! use pm_io
1250 : ! call disp%show("act")
1251 : ! call disp%show( act )
1252 : ! print *, sfc%nref
1253 : ! print *, sfc%nsam(0 : sfc%nref - 1)
1254 : ! !call disp%show("sfc%nref")
1255 : ! !call disp%show( sfc%nref )
1256 : ! call disp%show("sfc%nsam(0 : sfc%nref - 1)")
1257 : ! call disp%show( sfc%nsam(0 : sfc%nref - 1) )
1258 : !end block
1259 0 : if (act < 2._RKC) cycle loopRefinement ! no need for refinement. should happen only when transitioning from compact to verbose.
1260 39 : call setRefined(sfc%sampleLogFuncState(:, 1 : sfc%nsam(sfc%nref - 1)), 2_IK, sfc%sampleWeight(1 : sfc%nsam(sfc%nref - 1)), skip = int(act, IK), rsize = sfc%nsam(sfc%nref))
1261 :
1262 : end do loopRefinement
1263 :
1264 0 : elseif (0_IK < nsam) then blockSufficientChainSize
1265 :
1266 : ! Return only the last sampled state if the input sample is too small.
1267 :
1268 0 : sfc%nsam(sfc%nref) = size(sampleWeight, 1, IK)
1269 0 : call setResized(logFuncStateAxis, sfc%nsam(sfc%nref))
1270 0 : call setRebound(sfc%sampleLogFuncState, [0_IK, 1_IK], [ndim, nsam])
1271 0 : sfc%sampleLogFuncState(1 : ndim, :) = sampleState
1272 0 : sfc%sampleLogFuncState(0, :) = sampleLogFunc
1273 0 : sfc%sampleWeight = sampleWeight
1274 :
1275 : end if blockSufficientChainSize
1276 :
1277 39 : call setRebound(sfc%sampleLogFuncState, [0_IK, 1_IK], [ndim, sfc%nsam(sfc%nref)])
1278 13 : call setResized(sfc%sampleWeight, sfc%nsam(sfc%nref))
1279 :
1280 13 : end function getErrRefinement
1281 :
1282 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1283 :
1284 13 : function getErrRefinementMethod(method, outputSampleRefinementMethod) result(err)
1285 :
1286 : use pm_val2int, only: getInt
1287 : use pm_strASCII, only: getStrLower
1288 : use pm_arrayReplace, only: getReplaced
1289 : type(method_type), intent(inout) :: method
1290 : character(*, SK), parameter :: PROCEDURE_NAME = SK_"@getErrRefinementMethod()"
1291 : character(*, SK), intent(in), optional :: outputSampleRefinementMethod
1292 : character(:, SK), allocatable :: outputSampleRefinementMethod_def
1293 : type(err_type) :: err
1294 : integer(IK) :: idim
1295 :
1296 13 : err%msg = ""
1297 : err%occurred = .false._LK
1298 :
1299 : ! Define the ACT computation method.
1300 :
1301 13 : if (present(outputSampleRefinementMethod)) then
1302 13 : outputSampleRefinementMethod_def = getStrLower(outputSampleRefinementMethod)
1303 13 : if (0 < index(outputSampleRefinementMethod_def, SK_"batchmeansmax")) then
1304 0 : method%isBatchMeansMax = .true._LK
1305 13 : elseif (0 < index(outputSampleRefinementMethod_def, SK_"batchmeans")) then
1306 13 : method%isBatchMeans = .true._LK
1307 0 : elseif (0 < index(outputSampleRefinementMethod_def, SK_"cumsummax")) then
1308 0 : method%isCumSumMax = .true._LK
1309 0 : elseif (0 < index(outputSampleRefinementMethod_def, SK_"cumsum")) then
1310 0 : method%isCumSum = .true._LK
1311 : else
1312 : err%occurred = .true._LK ! LCOV_EXCL_LINE
1313 : err%msg = getFine(__FILE__, __LINE__)//PROCEDURE_NAME//SK_": Unknown unsupported ACT computation method name: "//outputSampleRefinementMethod_def ! LCOV_EXCL_LINE
1314 : return ! LCOV_EXCL_LINE
1315 : end if
1316 13 : method%isViaCompactChain = 0 < index(outputSampleRefinementMethod_def, SK_"compact")
1317 13 : method%isViaVerboseChain = 0 < index(outputSampleRefinementMethod_def, SK_"verbose")
1318 : else
1319 0 : outputSampleRefinementMethod_def = SK_"batchmeans"
1320 0 : method%isBatchMeans = .true._LK
1321 : end if
1322 :
1323 : ! Define the chain types to use for the ACT computation.
1324 :
1325 13 : if (.not. (method%isViaCompactChain .or. method%isViaVerboseChain)) then
1326 13 : method%isViaCompactChain = .true._LK
1327 13 : method%isViaVerboseChain = .true._LK
1328 : end if
1329 :
1330 : ! Define the statistic to use for the collective ACT computation.
1331 :
1332 13 : outputSampleRefinementMethod_def = getReplaced(outputSampleRefinementMethod_def, SK_" ", SK_"-")
1333 13 : method%isMed = index(outputSampleRefinementMethod_def, SK_"-med") > 0 .or. index(outputSampleRefinementMethod_def, SK_"-median") > 0
1334 13 : method%isAvg = index(outputSampleRefinementMethod_def, SK_"-avg") > 0 .or. index(outputSampleRefinementMethod_def, SK_"-average") > 0
1335 13 : method%isMin = index(outputSampleRefinementMethod_def, SK_"-min") > 0 .or. index(outputSampleRefinementMethod_def, SK_"-minimum") > 0
1336 13 : method%isMax = index(outputSampleRefinementMethod_def, SK_"-max") > 0 .or. index(outputSampleRefinementMethod_def, SK_"-maximum") > 0
1337 65 : idim = sum(getInt([method%isAvg, method%isMed, method%isMax, method%isMin]))
1338 13 : err%occurred = 1_IK < idim
1339 13 : if (err%occurred) then
1340 : err%msg = getFine(__FILE__, __LINE__)//PROCEDURE_NAME//SK_": avg, med, min, max cannot be simultaneously specified in outputSampleRefinementMethod: "//outputSampleRefinementMethod_def ! LCOV_EXCL_LINE
1341 : return ! LCOV_EXCL_LINE
1342 : end if
1343 13 : if (idim == 0_IK) method%isAvg = .true._LK ! default method of ACT summarization.
1344 :
1345 26 : end function getErrRefinementMethod
1346 :
1347 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
