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 : !> \brief
18 : !> This file contains the implementation details of the 2D routines under the generic interface [pm_sampleACT](@ref pm_sampleACT).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Saturday 4:40 PM, August 21, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if getACT_ENABLED && D1_ENABLED && DEF_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 1 : act = getACT(seq, batchMeans)
32 :
33 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34 : #elif getACT_ENABLED && D1_ENABLED && BMM_ENABLED
35 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 :
37 : integer(IK) :: smin, smax, step, isize, lenseq
38 2 : smin = method%smin
39 2 : step = method%step
40 2 : smax = method%smax
41 2 : lenseq = size(seq, 1, IK)
42 2 : CHECK_ASSERTION(__LINE__, 0 < lenseq, SK_"@getACT(): The condition `0 < size(seq)` must hold. size(seq) = "//getStr(lenseq))
43 2 : CHECK_ASSERTION(__LINE__, 1_IK < smin, SK_"@getACT(): The condition `1 < method%smin` must hold. method%smin = "//getStr(smin))
44 6 : CHECK_ASSERTION(__LINE__, smax == huge(0_IK) .or. (smax <= smin .and. smax <= lenseq / 2_IK), SK_"@getACT(): The condition `method%smin <= method%smax <= size(seq, 1) / 2` must hold. method%smax, size(seq) = "//getStr([smax, lenseq]))
45 2 : CHECK_ASSERTION(__LINE__, 0_IK < step, SK_"@getACT(): The condition `0 < method%step` must hold. method%step = "//getStr(step))
46 2 : if (1_IK < lenseq) then
47 2 : smax = min(smax, lenseq / 2_IK)
48 0 : act = -huge(act)
49 2 : do isize = smin, smax, step
50 : #if ONE_ENABLED
51 5085 : act = max(act, getACT(seq, batchMeans_type(isize)))
52 : #elif WTI_ENABLED
53 0 : act = max(act, getACT(seq, weight, batchMeans_type(isize)))
54 : #else
55 : #error "Unrecognized interface."
56 : #endif
57 : end do
58 : else
59 0 : act = lenseq
60 : end if
61 :
62 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 : #elif getACT_ENABLED && D1_ENABLED && BMD_ENABLED
64 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65 :
66 : integer(IK) :: batchSize, lenseq, weisum
67 : real(TKC), allocatable :: batchMean(:)
68 : real(TKC) :: batchSizeInv, avgSeq, sumDiffSq, varBatchMean, avgBatchMean
69 : integer(IK) :: nbatch, ibatch, seqLenEffective
70 : #if WTI_ENABLED
71 : real(TKC) :: diffSq
72 104 : integer(IK) :: cumSumWei(size(seq, 1, IK)), currentBatchEndLoc, iseq, iseqVerbose
73 52 : lenseq = size(seq, 1, IK)
74 52 : call setCumSum(cumSumWei, weight)
75 156 : CHECK_ASSERTION(__LINE__, size(weight, 1, IK) == lenseq, SK_"@getACT(): The condition `size(weight) == size(seq)` must hold. size(weight), size(seq) = "//getStr([size(weight, 1, IK), lenseq]))
76 52 : weisum = cumSumWei(lenseq)
77 : #elif ONE_ENABLED
78 : integer(IK) :: iseqBeg, iseqEnd
79 9822 : lenseq = size(seq, 1, IK)
80 : weisum = lenseq
81 : #else
82 : #error "Unrecognized interface."
83 : #endif
84 9874 : CHECK_ASSERTION(__LINE__, 0 < lenseq, SK_"@getACT(): The condition `0 < size(seq)` must hold. size(seq) = "//getStr(lenseq))
85 9874 : if (lenseq < 2_IK) then
86 0 : act = lenseq
87 0 : return
88 : end if
89 :
90 : ! Compute batch size and count.
91 :
92 0 : act = 1._TKC
93 9874 : batchSize = method%size
94 9874 : if (batchSize == 0_IK) batchSize = int(real(weisum)**real(2./3.), IK)
95 9874 : batchSizeInv = 1._TKC / real(batchSize, TKC)
96 9874 : nbatch = weisum / batchSize
97 9874 : CHECK_ASSERTION(__LINE__, 0_IK < batchSize, SK_"@getACT(): The condition `0 < method%size` must hold. method%size = "//getStr(batchSize))
98 9874 : if (nbatch < 2_IK) then
99 : !check_assertion(__LINE__, method%size == 0_IK, SK_"@getACT(): The condition `method%size * 2 < sum(weight)` must hold (there must be more than one batches to compute ACT). method%size, sum(weight) = "//getStr([batchSize, weisum]))
100 0 : act = weisum
101 0 : return
102 : end if
103 9874 : seqLenEffective = nbatch * batchSize
104 : ! xxx: here goes another GFortran 7.3 bug: `batchMean` is assumed already allocated, despite the first appearance here.
105 9874 : if (allocated(batchMean)) deallocate(batchMean)
106 9874 : allocate(batchMean(nbatch))
107 :
108 : ! \todo
109 : ! iterate from the end to the beginning of the chain to ignore initial points instead of the last points.
110 : ! This would be beneficial for MCMC samples.
111 :
112 : ! Compute the Batch means vector and mean of the whole (weighted) sequence.
113 :
114 0 : sumDiffSq = 0._TKC
115 0 : avgSeq = 0._TKC
116 : #if WTI_ENABLED
117 : iseq = 1
118 : ibatch = 1
119 : iseqVerbose = 0
120 : currentBatchEndLoc = batchSize
121 52 : batchMean(ibatch) = 0._TKC
122 381479 : loopOverWeight: do
123 381531 : iseqVerbose = iseqVerbose + 1
124 381531 : if (cumSumWei(iseq) < iseqVerbose) iseq = iseq + 1
125 381531 : if (currentBatchEndLoc < iseqVerbose) then ! we are done with the current batch
126 712 : avgSeq = avgSeq + batchMean(ibatch)
127 712 : if (seqLenEffective < iseqVerbose) exit loopOverWeight ! condition equivalent to currentBatchEndLoc == seqLenEffective.
128 660 : currentBatchEndLoc = currentBatchEndLoc + batchSize
129 660 : ibatch = ibatch + 1
130 660 : batchMean(ibatch) = 0._TKC
131 : end if
132 381531 : batchMean(ibatch) = batchMean(ibatch) + seq(iseq)
133 : end do loopOverWeight
134 764 : batchMean = batchMean * batchSizeInv
135 52 : avgSeq = avgSeq / real(seqLenEffective, TKC) ! whole sequence mean.
136 : ! compute whole sequence variance.
137 : iseq = 1
138 : iseqVerbose = 0
139 52 : diffSq = (seq(iseq) - avgSeq)**2
140 381479 : loopComputeVariance: do
141 381531 : iseqVerbose = iseqVerbose + 1
142 381531 : if (seqLenEffective < iseqVerbose) exit loopComputeVariance
143 381479 : if (cumSumWei(iseq) < iseqVerbose) then
144 365418 : iseq = iseq + 1 ! by definition, iseq never become > lenseq, otherwise it leads to disastrous errors
145 365418 : diffSq = (seq(iseq) - avgSeq)**2
146 : end if
147 381531 : sumDiffSq = sumDiffSq + diffSq ! whole sequence variance.
148 : end do loopComputeVariance
149 : #elif ONE_ENABLED
150 : iseqBeg = 1
151 : iseqEnd = 0
152 166727 : do ibatch = 1, nbatch
153 156905 : iseqEnd = iseqEnd + batchSize
154 82887268 : batchMean(ibatch) = sum(seq(iseqBeg : iseqEnd))
155 156905 : avgSeq = avgSeq + batchMean(ibatch)
156 156905 : batchMean(ibatch) = batchMean(ibatch) * batchSizeInv
157 166727 : iseqBeg = iseqEnd + 1_IK
158 : end do
159 9822 : avgSeq = avgSeq / real(seqLenEffective, TKC) ! whole sequence mean.
160 82740185 : sumDiffSq = sum((seq(1 : seqLenEffective) - avgSeq)**2) ! whole sequence variance.
161 :
162 : #else
163 : #error "Unrecognized interface."
164 : #endif
165 167491 : avgBatchMean = sum(batchMean) / real(nbatch, TKC) ! batch means
166 167491 : varBatchMean = sum((batchMean - avgBatchMean)**2) / real(nbatch - 1, TKC) ! batch variances
167 : !act = varBatchMean * seqLenEffective * (seqLenEffective - 1) / sumDiffSq
168 : !act = batchSize * varBatchMean * (seqLenEffective - 1) / sumDiffSq
169 9874 : act = varBatchMean * seqLenEffective * (seqLenEffective - 1) / sumDiffSq / real(nbatch, TKC)
170 : !print *, avgSeq, sumDiffSq
171 : !print *, sum(batchMean) / real(nbatch, TKC), varBatchMean
172 : !print *, nbatch, batchSize
173 : !print *, act
174 :
175 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
176 : #elif getACT_ENABLED && D1_ENABLED && (CSD_ENABLED || CSM_ENABLED)
177 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
178 :
179 8 : real(TKC) :: mean, sample(size(seq, 1, IK))
180 : integer(IK) :: lenseq, weisum
181 4 : lenseq = size(seq, 1, IK)
182 4 : CHECK_ASSERTION(__LINE__, 1 < lenseq, SK_"@getACT(): The condition `0 < size(seq)` must hold. size(seq) = "//getStr(lenseq))
183 4 : if (lenseq < 2_IK) then
184 0 : act = lenseq
185 0 : return
186 : end if
187 : #if WTI_ENABLED
188 0 : CHECK_ASSERTION(__LINE__, size(weight, 1, IK) == size(seq, 1, IK), SK_"@getACT(): The condition `size(weight) == size(seq)` must hold. size(weight), size(seq) = "//getStr([size(weight, 1, IK), size(seq, 1, IK)]))
189 0 : weisum = sum(weight)
190 0 : sample = getVerbose(seq, weight, weisum)
191 0 : mean = sum(sample * weight) / real(weisum, TKC)
192 0 : sample = sample - mean
193 : #elif ONE_ENABLED
194 : weisum = size(seq, 1, IK)
195 36992 : mean = sum(seq) / real(weisum, TKC)
196 36992 : sample = seq - mean
197 : #else
198 : #error "Unrecognized interface."
199 : #endif
200 : #if CSD_ENABLED
201 : block
202 : integer(IK) :: i
203 : real(TKC) :: signif
204 3 : signif = real(method%signif, TKC)
205 3 : CHECK_ASSERTION(__LINE__, 0 <= method%signif, SK_"@getACT(): The condition `0 <= method%signif` must hold. method%signif = "//getStr(signif))
206 27744 : sample = getACF(sample, norm = stdscale)
207 : ! For autocorrelation, under the assumption of a completely random series, the ACF standard error reduces to `sqrt(1 / ndata)`.
208 0 : act = 0._TKC
209 3 : signif = signif * sqrt(1._TKC / weisum)
210 275 : do i = 1, size(sample, 1, IK)
211 275 : if (sample(i) < signif) exit
212 275 : act = act + sample(i)
213 : end do
214 3 : act = 2 * act - 1
215 : end block
216 : #elif CSM_ENABLED
217 9248 : sample = getACF(sample, norm = stdscale)
218 9249 : act = 2 * maxval(getCumSum(sample)) - 1
219 : #else
220 : #error "Unrecognized interface."
221 : #endif
222 :
223 : #else
224 : !%%%%%%%%%%%%%%%%%%%%%%%%
225 : #error "Unrecognized interface."
226 : !%%%%%%%%%%%%%%%%%%%%%%%%
227 : #endif
228 : #undef TYPE_OF_SEQ
229 : #undef GET_ABSQ
230 : #undef GET_RE
|
