Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!!
4 : !!!! MIT License
5 : !!!!
6 : !!!! ParaMonte: plain powerful parallel Monte Carlo library.
7 : !!!!
8 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab
9 : !!!!
10 : !!!! This file is part of the ParaMonte library.
11 : !!!!
12 : !!!! Permission is hereby granted, free of charge, to any person obtaining a
13 : !!!! copy of this software and associated documentation files (the "Software"),
14 : !!!! to deal in the Software without restriction, including without limitation
15 : !!!! the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : !!!! and/or sell copies of the Software, and to permit persons to whom the
17 : !!!! Software is furnished to do so, subject to the following conditions:
18 : !!!!
19 : !!!! The above copyright notice and this permission notice shall be
20 : !!!! included in all copies or substantial portions of the Software.
21 : !!!!
22 : !!!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23 : !!!! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24 : !!!! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
25 : !!!! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
26 : !!!! DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
27 : !!!! OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
28 : !!!! OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
29 : !!!!
30 : !!!! ACKNOWLEDGMENT
31 : !!!!
32 : !!!! ParaMonte is an honor-ware and its currency is acknowledgment and citations.
33 : !!!! As per the ParaMonte library license agreement terms, if you use any parts of
34 : !!!! this library for any purposes, kindly acknowledge the use of ParaMonte in your
35 : !!!! work (education/research/industry/development/...) by citing the ParaMonte
36 : !!!! library as described on this page:
37 : !!!!
38 : !!!! https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
39 : !!!!
40 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42 :
43 : !> \brief This module contains the classes and procedures for refining MCMC output chains.
44 : !> \author Amir Shahmoradi
45 :
46 : module ParaMCMCRefinedChain_mod
47 :
48 : use SpecMCMC_SampleRefinementMethod_mod, only: BATCH_MEANS_METHOD_NAME
49 : use SpecMCMC_SampleRefinementMethod_mod, only: CUTOFF_AUTOCORR_METHOD_NAME
50 : use SpecMCMC_SampleRefinementMethod_mod, only: MAX_CUMSUM_AUTOCORR_METHOD_NAME
51 : use ParaMonteChainFileContents_mod, only: Count_type
52 : use JaggedArray_mod, only: CharVec_type
53 : use Constants_mod, only: IK, RK
54 : use Err_mod, only: Err_type
55 :
56 : implicit none
57 :
58 : character(*), parameter :: MODULE_NAME = "\paramCMCRefinedChain_mod"
59 :
60 : !> The `RefinedChain_type` class.
61 : type :: RefinedChain_type
62 : integer(IK) :: ndim = 0_IK ! number of sampling variables
63 : integer(IK) :: numRefinement = 0_IK ! number of refinements, zero if sample size is prescribed by the user
64 : type(Count_type) , allocatable :: Count(:) ! compact and verbose counts
65 : real(RK) , allocatable :: IAC(:,:) ! size of (ndim,0:numRefinement): The Integrated AutoCorrelation Time at each refinement stage
66 : real(RK) , allocatable :: LogFuncState(:,:) ! size of (ndim,Count%compact): LogFuncState is LogFunc + Variables
67 : integer(IK) , allocatable :: Weight(:) ! size of (Count%compact): Weight of each state
68 : type(CharVec_type) , allocatable :: ColHeader(:) ! refined sample column headers
69 : type(Err_type) :: Err
70 : contains
71 : procedure, pass :: get => getRefinedChain
72 : procedure, pass :: write => writeRefinedChain
73 : end type RefinedChain_type
74 :
75 : type Method_type
76 : logical :: isMax = .false.
77 : logical :: isMed = .false.
78 : logical :: isMin = .false.
79 : logical :: isAvg = .false.
80 : logical :: isBatchMeans = .false.
81 : logical :: isCutoffAutoCorr = .false.
82 : logical :: isViaCompactChain = .true.
83 : logical :: isViaVerboseChain = .true.
84 : logical :: isMaxCumSumAutoCorr = .false.
85 : end type Method_type
86 :
87 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88 :
89 : contains
90 :
91 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
92 :
93 : !> Return the refined Markov chain, given the input Markov chain and its specifications.
94 : !> This procedure is a method of the [RefinedChain_type](@ref refinedchain_type) class.
95 : !>
96 : !> \param[inout] RefinedChain : An object of class [RefinedChain_type](@ref refinedchain_type).
97 : !> \param[inout] CFC : An object of type [ChainFileContents_type](@ref paramontechainfilecontents_mod::chainfilecontents_type)
98 : !! containing the Markov chain.
99 : !> \param[out] Err : An object of class [Err_type](@ref err_mod::err_type) indicating whether any error has occurred or not.
100 : !> \param[in] burninLoc : The estimated location of burnin point in the Markov chain (**optional**).
101 : !! If not provided, it will be extracted from the components of the input `CFC`.
102 : !> \param[in] refinedChainSize : The requested refined sample size (**optional**). If the size of the refined sample is given as input,
103 : !! then the requested sample is directly generated based on the input size.
104 : !> \param[in] sampleRefinementCount : The maximum number of times the sample can be refined (**optional**, default = `Infinity`).
105 : !! : For example, if set to 1, then only one round of refinement will be performed on the Markov chain.
106 : !> \param[in] sampleRefinementMethod : The requested method of refining the sample (**optional**, default = "BatchMeans").
107 585 : subroutine getRefinedChain ( RefinedChain &
108 : , CFC &
109 : , Err &
110 : , burninLoc &
111 : , refinedChainSize &
112 : , sampleRefinementCount &
113 : , sampleRefinementMethod &
114 : )
115 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
116 : !DEC$ ATTRIBUTES DLLEXPORT :: getRefinedChain
117 : #endif
118 :
119 : use, intrinsic :: iso_fortran_env, only: output_unit
120 : use ParaMonteChainFileContents_mod, only: ChainFileContents_type, NUM_DEF_COL
121 : use CrossCorr_mod, only: getBatchMeansIAC, getCumSumIAC, getMaxCumSumIAC
122 : use String_mod, only: getLowerCase, replaceStr
123 : use Sort_mod, only: getMedian
124 :
125 : implicit none
126 :
127 : character(*), parameter :: PROCEDURE_NAME = MODULE_NAME//"@getRefinedChain()"
128 :
129 : class(RefinedChain_type) , intent(inout) :: RefinedChain
130 : type(ChainFileContents_type), intent(inout) :: CFC
131 : type(Err_type) , intent(out) :: Err
132 : character(*) , intent(in), optional :: sampleRefinementMethod
133 : integer(IK) , intent(in), optional :: burninLoc, refinedChainSize, sampleRefinementCount
134 : integer(IK) :: burninLocDefault, i, countCompact, ndimPlusOne
135 : integer(IK) :: maxRefinementCurrent, maxRefinementCount
136 354 : real(RK) :: integratedAutoCorrTime, ndimPlusOneInverse
137 : logical :: refinedChainSizeIsPresent
138 : logical :: maxRefinementCountIsReached
139 354 : character(:) , allocatable :: sampleRefinementMethodDefault, sampleRefinementMethodDefaultLowerCase, srmethod
140 : type(Count_type), allocatable :: DumCount(:)
141 : real(RK) , allocatable :: DumIAC(:,:), DumVec(:)
142 : integer(IK) , allocatable :: SampleWeight(:)
143 : type(Method_type) :: Method
144 354 : type(ChainFileContents_type) :: DumCFC
145 :
146 354 : Err%occurred = .false.
147 :
148 : ! if the size of the refined sample is given as input, then generate the requested sample straight
149 :
150 354 : refinedChainSizeIsPresent = present(refinedChainSize)
151 354 : if (refinedChainSizeIsPresent) then ! ignore sampleRefinementCount, even if it is given by the user
152 85 : maxRefinementCount = 1_IK
153 : else
154 269 : maxRefinementCount = 20_IK ! this is a temporary maximum value, to be increased later if needed
155 269 : if (present(sampleRefinementCount)) maxRefinementCount = sampleRefinementCount
156 : end if
157 :
158 : ! this is to avoid memory overflow due to extremely large maxRefinementCount requested by the user
159 :
160 354 : maxRefinementCurrent = min(2_IK, maxRefinementCount)
161 :
162 : ! compute ndim and the initial chain size
163 :
164 354 : RefinedChain%numRefinement = 0_IK
165 :
166 354 : if (CFC%ndim == 0_IK) then
167 6 : RefinedChain%ndim = size(CFC%State(:,1))
168 6 : CFC%ndim = RefinedChain%ndim
169 : else
170 348 : RefinedChain%ndim = CFC%ndim
171 : end if
172 :
173 : ! allocate components
174 :
175 354 : if (allocated(RefinedChain%IAC)) deallocate(RefinedChain%IAC)
176 354 : if (allocated(RefinedChain%Count)) deallocate(RefinedChain%Count)
177 354 : allocate(RefinedChain%IAC(0:RefinedChain%ndim,0:maxRefinementCurrent))
178 1327 : allocate(RefinedChain%Count(0:maxRefinementCurrent))
179 354 : if (CFC%Count%compact == 0_IK) then
180 6 : RefinedChain%Count(RefinedChain%numRefinement)%compact = size(CFC%State(1,:))
181 6 : CFC%Count%compact = RefinedChain%Count(RefinedChain%numRefinement)%compact
182 : else
183 348 : RefinedChain%Count(RefinedChain%numRefinement)%compact = CFC%Count%compact
184 : end if
185 354 : if (CFC%Count%verbose == 0_IK) CFC%Count%verbose = sum(CFC%Weight(1:CFC%Count%compact))
186 354 : RefinedChain%Count(RefinedChain%numRefinement)%verbose = CFC%Count%verbose
187 :
188 : ! define the AIC computation method
189 :
190 354 : sampleRefinementMethodDefault = BATCH_MEANS_METHOD_NAME; if (present(sampleRefinementMethod)) sampleRefinementMethodDefault = sampleRefinementMethod
191 354 : sampleRefinementMethodDefaultLowerCase = getLowerCase(sampleRefinementMethodDefault)
192 :
193 354 : if (index(sampleRefinementMethodDefaultLowerCase,getLowerCase(BATCH_MEANS_METHOD_NAME))>0) then
194 294 : Method%isBatchMeans = .true.
195 120 : elseif (index(sampleRefinementMethodDefaultLowerCase,getLowerCase(CUTOFF_AUTOCORR_METHOD_NAME))>0 .or. index(sampleRefinementMethodDefaultLowerCase,"cutoff")>0) then
196 28 : Method%isCutoffAutoCorr = .true.
197 64 : elseif (index(sampleRefinementMethodDefaultLowerCase,getLowerCase(MAX_CUMSUM_AUTOCORR_METHOD_NAME))>0 .or. index(sampleRefinementMethodDefaultLowerCase,"cumsum")>0) then
198 32 : Method%isMaxCumSumAutoCorr = .true.
199 : else
200 : ! LCOV_EXCL_START
201 : Err%occurred = .true.
202 : Err%msg = PROCEDURE_NAME // ": Unknown unsupported IAC computation method name: " // sampleRefinementMethodDefault
203 : return
204 : ! LCOV_EXCL_STOP
205 : end if
206 :
207 : ! define the chain types to use for the AIC computation
208 :
209 354 : Method%isViaCompactChain = index(sampleRefinementMethodDefaultLowerCase,"compact") > 0
210 354 : Method%isViaVerboseChain = index(sampleRefinementMethodDefaultLowerCase,"verbose") > 0
211 354 : if (.not.(Method%isViaCompactChain .or. Method%isViaVerboseChain)) then
212 354 : Method%isViaCompactChain = .true.
213 354 : Method%isViaVerboseChain = .true.
214 : end if
215 :
216 : ! define the statistic to use for the AIC computation
217 :
218 354 : srmethod = replaceStr(string=sampleRefinementMethodDefaultLowerCase,search=" ",substitute="-")
219 354 : Method%isAvg = index(srmethod,"-avg") > 0 .or. index(srmethod,"-average") > 0
220 354 : Method%isMed = index(srmethod,"-med") > 0 .or. index(srmethod,"-median") > 0
221 354 : Method%isMin = index(srmethod,"-min") > 0 .or. index(srmethod,"-minimum") > 0
222 354 : Method%isMax = index(srmethod,"-max") > 0 .or. index(srmethod,"-maximum") > 0
223 354 : if ( Method%isAvg .and. (Method%isMed .or. Method%isMax .or. Method%isMin) ) Err%occurred = .true.
224 354 : if ( Method%isMed .and. (Method%isMax .or. Method%isMin) ) Err%occurred = .true.
225 354 : if ( Method%isMax .and. Method%isMin ) Err%occurred = .true.
226 354 : if (.not.(Method%isAvg .or. Method%isMed .or. Method%isMin .or. Method%isMax)) Method%isAvg = .true. ! default method of AIC summarization.
227 :
228 354 : if (Method%isAvg) ndimPlusOneInverse = 1._RK / (RefinedChain%ndim + 1_IK)
229 354 : if (Method%isMed) ndimPlusOne = RefinedChain%ndim + 1_IK
230 :
231 : ! assign the column headers
232 :
233 354 : if (allocated(CFC%ColHeader)) then
234 768 : if (allocated(RefinedChain%ColHeader)) deallocate(RefinedChain%ColHeader)
235 1306 : allocate(RefinedChain%ColHeader(0:RefinedChain%ndim))
236 1306 : do i = 0, RefinedChain%ndim
237 1306 : RefinedChain%ColHeader(i)%record = CFC%ColHeader(i+NUM_DEF_COL)%record
238 : end do
239 : else
240 : ! LCOV_EXCL_START
241 : Err%occurred = .true.
242 : Err%msg = PROCEDURE_NAME // ": Internal error occurred. CFC%ColHeader is NULL."
243 : return
244 : ! LCOV_EXCL_STOP
245 : end if
246 :
247 354 : if (present(burninLoc)) then
248 348 : burninLocDefault = burninLoc
249 : else
250 6 : burninLocDefault = CFC%BurninLoc(CFC%Count%compact)
251 : end if
252 :
253 354 : RefinedChain%Count(RefinedChain%numRefinement)%compact = CFC%Count%compact - burninLocDefault + 1
254 354 : if (allocated(RefinedChain%LogFuncState)) deallocate(RefinedChain%LogFuncState)
255 354 : allocate(RefinedChain%LogFuncState(0:RefinedChain%ndim,RefinedChain%Count(RefinedChain%numRefinement)%compact))
256 139086 : RefinedChain%LogFuncState(0 ,1:RefinedChain%Count(RefinedChain%numRefinement)%compact) = CFC%LogFunc(burninLocDefault:CFC%Count%compact)
257 384932 : RefinedChain%LogFuncState(1:RefinedChain%ndim,1:RefinedChain%Count(RefinedChain%numRefinement)%compact) = CFC%State(1:RefinedChain%ndim, burninLocDefault:CFC%Count%compact)
258 :
259 : ! check if there are more than 1 sample points in the burnin-subtracted CFC
260 :
261 354 : if (RefinedChain%Count(RefinedChain%numRefinement)%compact==0_IK) then
262 : ! LCOV_EXCL_START
263 : Err%occurred = .true.
264 : Err%msg = PROCEDURE_NAME // ": The size of the refined sample is zero."
265 : return
266 : ! LCOV_EXCL_STOP
267 354 : elseif (RefinedChain%Count(RefinedChain%numRefinement)%compact==1_IK) then
268 6 : if (allocated(RefinedChain%Weight)) deallocate(RefinedChain%Weight)
269 6 : allocate(RefinedChain%Weight(RefinedChain%Count(RefinedChain%numRefinement)%compact))
270 6 : if (refinedChainSizeIsPresent) then
271 0 : RefinedChain%Weight = refinedChainSize
272 : else
273 12 : RefinedChain%Weight = 1
274 : end if
275 6 : if (allocated(RefinedChain%IAC)) deallocate(RefinedChain%IAC)
276 6 : if (allocated(RefinedChain%Count)) deallocate(RefinedChain%Count)
277 6 : allocate(RefinedChain%IAC(0:RefinedChain%ndim,0:0))
278 12 : allocate(RefinedChain%Count(0:0))
279 30 : RefinedChain%IAC = 0._RK
280 12 : RefinedChain%Count(RefinedChain%numRefinement)%verbose = sum(RefinedChain%Weight)
281 6 : return
282 : end if
283 :
284 139075 : RefinedChain%Weight = CFC%Weight(burninLocDefault:CFC%Count%compact) ! Weight is intentionally separately assigned from State here
285 :
286 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
287 :
288 515 : loopRefinement: do
289 :
290 1211 : countCompact = RefinedChain%Count(RefinedChain%numRefinement)%compact
291 :
292 : ! set the sample weight
293 :
294 1211 : if (Method%isViaCompactChain) then
295 722 : if (allocated(SampleWeight)) deallocate(SampleWeight)
296 489 : elseif (Method%isViaVerboseChain) then
297 95866 : SampleWeight = RefinedChain%Weight
298 : endif
299 :
300 : ! obtain the IAC for each individual variable
301 :
302 4470 : do i = 0, RefinedChain%ndim
303 793028 : DumVec = RefinedChain%LogFuncState(i,1:RefinedChain%Count(RefinedChain%numRefinement)%compact)
304 4470 : if (Method%isBatchMeans) then
305 10828 : RefinedChain%IAC(i,RefinedChain%numRefinement) = getBatchMeansIAC ( np = countCompact &
306 : , Point = DumVec & ! LCOV_EXCL_LINE
307 : , Weight = SampleWeight & ! LCOV_EXCL_LINE
308 2707 : )
309 552 : elseif (Method%isCutoffAutoCorr) then
310 1020 : RefinedChain%IAC(i,RefinedChain%numRefinement) = getCumSumIAC ( np = countCompact &
311 : , Point = DumVec & ! LCOV_EXCL_LINE
312 : , Weight = SampleWeight & ! LCOV_EXCL_LINE
313 255 : )
314 297 : elseif (Method%isMaxCumSumAutoCorr) then
315 1188 : RefinedChain%IAC(i,RefinedChain%numRefinement) = getMaxCumSumIAC ( np = countCompact &
316 : , Point = DumVec & ! LCOV_EXCL_LINE
317 : , Weight = SampleWeight & ! LCOV_EXCL_LINE
318 297 : )
319 : end if
320 : end do
321 :
322 1211 : if (refinedChainSizeIsPresent) then
323 83151 : integratedAutoCorrTime = real(sum(RefinedChain%Weight),kind=RK) / real(refinedChainSize,kind=RK)
324 : else
325 980 : if (Method%isAvg) then
326 2842 : integratedAutoCorrTime = sum( RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement) ) * ndimPlusOneInverse
327 179 : elseif (Method%isMed) then
328 90 : call getMedian(lenArray=ndimPlusOne,Array=RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement),median=integratedAutoCorrTime,Err=Err)
329 90 : if (Err%occurred) then
330 : ! LCOV_EXCL_START
331 : Err%msg = PROCEDURE_NAME//Err%msg
332 : return
333 : ! LCOV_EXCL_STOP
334 : end if
335 89 : elseif (Method%isMax) then
336 320 : integratedAutoCorrTime = maxval( RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement) )
337 25 : elseif (Method%isMin) then
338 125 : integratedAutoCorrTime = minval( RefinedChain%IAC(0:RefinedChain%ndim,RefinedChain%numRefinement) )
339 : end if
340 : end if
341 :
342 : ! so far, we have computed the max IAC of the sample, but no refinement. Refine the sample only if needed.
343 :
344 1211 : maxRefinementCountIsReached = RefinedChain%numRefinement==maxRefinementCount
345 1211 : if (integratedAutoCorrTime<2._RK .or. maxRefinementCountIsReached) then
346 696 : if (Method%isViaCompactChain .and. Method%isViaVerboseChain) then
347 : !if (maxRefinementCountIsReached) maxRefinementCount = maxRefinementCount * 2_IK
348 348 : maxRefinementCount = maxRefinementCount * 2_IK
349 348 : Method%isViaCompactChain = .false.
350 348 : cycle loopRefinement
351 : end if
352 348 : exit loopRefinement
353 : end if
354 :
355 : ! generate the refined sample, dump it in CFC, then put it back into RefinedChain to start over again
356 :
357 515 : RefinedChain%numRefinement = RefinedChain%numRefinement + 1_IK
358 :
359 : ! reallocate to bigger array if nedded
360 :
361 515 : if (RefinedChain%numRefinement>maxRefinementCurrent) then
362 45 : call move_alloc( from = RefinedChain%IAC, to = DumIAC )
363 45 : call move_alloc( from = RefinedChain%Count, to = DumCount )
364 45 : maxRefinementCurrent = min( maxRefinementCurrent*2 , maxRefinementCount )
365 45 : allocate( RefinedChain%IAC( 0:RefinedChain%ndim , 0:maxRefinementCurrent ) )
366 274 : allocate( RefinedChain%Count( 0:maxRefinementCurrent ) )
367 575 : RefinedChain%IAC(0:RefinedChain%ndim,0:RefinedChain%numRefinement-1) = DumIAC
368 182 : RefinedChain%Count(0:RefinedChain%numRefinement-1) = DumCount
369 : end if
370 :
371 515 : if (integratedAutoCorrTime<2._RK) cycle loopRefinement ! no need for refinement. should happen only when transitioning from compact to verbose
372 :
373 : call refineWeightedSample ( nd = RefinedChain%ndim & ! LCOV_EXCL_LINE
374 : , np = countCompact & ! LCOV_EXCL_LINE
375 : , skip = integratedAutoCorrTime & ! LCOV_EXCL_LINE
376 : , Sample = RefinedChain%LogFuncState & ! LCOV_EXCL_LINE
377 : , Weight = RefinedChain%Weight & ! LCOV_EXCL_LINE
378 : , RefinedChain = DumCFC%State & ! LCOV_EXCL_LINE
379 : , RefinedWeight = DumCFC%Weight & ! LCOV_EXCL_LINE
380 1030 : , PointCount = RefinedChain%Count(RefinedChain%numRefinement) &
381 : , refinedChainSize = refinedChainSize & ! LCOV_EXCL_LINE
382 515 : )
383 70036 : RefinedChain%Weight = DumCFC%Weight
384 260061 : RefinedChain%LogFuncState = DumCFC%State
385 :
386 : end do loopRefinement
387 :
388 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
389 :
390 354 : end subroutine getRefinedChain
391 :
392 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
393 :
394 : !> Return the refined vector of weights of the vector of weights of a weighted Markov chain.
395 : !>
396 : !> \param[in] np : The number of elements of the `Weight` vector.
397 : !> \param[in] Weight : The input vector of weights.
398 : !> \param[in] skip : The size of the jumps that have to be made through the weighted Markov chain.
399 : !> \param[in] refinedChainSize : The requested refined sample size (**optional**). If present, then the refined chain (represented by the
400 : !> : vector `Weight`) will be refined such that the resulting refined chain has the size `refinedChainSize`.
401 : !>
402 : !> \return
403 : !> `RefinedWeight` : An array of size `np`, whose elements indicate which points are present in the final refined chain.\n
404 : !> Examples:
405 : !> ```
406 : !> skip: 1
407 : !> Weight: 5, 0, 1, 3, 1
408 : !> RefinedWeight: 5, 0, 1, 3, 1
409 : !> skip: 2
410 : !> Weight: 5, 0, 1, 3, 1
411 : !> RefinedWeight: 3, 0, 0, 2, 0
412 : !> skip: 3
413 : !> Weight: 5, 0, 1, 3, 1
414 : !> RefinedWeight: 2, 0, 0, 1, 1
415 : !> ```
416 1030 : pure function getRefinedWeight(np,Weight,skip,refinedChainSize) result(RefinedWeight)
417 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
418 : !DEC$ ATTRIBUTES DLLEXPORT :: getRefinedWeight
419 : #endif
420 354 : use Math_mod, only: getCumSum
421 : use Constants_mod, only: IK, RK
422 : implicit none
423 : real(RK), intent(in) :: skip
424 : integer(IK), intent(in) :: np, Weight(np)
425 : integer(IK) , intent(in) , optional :: refinedChainSize
426 : integer(IK) :: RefinedWeight(np)
427 : integer(IK) :: ip, refinedChainSizeCounter, offset
428 154616 : real(RK) :: sumSkips, CumSumWeight(np)
429 : logical :: refinedChainSizeIsPresent
430 515 : refinedChainSizeIsPresent = present(refinedChainSize)
431 515 : if (refinedChainSizeIsPresent) refinedChainSizeCounter = 0_IK
432 154101 : CumSumWeight = real(getCumSum(np,Weight),kind=RK)
433 515 : sumSkips = skip
434 515 : offset = 1_IK
435 515 : ip = offset
436 154101 : RefinedWeight = 0_IK
437 153408 : loopOverAllSample: do
438 108744 : loopOverCurrentSample: do
439 262667 : if (sumSkips>CumSumWeight(ip)) then
440 153894 : exit loopOverCurrentSample
441 108773 : elseif (refinedChainSizeIsPresent) then
442 11605 : if (refinedChainSizeCounter==refinedChainSize) exit loopOverAllSample
443 11576 : refinedChainSizeCounter = refinedChainSizeCounter + 1_IK
444 : end if
445 108744 : RefinedWeight(ip) = RefinedWeight(ip) + 1_IK
446 108744 : sumSkips = sumSkips + skip
447 108744 : cycle loopOverCurrentSample
448 : end do loopOverCurrentSample
449 153894 : if (ip==np) then
450 515 : if (refinedChainSizeIsPresent) then
451 61 : if (refinedChainSizeCounter<refinedChainSize) then
452 29 : offset = offset + 1_IK
453 29 : if (offset==np) offset = 1_IK
454 29 : ip = offset
455 29 : sumSkips = skip
456 29 : if (offset/=1_IK) sumSkips = sumSkips + CumSumWeight(ip-1)
457 : cycle loopOverAllSample ! LCOV_EXCL_LINE
458 : end if
459 : end if
460 486 : exit loopOverAllSample
461 : end if
462 153379 : ip = ip + 1_IK
463 : end do loopOverAllSample
464 515 : end function getRefinedWeight
465 :
466 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
467 :
468 : !> Refined an input weighted sample according to the new requested weights.
469 : !>
470 : !> \param[in] nd : The number of dimensions of the input `Sample(0:nd,np)`.
471 : !> \param[in] np : The number of points in the input `Sample(0:nd,np)`.
472 : !> \param[in] skip : The jump size with which the input chain has to be refined.
473 : !> \param[in] Sample : The input 2-dimensional array of sampled states which has to be refined.
474 : !> \param[in] Weight : The weights of the sampled points.
475 : !> \param[out] RefinedChain : The refined array.
476 : !> \param[out] RefinedWeight : The vector of refined weights corresponding to the output refined array.
477 : !> \param[out] PointCount : An object of derived type [Count_type](@ref paramontechainfilecontents_mod::count_type)
478 : !> containing the number of points in the refined sample.
479 : !> \param[in] refinedChainSize : The requested refined sample size (**optional**). If the size of the refined sample is given as input,
480 : !> then the requested sample is directly generated based on the input size.
481 515 : pure subroutine refineWeightedSample(nd,np,skip,Sample,Weight,RefinedChain,RefinedWeight,PointCount,refinedChainSize)
482 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
483 : !DEC$ ATTRIBUTES DLLEXPORT :: refineWeightedSample
484 : #endif
485 515 : use Constants_mod, only: IK, RK
486 : implicit none
487 : integer(IK) , intent(in) :: nd, np, Weight(np) !, skip
488 : real(RK) , intent(in) :: Sample(0:nd,np), skip
489 : real(RK) , intent(out), allocatable :: RefinedChain(:,:)
490 : integer(IK) , intent(out), allocatable :: RefinedWeight(:)
491 : integer(IK) , intent(in) , optional :: refinedChainSize
492 : type(Count_type), intent(out) :: PointCount
493 : integer(IK) :: ip, ipRefined, npRefined, UpdatedWeight(np) ! LCOV_EXCL_LINE
494 515 : UpdatedWeight = getRefinedWeight(np,Weight,skip,refinedChainSize)
495 154101 : npRefined = count(UpdatedWeight>0)
496 515 : allocate( RefinedChain(0:nd,npRefined) , RefinedWeight(npRefined) )
497 515 : ipRefined = 0
498 515 : PointCount%verbose = 0
499 154101 : do ip = 1, np
500 154101 : if (UpdatedWeight(ip)>0) then
501 69520 : ipRefined = ipRefined + 1
502 259545 : RefinedChain(0:nd,ipRefined) = Sample(0:nd,ip)
503 69520 : RefinedWeight(ipRefined) = UpdatedWeight(ip)
504 69520 : PointCount%verbose = PointCount%verbose + RefinedWeight(ipRefined)
505 : end if
506 : end do
507 515 : PointCount%compact = npRefined
508 515 : end subroutine refineWeightedSample
509 :
510 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
511 :
512 : !> Return the best skip size through a Markov chain to refined it to the optimal requested size.
513 : !>
514 : !> \param[in] oldSampleSize : The original size of the Markov chain.
515 : !> \param[in] newSampleSize : The final desired size of the refined sample.
516 : !>
517 : !> \return
518 : !> `skip4NewSampleSize` : The computed skip size.
519 : !>
520 : !> \warning
521 : !> The condition `oldSampleSize >= newSampleSize` must always hold,
522 : !> otherwise a negative value for `skip4NewSampleSize` will be returned to indicate the occurrence of an error.
523 12 : pure function getSkip4NewSampleSize(oldSampleSize,newSampleSize) result(skip4NewSampleSize)
524 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
525 : !DEC$ ATTRIBUTES DLLEXPORT :: getSkip4NewSampleSize
526 : #endif
527 515 : use Constants_mod, only: IK, RK
528 : implicit none
529 : integer(IK) , intent(in) :: oldSampleSize,newSampleSize
530 : integer(IK) :: skip4NewSampleSize, addition, quotient
531 12 : if (oldSampleSize < newSampleSize) then
532 3 : skip4NewSampleSize = -1_IK
533 3 : return
534 : end if
535 9 : addition = 1
536 9 : quotient = oldSampleSize / newSampleSize
537 9 : if (mod(oldSampleSize,newSampleSize)==0) addition = 0
538 9 : skip4NewSampleSize = quotient + addition
539 12 : end function getSkip4NewSampleSize
540 :
541 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
542 :
543 : !> Write the computed refined chain to the specified output file.
544 : !>
545 : !> \param[in] RefinedChain : An object of class [RefinedChain_type](@ref refinedchain_type)
546 : !> containing the refined sample to be written to the output file.
547 : !> \param[in] sampleFileUnit : The unit of the file to which the sample must be written.
548 : !> \param[in] sampleFileHeaderFormat : The IO format of the header of the sample file.
549 : !> \param[in] sampleFileContentsFormat : The IO format of the contents (sampled states) in the sample file.
550 261 : subroutine writeRefinedChain(RefinedChain,sampleFileUnit,sampleFileHeaderFormat,sampleFileContentsFormat)
551 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
552 : !DEC$ ATTRIBUTES DLLEXPORT :: writeRefinedChain
553 : #endif
554 : implicit none
555 : class(RefinedChain_type), intent(in) :: RefinedChain
556 : integer(IK) , intent(in) :: sampleFileUnit
557 : character(*), intent(in) :: sampleFileHeaderFormat, sampleFileContentsFormat
558 : integer(IK) :: isample, iweight, i
559 940 : write(sampleFileUnit, sampleFileHeaderFormat) (trim(adjustl(RefinedChain%ColHeader(i)%record)),i=0,RefinedChain%ndim)
560 44794 : do isample = 1, RefinedChain%Count(RefinedChain%numRefinement)%compact
561 154870 : do iweight = 1, RefinedChain%Weight(isample)
562 154609 : write(sampleFileUnit,sampleFileContentsFormat) RefinedChain%LogFuncState(0:RefinedChain%ndim,isample)
563 : end do
564 : end do
565 261 : flush(sampleFileUnit)
566 : end subroutine writeRefinedChain ! LCOV_EXCL_LINE
567 :
568 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
569 :
570 : !> Write the computed refined chain to the specified output file.
571 : !>
572 : !> \param[in] sampleFilePath : The path to the input chain file that must be read.
573 : !> \param[in] delimiter : The delimiter used in the file.
574 : !> \param[in] ndim : The number of dimensions of the sampled states in the sample file.
575 : !> This is basically, the size of the domain of the objective function.
576 : !> \param[in] tenabled : An optional input logical value standing for `transpose-enabled`. If `.false.`,
577 : !> the input data will be naturally read according to Fortran column-wise data storage
578 : !> rule as a matrix of rank `0:nd * 1:np`. If `.false.`, the input sample file will be
579 : !> read as a matrix of rank `1:np * 0:nd`. Note that `np` represents the number of rows
580 : !> in the files (that is, the number of sampled points, whereas `nd` represents the
581 : !> number of columns in the input file (**optional**, default = `.false.`).
582 : !>
583 : !> \return
584 : !> `RefinedChain` : An object of class [RefinedChain_type](@ref refinedchain_type) containing
585 : !> the sampled states read from the specified input file.
586 176 : function readRefinedChain(sampleFilePath,delimiter,ndim,tenabled) result(RefinedChain)
587 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
588 : !DEC$ ATTRIBUTES DLLEXPORT :: readRefinedChain
589 : #endif
590 261 : use FileContents_mod, only: getNumRecordInFile
591 : use String_mod, only: String_type, num2str
592 : implicit none
593 : integer(IK) , intent(in) :: ndim
594 : character(*), intent(in) :: sampleFilePath, delimiter
595 : logical , intent(in), optional :: tenabled
596 : type(RefinedChain_type) :: RefinedChain
597 : character(*), parameter :: PROCEDURE_NAME = MODULE_NAME//"@readRefinedChain()"
598 : integer(IK) :: sampleFileUnit, isample, i
599 : logical :: tenabledDefault
600 176 : type(String_type) :: Record
601 :
602 : if (allocated(Record%value)) deallocate(Record%value) ! LCOV_EXCL_LINE
603 176 : allocate( character(99999) :: Record%value )
604 :
605 176 : RefinedChain%ndim = ndim
606 176 : RefinedChain%numRefinement = 0_IK
607 352 : allocate(RefinedChain%Count(RefinedChain%numRefinement:RefinedChain%numRefinement)) ! allocate just the zeroth level of `RefinedChain`.
608 :
609 : ! find the number of lines in the sample file
610 :
611 176 : call getNumRecordInFile(filePath=sampleFilePath, numRecord=RefinedChain%Count(RefinedChain%numRefinement)%verbose, Err=RefinedChain%Err)
612 176 : if (RefinedChain%Err%occurred) then
613 : ! LCOV_EXCL_START
614 : RefinedChain%Err%msg = PROCEDURE_NAME // RefinedChain%Err%msg
615 : return
616 : ! LCOV_EXCL_STOP
617 : end if
618 :
619 176 : RefinedChain%Count(RefinedChain%numRefinement)%verbose = RefinedChain%Count(RefinedChain%numRefinement)%verbose - 1_IK ! remove header from the count
620 :
621 : open( newunit = sampleFileUnit &
622 : , file = sampleFilePath &
623 : , status = "old" &
624 : #if defined INTEL_COMPILER_ENABLED && defined OS_IS_WINDOWS
625 : , SHARED &
626 : #endif
627 176 : )
628 :
629 : ! read header
630 :
631 176 : read(sampleFileUnit,"(A)") Record%value
632 1882 : Record%Parts = Record%split(Record%value,delimiter,Record%nPart)
633 176 : if (Record%nPart/=ndim+1_IK) then
634 : ! LCOV_EXCL_START
635 : RefinedChain%Err%occurred = .true.
636 : RefinedChain%Err%msg = PROCEDURE_NAME // ": The number of header columns ("//num2str(Record%nPart)//") is not equal to ndim + 1: "//num2str(ndim+1_IK)
637 : return
638 : ! LCOV_EXCL_STOP
639 : end if
640 :
641 686 : allocate(RefinedChain%ColHeader(0:ndim))
642 686 : do i = 0, ndim
643 686 : RefinedChain%ColHeader(i)%record = Record%Parts(i+1)%record
644 : end do
645 :
646 : ! read contents
647 :
648 176 : tenabledDefault = .false.
649 176 : if (present(tenabled)) tenabledDefault = tenabled
650 :
651 176 : if (tenabledDefault) then
652 :
653 168 : allocate( RefinedChain%LogFuncState(RefinedChain%Count(RefinedChain%numRefinement)%verbose, 0:ndim) )
654 248454 : do isample = 1, RefinedChain%Count(RefinedChain%numRefinement)%verbose
655 248286 : read(sampleFileUnit, "(A)") Record%value
656 1981090 : Record%Parts = Record%split(trim(adjustl(Record%value)),delimiter)
657 990714 : do i = 0, ndim
658 990546 : read(Record%Parts(i+1)%record,*) RefinedChain%LogFuncState(isample,i)
659 : end do
660 : end do
661 :
662 : else
663 :
664 8 : allocate( RefinedChain%LogFuncState(0:ndim, RefinedChain%Count(RefinedChain%numRefinement)%verbose) )
665 88 : do isample = 1, RefinedChain%Count(RefinedChain%numRefinement)%verbose
666 80 : read(sampleFileUnit, "(A)") Record%value
667 640 : Record%Parts = Record%split(trim(adjustl(Record%value)),delimiter)
668 328 : do i = 0, ndim
669 320 : read(Record%Parts(i+1)%record,*) RefinedChain%LogFuncState(i,isample)
670 : end do
671 : end do
672 :
673 : end if
674 :
675 176 : close(sampleFileUnit)
676 :
677 862 : end function readRefinedChain
678 :
679 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
680 :
681 : end module ParaMCMCRefinedChain_mod ! LCOV_EXCL_LINE
|
