Line data Source code
1 : !%%%%%%%%%%%%%%%
2 : #if namelist_ENABLED
3 : !%%%%%%%%%%%%%%%
4 :
5 : ! specbase
6 : namelist /NAMELIST/ description
7 : namelist /NAMELIST/ domain
8 : namelist /NAMELIST/ domainAxisName
9 : namelist /NAMELIST/ domainBallCenter
10 : namelist /NAMELIST/ domainBallCorMat
11 : namelist /NAMELIST/ domainBallCovMat
12 : namelist /NAMELIST/ domainBallStdVec
13 : namelist /NAMELIST/ domainCubeLimitLower
14 : namelist /NAMELIST/ domainCubeLimitUpper
15 : namelist /NAMELIST/ domainErrCount
16 : namelist /NAMELIST/ domainErrCountMax
17 : namelist /NAMELIST/ inputFileHasPriority
18 : namelist /NAMELIST/ outputChainFileFormat
19 : namelist /NAMELIST/ outputColumnWidth
20 : namelist /NAMELIST/ outputFileName
21 : namelist /NAMELIST/ outputStatus
22 : namelist /NAMELIST/ outputPrecision
23 : namelist /NAMELIST/ outputReportPeriod
24 : namelist /NAMELIST/ outputRestartFileFormat
25 : namelist /NAMELIST/ outputSampleSize
26 : namelist /NAMELIST/ outputSeparator
27 : namelist /NAMELIST/ outputSplashMode
28 : namelist /NAMELIST/ parallelism
29 : namelist /NAMELIST/ parallelismMpiFinalizeEnabled
30 : namelist /NAMELIST/ parallelismNumThread
31 : namelist /NAMELIST/ randomSeed
32 : namelist /NAMELIST/ sysInfoFilePath
33 : namelist /NAMELIST/ targetAcceptanceRate
34 : ! specmcmc
35 : namelist /NAMELIST/ outputChainSize
36 : namelist /NAMELIST/ outputSampleRefinementCount
37 : namelist /NAMELIST/ outputSampleRefinementMethod
38 : namelist /NAMELIST/ proposal
39 : namelist /NAMELIST/ proposalCorMat
40 : namelist /NAMELIST/ proposalCovMat
41 : namelist /NAMELIST/ proposalScaleFactor
42 : namelist /NAMELIST/ proposalStart
43 : namelist /NAMELIST/ proposalStartDomainCubeLimitLower
44 : namelist /NAMELIST/ proposalStartDomainCubeLimitUpper
45 : namelist /NAMELIST/ proposalStartRandomized
46 : namelist /NAMELIST/ proposalStdVec
47 : ! specdram
48 : namelist /NAMELIST/ burninAdaptationMeasure
49 : namelist /NAMELIST/ proposalAdaptationCount
50 : namelist /NAMELIST/ proposalAdaptationCountGreedy
51 : namelist /NAMELIST/ proposalAdaptationPeriod
52 : namelist /NAMELIST/ proposalDelayedRejectionCount
53 : namelist /NAMELIST/ proposalDelayedRejectionScaleFactor
54 : ! specnest
55 : namelist /NAMELIST/ domainPartitionAdaptationCount
56 : namelist /NAMELIST/ domainPartitionAdaptationPeriod
57 : namelist /NAMELIST/ domainPartitionBiasCorrectionEnabled
58 : namelist /NAMELIST/ domainPartitionCountMax
59 : namelist /NAMELIST/ domainPartitionFactorExpansion
60 : namelist /NAMELIST/ domainPartitionFactorShrinkage
61 : namelist /NAMELIST/ domainPartitionKmeansClusterCountMax
62 : namelist /NAMELIST/ domainPartitionKmeansClusterSizeMin
63 : namelist /NAMELIST/ domainPartitionKmeansNormalizationEnabled
64 : namelist /NAMELIST/ domainPartitionKmeansNumFailMax
65 : namelist /NAMELIST/ domainPartitionKmeansNumRecursionMax
66 : namelist /NAMELIST/ domainPartitionKmeansNumTry
67 : namelist /NAMELIST/ domainPartitionKvolumeNumRecursionMax
68 : namelist /NAMELIST/ domainPartitionKvolumeWeightExponent
69 : namelist /NAMELIST/ domainPartitionMethod
70 : namelist /NAMELIST/ domainPartitionObject
71 : namelist /NAMELIST/ domainPartitionOptimizationScaleEnabled
72 : namelist /NAMELIST/ domainPartitionOptimizationShapeEnabled
73 : namelist /NAMELIST/ domainPartitionOptimizationShapeScaleEnabled
74 : namelist /NAMELIST/ domainPartitionScaleFactor
75 : namelist /NAMELIST/ domainSampler
76 : namelist /NAMELIST/ liveSampleSize
77 : namelist /NAMELIST/ tolerance
78 :
79 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 : #elif ParaDISE_ENABLED || ParaDRAM_ENABLED || ParaNest_ENABLED
81 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82 :
83 : #define RETURN_IF_FAILED(LINE,FAILED,MSG) \
84 : if (FAILED) then; \
85 : err%occurred = .true._LK; \
86 : err%msg = PROCEDURE_NAME//getFine(__FILE__, LINE)//SK_": "//trim(MSG); \
87 : if (.not. opened) close(unit, iostat = err%stat); \
88 : return; \
89 : end if;
90 : #if ParaDISE_ENABLED
91 : #define CFC_TYPE cfcdise_type
92 : #define getErrChainRead getErrChainReadDISE
93 : #define getErrChainWrite getErrChainWriteDISE
94 : #define chainFileColName chainFileColNameDISE
95 : #define isFailedChainResize isFailedChainResizeDISE
96 : #elif ParaDRAM_ENABLED
97 : #define CFC_TYPE cfcdram_type
98 : #define getErrChainRead getErrChainReadDRAM
99 : #define getErrChainWrite getErrChainWriteDRAM
100 : #define chainFileColName chainFileColNameDRAM
101 : #define isFailedChainResize isFailedChainResizeDRAM
102 : #elif ParaNest_ENABLED
103 : #define CFC_TYPE cfcnest_type
104 : #define getErrChainRead getErrChainReadNest
105 : #define getErrChainWrite getErrChainWriteNest
106 : #define chainFileColName chainFileColNameNest
107 : #define isFailedChainResize isFailedChainResizeNest
108 : #else
109 : #error "Unrecognized interface."
110 : #endif
111 : !> \brief
112 : !> Generate and return an object of type [err_type](@ref pm_err:err_type) whose components inform
113 : !> the occurrence of an during reading and returning the contents of the specified input chain `file` of a [ParaMonte sampler](@ref pm_sampling).<br>
114 : !>
115 : !> \details
116 : !> This generic interface has two primary calling formats:<br>
117 : !> <ol>
118 : !> <li> One in which all chain file contents are stored implicitly in the `chain` argument.
119 : !> <li> One in which all chain file contents are stored explicitly in the respectively arguments following `chain`.
120 : !> </ol>
121 : !> The latter exists solely because of the lack of PDT support in \gfortran{13.1} and earlier versions.<br>
122 : !>
123 : !> \param[inout] cfc : The input/output scalar object that can be of,<br>
124 : !> <ol>
125 : !> <li> type [cfcdise_type](@ref pm_sampling_dise::cfcdise_type) indicating the chain to read is an output of the ParaDRAM sampler,<br>
126 : !> <li> type [cfcdram_type](@ref pm_sampling_dram::cfcdram_type) indicating the chain to read is an output of the ParaDISE sampler,<br>
127 : !> <li> type [cfcnest_type](@ref pm_sampling_nest::cfcnest_type) indicating the chain to read is an output of the ParaNest sampler,<br>
128 : !> </ol>
129 : !> If any optional arguments below are present, then `chain` has `intent(in)` and will remain untouched.<br>
130 : !> If all optional arguments below are missing, then `chain` has `intent(inout)` and it components will be filled with chain file contents upon return.<br>
131 : !> \param[in] file : The input scalar `character` of default kind \SK, containing the path to the chain file whose contents is to be read.<br>
132 : !> \param[in] form : The input scalar `character` of default kind \SK, containing the `form` of the chain file whose contents is to be read.<br>
133 : !> This argument is typically the value of `outputChainFileFormat` in the [base settings](@ref pm_sampling_base) of a ParaMonte sampler.<br>
134 : !> But in general, the following values are accepted and recognized:<br>
135 : !> <ol>
136 : !> <li> The string `"unformatted"`, implying that the chain file form is `unformatted` (i.e., binary, non-humane-readable) whose contents must be read with `stream` access.<br>
137 : !> <li> The string `"formatted"`, implying that the chain file form is `formatted` (i.e., ASCII, humane-readable).<br>
138 : !> <li> The string `"compact"`, implying that the chain file form is `formatted` (i.e., ASCII, humane-readable).<br>
139 : !> Additionally, it implies that the individual states within the chain can be potentially unequally weighted.<br>
140 : !> This option is solely relevant to MCMC type of chain files (returned by [ParaDISE](pm_sampling::paradise_type) and [ParaDRAM](pm_sampling::paradram_type) samplers.<br>
141 : !> For all other samplers, it is functionally equivalent to specifying `"formatted"` for the input argument `form`.<br>
142 : !> <li> The string `"verbose"`, implying that the chain file form is `formatted` (i.e., ASCII, humane-readable).<br>
143 : !> Additionally, any two consecutive states in the chain file that are identical will be merged and the corresponding `sampleWeight` will be incremented by one.<br>
144 : !> This option is solely relevant to MCMC type of chain files (returned by [ParaDISE](pm_sampling::paradise_type) and [ParaDRAM](pm_sampling::paradram_type) samplers.<br>
145 : !> For all other samplers, it is functionally equivalent to specifying `"formatted"` for the input argument `form`.<br>
146 : !> </ol>
147 : !> The specified value for `form` must **not** be `"unformatted"` (binary) if `chain` has `intent(inout)`.<br>
148 : !> Furthermore, if the contents of the chain to read are in binary format, the kinds of the column storage space vectors,
149 : !> particularly, for columns holding `real` value must be correct, otherwise, the IO process is bound to fail, as the procedure has no way of discerning the kinds of the values from the file.<br>
150 : !> \param[in] pre : The input scalar `integer` of default kind \IK, containing the pre-specified chain size with which all output components used for storing the chain file contents must be allocated.<br>
151 : !> This is relevant to the case where further samples are to be added to the components, for example, in a simulation restart.<br>
152 : !> Specifying `pre` will also lead to faster chain read action, because it requires only one memory allocation per chain file column as opposed to repeated resizing.<br>
153 : !> However, if specified, `pre` must be at least as large as the total number rows in the chain file (minus `1`, the header line).<br>
154 : !> (**optional**. If missing, data will be written to the corresponding component of `chain`. It can be present **only if** the specified `chain` object contains component with the same name.)
155 : !> \param[out] pos : The output `allocatable` vector of the same size as the output (weighted) chain (that is `nsam`) of type `integer` of default **processor** kind (which may not necessarily be the same as \IK),
156 : !> containing the positions of the beginning of the row records of the binary chain file as returned by the `pos` specifier of the Fortran intrinsic `read()` statement.<br>
157 : !> Note that `pos` does not contain the starting position of the file header as it is `1` by definition.<br>
158 : !> This vector can be later used for repositioning the chain file for subsequent writes to the file.<br>
159 : !> (**optional**. Its presence is relevant only if the input `form` is `"binary"` or `"unformatted"`.)
160 : !>
161 : !> \return
162 : !> `err` : The output scalar of type [err_type](@ref pm_err::err_type) whose component `occur` is set to the `.true.` **if and only if** the
163 : !> procedure fails to accomplish the task, in which case, the `msg` component of `err` is set to a description of the error origin.<br>
164 : !> On output, the `msg` component of `err` is an empty string if the task is accomplished successfully.<br>
165 : !> Possible causes of read action failure include (but are not limited to):<br>
166 : !> <ol>
167 : !> <li> The file header contents have been changed.
168 : !> <li> The file is read with a wrong specified `form`.
169 : !> <li> The file contents have been compromised by the user.
170 : !> </ol>
171 : !>
172 : !> \interface{getErrChainRead}
173 : !> \code{.F90}
174 : !>
175 : !> use pm_sampling_scio, only: getErrChainRead
176 : !> use pm_err, only: err_type
177 : !> type(err_type) :: err
178 : !>
179 : !> err = getErrChainRead(chain, file, form, pre = pre, pos = pos)
180 : !> !
181 : !> \endcode
182 : !>
183 : !> \impure
184 : !>
185 : !> \devnote
186 : !> Upon entry, it file is closed already, it will be opened and rewound and will be closed upon return.<br>
187 : !> If file is open, it is closed and reopened and rewound, however, it will remain open and repositioned correctly such that further write to the file can resume upon return.<br>
188 : !>
189 : !> \example{getErrChainRead}
190 : !> \include{lineno} example/pm_sampling_scio/getErrChainRead/main.F90
191 : !> \compilef{getErrChainRead}
192 : !> \output{getErrChainRead}
193 : !> \include{lineno} example/pm_sampling_scio/getErrChainRead/main.out.F90
194 : !>
195 : !> \test
196 : !> [test_pm_sampling_scio](@ref test_pm_sampling_scio)
197 : !>
198 : !> \todo
199 : !> \pvhigh
200 : !> The madness seen in the current interface of the procedures of this generic interface is due to the lack of PDT support in \gfortran{13.1}.<br>
201 : !> Ideally, the derived type [chainParaDRAM_type](@ref pm_sampling_scio::chainParaDRAM_type) must be converted to a parameterized derived type (PDT) once supported by \gfortran.<br>
202 : !> Consequently, all chain contents components that are explicitly passed to the procedures can be implicitly returned as components of the output `chain` argument.<br>
203 : !>
204 : !> \finmain{getErrChainRead}
205 : !>
206 : !> \author
207 : !> \AmirShahmoradi, September 1, 2012, 12:00 AM, National Institute for Fusion Studies, The University of Texas at Austin
208 0 : function getErrChainRead(cfc, file, form, pre, pos) result(err)
209 : #if __INTEL_COMPILER && DLL_ENABLED && (_WIN32 || _WIN64)
210 : !DEC$ ATTRIBUTES DLLEXPORT :: getErrChainRead
211 : #endif
212 : use pm_io, only: setContentsFrom
213 : use pm_strASCII, only: getStrLower
214 : !use pm_arrayRebind, only: setRebound
215 : use pm_arrayResize, only: setResized
216 : use pm_arrayFind, only: getCountLoc
217 : use pm_arraySplit, only: setSplit
218 : use pm_arrayFind, only: discrete
219 : use pm_io, only: getCountRecord
220 : use pm_io, only: setRecordFrom
221 : use pm_sysPath, only: isExtant
222 : use pm_val2str, only: getStr
223 : use pm_io, only: isOpen
224 :
225 : type(CFC_TYPE), intent(inout) :: cfc
226 : character(*, SK), intent(in) :: file, form
227 : integer, intent(out), allocatable, optional :: pos(:)
228 : integer(IK), intent(in), optional :: pre
229 0 : type(CFC_TYPE) :: row
230 : type(err_type) :: err
231 :
232 : integer :: unit, posCurrent
233 0 : character(:, SK), allocatable :: record
234 : character(*, SK), parameter :: PROCEDURE_NAME = MODULE_NAME//SK_"@getErrChainRead()"
235 : logical(LK) :: opened, isBinary, isASCII, isVerbose
236 : integer(IK) :: lenrec, ibeg, iend, ndim, pre_def, nsep, seplen
237 : #if ParaDISE_ENABLED || ParaDRAM_ENABLED
238 : logical(LK) :: isPastFirstRead
239 : isPastFirstRead = .false._LK
240 : #elif !ParaNest_ENABLED
241 : #error "Unrecognized interface."
242 : #endif
243 0 : record = getStrLower(form)
244 0 : isVerbose = record == SK_"verbose"
245 0 : isBinary = record == SK_"unformatted" .or. record == SK_"binary"
246 0 : isASCII = record == SK_"formatted" .or. record == SK_"compact" .or. record == SK_"ascii" .or. isVerbose
247 0 : err = err_type(msg = repeat(SK_" ", 255))
248 :
249 : ! First check if file is open and close it.
250 :
251 0 : inquire(file = file, exist = err%occurred, opened = opened, number = unit, iostat = err%stat, iomsg = err%msg)
252 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0 .or. .not. err%occurred,err%msg)
253 : !if (opened) close(unit, iostat = err%stat, iomsg = err%msg)
254 : !RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
255 0 : err%occurred = .false._LK ! THIS MUST BE HERE.
256 :
257 : ! open file with the right access and read the header.
258 :
259 0 : if (.not. opened) then
260 0 : if (isBinary) then
261 0 : open(newunit = unit, file = file, form = "unformatted", access = "stream", status = "old", iostat = err%stat, iomsg = err%msg SHARED)
262 0 : elseif (isASCII) then
263 0 : open(newunit = unit, file = file, status = "old", iostat = err%stat, iomsg = err%msg SHARED)
264 : end if
265 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
266 : end if
267 :
268 0 : if (isBinary) then
269 : ! header is assumed to terminate with a null character in unformatted files.
270 0 : call setResized(record, 132_IK)
271 0 : lenrec = 0_IK
272 : do
273 0 : lenrec = lenrec + 1_IK
274 0 : if (len(record, IK) < lenrec) call setResized(record)
275 0 : read(unit, pos = lenrec, iostat = err%stat, iomsg = err%msg) record(lenrec : lenrec)
276 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
277 0 : if (record(lenrec : lenrec) /= achar(0)) cycle
278 0 : lenrec = lenrec - 1_IK
279 0 : exit
280 : end do
281 0 : elseif (isASCII) then
282 0 : call setRecordFrom(unit, record, iostat = err%stat, iomsg = err%msg, ub = lenrec)
283 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
284 : else
285 0 : err%occurred = .true._LK
286 0 : err%msg = PROCEDURE_NAME//SK_"Unrecognized file `form`: """//form//SK_""""
287 0 : return
288 : end if
289 :
290 : ! Parse the header to infer the separator and the column names.
291 :
292 : ibeg = len_trim(chainFileColName(1)) + 1 ! beginning of the separator.
293 0 : seplen = index(record(ibeg : lenrec), trim(chainFileColName(2)), kind = IK) - 1_IK
294 0 : cfc%sep = record(ibeg : ibeg + seplen - 1)
295 0 : cfc%header = record(1 : lenrec)
296 0 : RETURN_IF_FAILED(__LINE__,seplen < 1_IK,SK_"Failed to infer the chain file column separator: """//cfc%header//SK_"""")
297 0 : nsep = getCountLoc(cfc%header, cfc%sep, discrete, blindness = len(cfc%sep, IK))
298 0 : ndim = nsep - size(chainFileColName, 1, IK) + 1_IK
299 0 : RETURN_IF_FAILED(__LINE__,ndim < 1_IK,SK_"Failed to infer the chain file column separator: """//cfc%header//SK_"""")
300 :
301 : ! \todo
302 : ! Additional check could be performed here to ensure the first default column names are indeed correctly inferred from the chain file.
303 :
304 : ! Allocate components.
305 :
306 0 : pre_def = 100000_IK
307 0 : if (present(pre)) pre_def = pre
308 0 : RETURN_IF_FAILED(__LINE__,isFailedChainResize(row, ndim, 1_IK, err%msg),SK_"Failed to allocate storage for the chain row.")
309 0 : RETURN_IF_FAILED(__LINE__,isFailedChainResize(cfc, ndim, pre_def, err%msg),SK_"Failed to allocate storage for the chain columns.")
310 0 : if (present(pos)) then
311 0 : call setResized(pos, pre_def, failed = err%occurred, errmsg = err%msg)
312 0 : RETURN_IF_FAILED(__LINE__,err%occurred,err%msg)
313 : end if
314 :
315 0 : cfc%nsam = 0_IK
316 : loopOverRows: do
317 :
318 0 : blockFileForm: if (isBinary) then
319 : ! \todo
320 : ! The following reading approach cannot capture end of record errors in binary mode.
321 : ! If the record is incomplete, an end of file error will be returned instead.
322 : ! As such, restart with compromised binary chain file may fail.
323 : ! This could be fixed by separating the individual field reads from each other and
324 : ! setting end of file error for all fields (except the last) to `iostat_eor` for later correction outside the read block.
325 : ! This was currently deemed too much work with little obvious gain. The future developer can decide to proceed with this or not if needed.
326 : ! There may be a remote possibility of natural incomplete record to occur under natural premature simulation stop (without manual chain file compromise).
327 : ! But the odds are extremely low and the effort to check for this not worth currently.
328 : ! An alternative solution is to always ignore the last chain row read which is followed by an end of file record.
329 : ! But this has to be done within the calling routine and not here to guarantee the output of this procedure is consistent for all chain form types.
330 : ! Update Nov 2023: After extensive tests, the current approach appears to work even for compromised binary chain files.
331 : #if ParaDRAM_ENABLED || ParaDISE_ENABLED
332 0 : inquire(unit = unit, pos = posCurrent, iostat = err%stat, iomsg = err%msg)
333 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
334 : read(unit, iostat = err%stat, iomsg = err%msg ) row%processID (1) & ! LCOV_EXCL_LINE
335 : , row%delayedRejectionStage (1) & ! LCOV_EXCL_LINE
336 : , row%meanAcceptanceRate (1) & ! LCOV_EXCL_LINE
337 : , row%adaptationMeasure (1) & ! LCOV_EXCL_LINE
338 : , row%burninLocation (1) & ! LCOV_EXCL_LINE
339 : , row%sampleWeight (1) & ! LCOV_EXCL_LINE
340 : , row%sampleLogFunc (1) & ! LCOV_EXCL_LINE
341 0 : , row%sampleState (1:ndim,1)
342 : #elif ParaNest_ENABLED
343 : read(unit, iostat = err%stat, iomsg = err%msg ) row%processID (1) & ! LCOV_EXCL_LINE
344 : , row%meanAcceptanceRate (1) & ! LCOV_EXCL_LINE
345 : , row%domainLogVol (1) & ! LCOV_EXCL_LINE
346 : , row%funcLogIntegral (1) & ! LCOV_EXCL_LINE
347 : , row%logMaxRelErr (1) & ! LCOV_EXCL_LINE
348 : , row%sampleLogWeight (1) & ! LCOV_EXCL_LINE
349 : , row%sampleLogFunc (1) & ! LCOV_EXCL_LINE
350 0 : , row%sampleState (1:ndim,1)
351 : #else
352 : #error "Unrecognized interface."
353 : #endif
354 : else blockFileForm ! only if compact or verbose.
355 : ! Read the row as a string record.
356 : ! We do so, because we do not know a priori what sort of `sep` we are dealing with.
357 0 : call setRecordFrom(unit, record, iostat = err%stat, iomsg = err%msg, ub = lenrec)
358 0 : blockParseRecord: if (err%stat == 0) then
359 : ! ibeg always refers to the position of the beginning of a column.
360 : iend = 1_IK - seplen
361 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%processID ; if (iserf(err%stat)) exit blockParseRecord
362 : #if ParaDRAM_ENABLED || ParaDISE_ENABLED
363 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%delayedRejectionStage; if (iserf(err%stat)) exit blockParseRecord
364 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%meanAcceptanceRate ; if (iserf(err%stat)) exit blockParseRecord
365 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%adaptationMeasure ; if (iserf(err%stat)) exit blockParseRecord
366 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%burninLocation ; if (iserf(err%stat)) exit blockParseRecord
367 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%sampleWeight ; if (iserf(err%stat)) exit blockParseRecord
368 : #elif ParaNest_ENABLED
369 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%meanAcceptanceRate ; if (iserf(err%stat)) exit blockParseRecord
370 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%domainLogVol ; if (iserf(err%stat)) exit blockParseRecord
371 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%funcLogIntegral ; if (iserf(err%stat)) exit blockParseRecord
372 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%logMaxRelErr ; if (iserf(err%stat)) exit blockParseRecord
373 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%sampleLogWeight ; if (iserf(err%stat)) exit blockParseRecord
374 : #else
375 : #error "Unrecognized interface."
376 : #endif
377 0 : ibeg = iend + seplen; iend = ibeg + index(record(ibeg:lenrec), cfc%sep, kind = IK) - 1; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%sampleLogFunc ; if (iserf(err%stat)) exit blockParseRecord
378 0 : ibeg = iend + seplen; iend = lenrec; read(record(ibeg:iend), *, iostat = err%stat, iomsg = err%msg) row%sampleState ; if (iserf(err%stat)) exit blockParseRecord
379 : end if blockParseRecord
380 : end if blockFileForm
381 :
382 0 : blockParseSuccess: if (err%stat == 0) then
383 0 : cfc%nsam = cfc%nsam + 1_IK
384 0 : if (pre_def < cfc%nsam) then
385 0 : pre_def = pre_def * 2_IK
386 0 : RETURN_IF_FAILED(__LINE__,present(pre),SK_"The chain file contains more rows ("//getStr(cfc%nsam)//SK_") than the specified input `pre` ("//getStr(pre_def/2)//SK_"). If you are unsure of the right size, do not specify the optional input argument `pre`.")
387 0 : RETURN_IF_FAILED(__LINE__,isFailedChainResize(cfc, ndim, pre_def, err%msg),SK_"Failed to allocate storage for the chain columns.")
388 0 : if (present(pos)) then
389 0 : call setResized(pos, pre_def, failed = err%occurred, errmsg = err%msg)
390 0 : RETURN_IF_FAILED(__LINE__,err%occurred,err%msg)
391 : end if
392 : end if
393 0 : if (present(pos)) pos(cfc%nsam) = posCurrent
394 : #if ParaDRAM_ENABLED || ParaDISE_ENABLED
395 0 : if (isVerbose) then
396 0 : if (isPastFirstRead) then
397 0 : if (1_IK < cfc%nsam) then
398 0 : if (all(row%sampleState(1 : ndim, 1) == cfc%sampleState(1 : ndim, cfc%nsam - 1))) then
399 0 : row%adaptationMeasure(1) = max(row%adaptationMeasure(1), cfc%adaptationMeasure(cfc%nsam - 1))
400 0 : row%sampleWeight(1) = row%sampleWeight(1) + cfc%sampleWeight(cfc%nsam - 1)
401 0 : cfc%nsam = cfc%nsam - 1_IK
402 : end if
403 : end if
404 : else ! happens once only after reading the first line (after header) in the chain file.
405 : isPastFirstRead = .true._LK
406 : end if
407 : end if
408 0 : cfc%processID (cfc%nsam) = row%processID (1)
409 0 : cfc%delayedRejectionStage (cfc%nsam) = row%delayedRejectionStage (1)
410 0 : cfc%meanAcceptanceRate (cfc%nsam) = row%meanAcceptanceRate (1)
411 0 : cfc%adaptationMeasure (cfc%nsam) = row%adaptationMeasure (1)
412 0 : cfc%burninLocation (cfc%nsam) = row%burninLocation (1)
413 0 : cfc%sampleWeight (cfc%nsam) = row%sampleWeight (1)
414 0 : cfc%sampleLogFunc (cfc%nsam) = row%sampleLogFunc (1)
415 0 : cfc%sampleState (1 : ndim, cfc%nsam) = row%sampleState (1:ndim, 1)
416 : #elif ParaNest_ENABLED
417 0 : cfc%processID (cfc%nsam) = row%processID (1)
418 0 : cfc%meanAcceptanceRate (cfc%nsam) = row%meanAcceptanceRate (1)
419 0 : cfc%domainLogVol (cfc%nsam) = row%domainLogVol (1)
420 0 : cfc%funcLogIntegral (cfc%nsam) = row%funcLogIntegral (1)
421 0 : cfc%logMaxRelErr (cfc%nsam) = row%logMaxRelErr (1)
422 0 : cfc%sampleLogWeight (cfc%nsam) = row%sampleLogWeight (1)
423 0 : cfc%sampleLogFunc (cfc%nsam) = row%sampleLogFunc (1)
424 0 : cfc%sampleState (1 : ndim, cfc%nsam) = row%sampleState (1:ndim, 1)
425 : #else
426 : #error "Unrecognized interface."
427 : #endif
428 : cycle loopOverRows
429 : end if blockParseSuccess
430 :
431 : exit loopOverRows
432 :
433 : end do loopOverRows
434 :
435 : ! last read is always end-of-file (normal) or end-of-record error (corrupt file).
436 :
437 0 : RETURN_IF_FAILED(__LINE__,.not. iserf(err%stat),err%msg)
438 :
439 : ! End of record condition requires ignoring the last read sample. This gets tricky for verbose chain where the IO could fail amid writing the same unweighted state repeatedly.
440 :
441 0 : if (is_iostat_eor(err%stat)) then
442 0 : err%iswarned = .true._LK
443 : err%msg = SKC_"An end-of-record condition occurred while parsing the contents of the chain file at chain row "//getStr(cfc%nsam)//SKC_": "//trim(err%msg)//& ! LCOV_EXCL_LINE
444 0 : SKC_". Assuming the previous chain sampleState as the last in the chain file with the corresponding values: "
445 0 : if (0 == cfc%nsam) then
446 0 : err%msg = err%msg//SKC_"Nothing. The first row of the chain file is corrupt (incomplete)."
447 : else
448 0 : err%msg = err%msg//getStr(cfc%sampleState(1 : ndim, cfc%nsam - 1))
449 : ! return the position to the end of the last successful row read.
450 0 : if (opened) then
451 0 : if (isBinary) then
452 0 : read(unit, pos = posCurrent, iostat = err%stat, iomsg = err%msg)
453 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
454 : else
455 0 : backspace(unit, iostat = err%stat, iomsg = err%msg)
456 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
457 : #if ParaDRAM_ENABLED || ParaDISE_ENABLED
458 0 : if (isVerbose) then
459 0 : do ibeg = 1, cfc%sampleWeight(cfc%nsam)
460 0 : backspace(unit, iostat = err%stat, iomsg = err%msg)
461 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
462 : end do
463 : end if
464 : #elif !ParaNest_ENABLED
465 : #error "Unrecognized interface."
466 : #endif
467 : end if
468 : end if
469 : end if
470 0 : elseif (is_iostat_end(err%stat)) then
471 0 : if (isASCII) then
472 : ! End of file condition for formatted files requires a backspace to go before the end of file record.
473 0 : backspace(unit, iostat = err%stat, iomsg = err%msg)
474 0 : elseif (isBinary) then
475 : ! Ignore the last read row before hitting end of file, as the record may have been incomplete.
476 0 : read(unit, pos = posCurrent, iostat = err%stat, iomsg = err%msg)
477 : end if
478 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
479 : end if
480 :
481 : ! Clean up the output.
482 :
483 0 : if (.not. present(pre)) then
484 0 : RETURN_IF_FAILED(__LINE__,isFailedChainResize(cfc, ndim, cfc%nsam, err%msg),SK_"Failed to resize the storage for the chain columns.")
485 0 : if (present(pos)) then
486 0 : call setResized(pos, cfc%nsam, failed = err%occurred, errmsg = err%msg)
487 0 : RETURN_IF_FAILED(__LINE__,err%occurred,err%msg)
488 : end if
489 : end if
490 :
491 0 : if (.not. opened) then
492 0 : close(unit, iostat = err%stat, iomsg = err%msg)
493 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
494 : end if
495 :
496 0 : end function
497 :
498 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
499 :
500 : ! \todo
501 : ! The performance of the generic interface might slightly (perhaps negligibly) improve
502 : ! if the array initializations are only done for the left-over elements when the input argument `pre` is specified.
503 : ! For now, code clarity by avoiding code duplication wins over the slight performance improvement.
504 13 : function isFailedChainResize(cfc, ndim, nrow, errmsg) result(failed)
505 : use pm_arrayRefill, only: setRefilled
506 : real(RKC), parameter :: NULL_RK = -huge(0._RKC)
507 : integer(IK), parameter :: NULL_IK = -huge(0_IK)
508 : type(CFC_TYPE), intent(inout) :: cfc
509 : integer(IK), intent(in) :: ndim, nrow
510 : character(*, SK), intent(out) :: errmsg
511 : logical(LK) :: failed
512 13 : call setRefilled(cfc%processID , NULL_IK , nrow , failed = failed, errmsg = errmsg)
513 13 : call setRefilled(cfc%meanAcceptanceRate , NULL_RK , nrow , failed = failed, errmsg = errmsg)
514 : #if ParaDISE_ENABLED || ParaDRAM_ENABLED
515 13 : call setRefilled(cfc%delayedRejectionStage , NULL_IK , nrow , failed = failed, errmsg = errmsg)
516 13 : call setRefilled(cfc%adaptationMeasure , NULL_RK , nrow , failed = failed, errmsg = errmsg)
517 13 : call setRefilled(cfc%burninLocation , NULL_IK , nrow , failed = failed, errmsg = errmsg)
518 13 : call setRefilled(cfc%sampleWeight , NULL_IK , nrow , failed = failed, errmsg = errmsg)
519 : #elif ParaNest_ENABLED
520 0 : call setRefilled(cfc%domainLogVol , NULL_RK , nrow , failed = failed, errmsg = errmsg)
521 0 : call setRefilled(cfc%funcLogIntegral , NULL_RK , nrow , failed = failed, errmsg = errmsg)
522 0 : call setRefilled(cfc%logMaxRelErr , NULL_RK , nrow , failed = failed, errmsg = errmsg)
523 0 : call setRefilled(cfc%sampleLogWeight , NULL_RK , nrow , failed = failed, errmsg = errmsg)
524 : #else
525 : #error "Unrecognized interface."
526 : #endif
527 13 : call setRefilled(cfc%sampleLogFunc , NULL_RK , nrow , failed = failed, errmsg = errmsg)
528 39 : call setRefilled(cfc%sampleState , NULL_RK , [ndim, nrow] , failed = failed, errmsg = errmsg)
529 13 : end function
530 :
531 : !> \brief
532 : !> Write the chain properties to the chain file.
533 : !>
534 : !> \param[in] cfc : The object of class [ChainFileContents_type](@ref ChainFileContents_type).
535 : !> \param[in] ibeg : The beginning index of the compact chain beyond which the elements of the chain will be written to the output file.
536 : !> \param[in] iend : The ending index of the compact chain below which the elements of the chain will be written to the output file.
537 : !> \param[in] unit : The unit ID of the chain file to which the header should be written.
538 : !> \param[in] form : The file format of the chain file (`"binary"` vs. `"compact"` vs. `"verbose"`).
539 : !> \param[in] format : The Fortran IO formatting string to be used to read the contents of the chain file (**optional**).
540 : !> This argument is only required with a non-binary chain file, i.e., when `isBinary = .false.`.
541 : !> \param[in] proposalAdaptationPeriod : The adaptive update period (**optional**). It must be provided if `form = "verbose"`.
542 : !>
543 : !> \warning
544 : !> The file header is written only if the input is not already open.
545 : !>
546 : !> \devnote
547 : !> Upon entry, it file is closed already, it will be opened and rewound and will be closed upon return.
548 : !> If file is open, it is closed and reopened and rewound, however, it will remain open at the last position upon return.
549 0 : function getErrChainWrite(cfc, file, form, format, ibeg, iend, proposalAdaptationPeriod) result(err)
550 : #if __INTEL_COMPILER && DLL_ENABLED && (_WIN32 || _WIN64)
551 : !DEC$ ATTRIBUTES DLLEXPORT :: getErrChainWrite
552 : #endif
553 : use pm_strASCII, only: setStrLower
554 : type(CFC_TYPE) , intent(in) :: cfc
555 : character(*, SK), intent(in) :: file, form
556 : integer(IK) , intent(in), optional :: ibeg, iend
557 : character(*, SK), intent(in), optional :: format
558 : integer(IK) , intent(in), optional :: proposalAdaptationPeriod
559 : type(err_type) :: err
560 :
561 : character(*, SK), parameter :: PROCEDURE_NAME = MODULE_NAME//SK_"@getErrChainWrite()"
562 : logical(LK) :: opened, isBinary, isCompact, isVerbose, isASCII!, exist
563 : integer(IK) :: i, ibeg_def, iend_def
564 0 : character((len(form, IK)), SK) :: lcform
565 : integer :: unit
566 :
567 : ibeg_def = 1_IK
568 0 : iend_def = cfc%nsam
569 0 : if (present(ibeg)) ibeg_def = ibeg
570 0 : if (present(iend)) iend_def = iend
571 0 : RETURN_IF_FAILED(__LINE__,ibeg_def < 1_IK .or. iend_def < ibeg_def,SK_": The condition `0 < ibeg <= iend <= cfc%nsam` must hold. ibeg, iend, cfc%nsam = "//getStr([ibeg_def, iend_def, cfc%nsam]))
572 :
573 : isBinary = .false._LK
574 : isCompact = .false._LK
575 : isVerbose = .false._LK
576 0 : lcform = form
577 0 : call setStrLower(lcform)
578 0 : isVerbose = lcform == SK_"verbose"
579 0 : isBinary = lcform == SK_"unformatted" .or. lcform == SK_"binary"
580 0 : isASCII = lcform == SK_"formatted" .or. lcform == SK_"compact" .or. isVerbose
581 0 : RETURN_IF_FAILED(__LINE__,isASCII .and. .not. present(format),SK_": The optional argument `format` must be present when the file is formatted. form = """//form//SK_"""")
582 0 : RETURN_IF_FAILED(__LINE__,isVerbose .and. .not. present(proposalAdaptationPeriod),SK_": The optional argument `proposalAdaptationPeriod` must be present when the file is verbose. form = """//form//SK_"""")
583 : #if ParaNest_ENABLED
584 : isCompact = isCompact .or. isVerbose
585 : isVerbose = .false._LK
586 : #endif
587 0 : err%msg = repeat(SKC_" ", 255)
588 0 : inquire(file = file, opened = opened, iostat = err%stat, iomsg = err%msg)!, exist = exist
589 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)!.not. exist .or.
590 :
591 0 : if (.not. opened) then
592 0 : if (isBinary) then
593 0 : open(newunit = unit, file = file, form = "unformatted", access = "stream", iostat = err%stat, iomsg = err%msg SHARED)
594 : else
595 0 : open(newunit = unit, file = file, form = "formatted", access = "sequential", iostat = err%stat, iomsg = err%msg SHARED)
596 : end if
597 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
598 0 : if (isBinary) then
599 0 : write(unit, iostat = err%stat, iomsg = err%msg) cfc%header//achar(0)
600 : else
601 0 : write(unit, "(A)", iostat = err%stat, iomsg = err%msg) cfc%header
602 : end if
603 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
604 : end if
605 :
606 0 : blockWrite: if (isBinary) then
607 0 : do i = ibeg_def, iend_def
608 0 : write(unit, iostat = err%stat, iomsg = err%msg ) cfc%processID(i) &
609 : #if ParaDISE_ENABLED || ParaDRAM_ENABLED
610 0 : , cfc%delayedRejectionStage(i) &
611 0 : , cfc%meanAcceptanceRate(i) &
612 0 : , cfc%adaptationMeasure(i) &
613 0 : , cfc%burninLocation(i) &
614 0 : , cfc%sampleWeight(i) &
615 : #elif ParaNest_ENABLED
616 0 : , cfc%meanAcceptanceRate(i) &
617 0 : , cfc%domainLogVol(i) &
618 0 : , cfc%funcLogIntegral(i) &
619 0 : , cfc%sampleLogWeight(i) &
620 : #else
621 : #error "Unrecognized interface."
622 : #endif
623 0 : , cfc%sampleLogFunc(i) &
624 0 : , cfc%sampleState(:, i)
625 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
626 : end do
627 0 : elseif (isASCII) then blockWrite
628 0 : do i = ibeg_def, iend_def
629 0 : write(unit, format, iostat = err%stat, iomsg = err%msg ) cfc%processID(i) &
630 : #if ParaDISE_ENABLED || ParaDRAM_ENABLED
631 0 : , cfc%delayedRejectionStage(i) &
632 0 : , cfc%meanAcceptanceRate(i) &
633 0 : , cfc%adaptationMeasure(i) &
634 0 : , cfc%burninLocation(i) &
635 0 : , cfc%sampleWeight(i) &
636 : #elif ParaNest_ENABLED
637 0 : , cfc%meanAcceptanceRate(i) &
638 0 : , cfc%domainLogVol(i) &
639 0 : , cfc%funcLogIntegral(i) &
640 0 : , cfc%sampleLogWeight(i) &
641 : #else
642 : #error "Unrecognized interface."
643 : #endif
644 0 : , cfc%sampleLogFunc(i) &
645 0 : , cfc%sampleState(:, i)
646 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
647 : end do
648 : else
649 : #if ParaDISE_ENABLED || ParaDRAM_ENABLED
650 0 : if (isVerbose) then
651 : block
652 : integer(IK) :: j, counter
653 : real(RKC) :: adaptation
654 : counter = ibeg_def
655 0 : do i = ibeg_def, iend_def
656 0 : do j = 1, cfc%sampleWeight(i)
657 0 : if (mod(counter, proposalAdaptationPeriod) == 0_IK) then
658 0 : adaptation = cfc%adaptationMeasure(i)
659 : else
660 0 : adaptation = 0._RKC
661 : end if
662 0 : write(unit, format, iostat = err%stat, iomsg = err%msg ) cfc%processID(i) &
663 0 : , cfc%delayedRejectionStage(i) &
664 0 : , cfc%meanAcceptanceRate(i) &
665 0 : , cfc%adaptationMeasure(i) &
666 0 : , cfc%burninLocation(i) &
667 0 : , 1_IK &
668 0 : , cfc%sampleLogFunc(i) &
669 0 : , cfc%sampleState(:, i)
670 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
671 0 : counter = counter + 1
672 : end do
673 : end do
674 : exit blockWrite
675 : end block
676 : end if
677 : #endif
678 0 : RETURN_IF_FAILED(__LINE__,.true.,SK_": Unrecognized input value for the input argument `form`. form = """//form//SK_"""")
679 : end if blockWrite
680 :
681 0 : flush(unit)
682 0 : if (.not. opened) close(unit, iostat = err%stat, iomsg = err%msg)
683 0 : RETURN_IF_FAILED(__LINE__,err%stat /= 0,err%msg)
684 :
685 0 : end function getErrChainWrite
686 :
687 : #undef isFailedChainResize
688 : #undef RETURN_IF_FAILED
689 : #undef chainFileColName
690 : #undef getErrChainWrite
691 : #undef getErrChainRead
692 : #undef CFC_TYPE
693 :
694 : #else
695 : !%%%%%%%%%%%%%%%%%%%%%%%%
696 : #error "Unrecognized interface."
697 : !%%%%%%%%%%%%%%%%%%%%%%%%
698 : #endif
|
