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 implementations of procedures [pm_clustering](@ref pm_clustering).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #define INTRINSIC_ENABLED 0
28 :
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 : #if setMember_ENABLED && D0_D1_ENABLED
31 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32 :
33 : real(RKC) :: temp
34 : integer(IK) :: icls
35 1 : CHECK_ASSERTION(__LINE__, 0_IK < size(center, 1, IK), SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
36 : ! First cluster.
37 : disq = 0._RKC
38 1 : disq = (sample - center(1))**2
39 1 : membership = 1
40 : ! All the rest of clusters.
41 3 : do icls = 2, size(center, 1, IK)
42 2 : temp = (sample - center(icls))**2
43 3 : if (temp < disq) then
44 0 : disq = temp
45 0 : membership = icls
46 : end if
47 : end do
48 :
49 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50 : #elif setMember_ENABLED && D1_D1_ENABLED
51 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 :
53 : real(RKC) :: temp
54 : integer(IK) :: icls, isam, nsam
55 1 : nsam = size(sample, 1, IK)
56 1 : CHECK_ASSERTION(__LINE__, 0_IK < size(center, 1, IK), SK_"@setMember(): The condition `0 < size(center, 2)` must hold. size(center, 1) = "//getStr(size(center, 1, IK)))
57 4 : CHECK_ASSERTION(__LINE__, nsam == size(disq, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), size(disq, 1, IK)]))
58 4 : CHECK_ASSERTION(__LINE__, nsam == size(membership, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(membership)` must hold. shape(sample), size(membership) = "//getStr([shape(sample, IK), size(membership, 1, IK)]))
59 : ! First cluster.
60 5 : do isam = 1, nsam
61 4 : disq(isam) = (sample(isam) - center(1))**2
62 5 : membership(isam) = 1
63 : end do
64 : ! All the rest of clusters.
65 3 : do icls = 2, size(center, 1, IK)
66 11 : do isam = 1, nsam
67 8 : temp = (sample(isam) - center(icls))**2
68 10 : if (temp < disq(isam)) then
69 2 : disq(isam) = temp
70 2 : membership(isam) = icls
71 : end if
72 : end do
73 : end do
74 :
75 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
76 : #elif setMember_ENABLED && D1_D2_ENABLED
77 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78 :
79 : real(RKC) :: temp
80 : integer(IK) :: icls, idim, ndim
81 1 : ndim = size(sample, 1, IK)
82 1 : CHECK_ASSERTION(__LINE__, 0_IK < size(center, rank(center), IK), SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
83 3 : CHECK_ASSERTION(__LINE__, ndim == size(center, 1, IK), SK_"@setMember(): The condition `size(sample, 1) == size(center, 1)` must hold. size(sample, 1), size(center, 1) = "//getStr([size(sample, 1, IK), size(center, 1, IK)]))
84 : ! First cluster.
85 : !disq = sum((sample - center(1 : ndim, 1))**2)
86 1 : disq = 0._RKC
87 4 : do idim = 1, ndim
88 4 : disq = disq + (sample(idim) - center(idim, 1))**2
89 : end do
90 1 : membership = 1
91 : ! All the rest of clusters.
92 3 : do icls = 2, size(center, 2, IK)
93 : !temp = sum((sample - center(1 : ndim, icls))**2)
94 0 : temp = 0._RKC
95 8 : do idim = 1, ndim
96 8 : temp = temp + (sample(idim) - center(idim, icls))**2
97 : end do
98 3 : if (temp < disq) then
99 1 : disq = temp
100 1 : membership = icls
101 : end if
102 : end do
103 :
104 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105 : #elif setMember_ENABLED && D2_D2_ENABLED
106 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107 :
108 : real(RKC) :: temp
109 : integer(IK) :: icls, idim, ndim, isam, nsam
110 38 : ndim = size(sample, 1, IK)
111 38 : nsam = size(sample, 2, IK)
112 38 : CHECK_ASSERTION(__LINE__, 0_IK < size(center, 2, IK), SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
113 228 : CHECK_ASSERTION(__LINE__, nsam == size(disq, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), size(disq, 1, IK)]))
114 114 : CHECK_ASSERTION(__LINE__, ndim == size(center, 1, IK), SK_"@setMember(): The condition `size(sample, 1) == size(center, 1)` must hold. size(sample, 1), size(center, 1) = "//getStr([size(sample, 1, IK), size(center, 1, IK)]))
115 228 : CHECK_ASSERTION(__LINE__, nsam == size(membership, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(membership, 1)` must hold. shape(sample), size(membership, 1) = "//getStr([shape(sample, IK), size(membership, 1, IK)]))
116 : ! First cluster.
117 : icls = 1
118 130086 : do isam = 1, nsam
119 : !disq(isam) = sum((center(1 : ndim, icls) - sample(1 : ndim, isam))**2)
120 130048 : disq(isam) = 0._RKC
121 390180 : do idim = 1, ndim
122 390180 : disq(isam) = disq(isam) + (sample(idim, isam) - center(idim, icls))**2
123 : end do
124 130086 : membership(isam) = icls
125 : end do
126 : ! All the rest of clusters.
127 165 : do icls = 2, size(center, 2, IK)
128 520283 : do isam = 1, nsam
129 : !temp = sum((center(1 : ndim, icls) - sample(1 : ndim, isam))**2)
130 0 : temp = 0._RKC
131 1560452 : do idim = 1, ndim
132 1560452 : temp = temp + (sample(idim, isam) - center(idim, icls))**2
133 : end do
134 520245 : if (temp < disq(isam)) then
135 184792 : disq(isam) = temp
136 184792 : membership(isam) = icls
137 : end if
138 : end do
139 : end do
140 :
141 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142 : #elif setCenter_ENABLED && D1_D1_ENABLED
143 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
144 :
145 : integer(IK) :: icls, ncls, isam, nsam
146 :
147 1 : nsam = ubound(sample, rank(sample), IK)
148 1 : ncls = ubound(center, rank(center), IK)
149 :
150 5 : CHECK_ASSERTION(__LINE__, all(0._RKC <= disq), SK_"@setMember(): The condition `all(0 <= disq)` must hold. disq = "//getStr(disq))
151 5 : CHECK_ASSERTION(__LINE__, nsam == ubound(disq, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
152 : !CHECK_ASSERTION(__LINE__, nsam >= ubound(center, rank(center), IK), SK_"@setMember(): The condition `size(sample, rank(sample)) >= size(center, rank(center))` must hold. shape(sample), shape(center) = "//getStr([shape(sample, IK), shape(center, IK)]))
153 5 : CHECK_ASSERTION(__LINE__, nsam == ubound(membership, 1, IK), SK_"@setMember(): The condition `size(sample, rank(sample)) == size(membership)` must hold. shape(sample), size(membership) = "//getStr([shape(sample, IK), ubound(membership, 1, IK)]))
154 5 : CHECK_ASSERTION(__LINE__, ncls == ubound(potential, 1, IK), SK_"@setMember(): The condition `size(center, rank(center)) == size(potential)` must hold. shape(center), size(potential) = "//getStr([shape(center, IK), ubound(potential, 1, IK)]))
155 5 : CHECK_ASSERTION(__LINE__, ncls == ubound(size, 1, IK), SK_"@setMember(): The condition `size(center, rank(center)) == size(size)` must hold. shape(center), size(size) = "//getStr([shape(center, IK), ubound(size, 1, IK)]))
156 1 : CHECK_ASSERTION(__LINE__, 0_IK < ncls, SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
157 15 : CHECK_ASSERTION(__LINE__, all(0_IK < membership .and. membership <= ncls), SK_"@setMember(): The condition `all(0 < membership .and. membership <= size(center, rank(center)))` must hold. size(center, rank(center)), membership = "//getStr([ncls, membership]))
158 :
159 : ! Initialize.
160 :
161 : do concurrent(icls = 1 : ncls)
162 2 : potential(icls) = 0._RKC
163 2 : center(icls) = 0._RKC
164 3 : size(icls) = 0_IK
165 : end do
166 :
167 : ! compute new centers.
168 :
169 5 : do isam = 1, nsam
170 4 : center(membership(isam)) = center(membership(isam)) + sample(isam)
171 4 : potential(membership(isam)) = potential(membership(isam)) + disq(isam)
172 5 : size(membership(isam)) = size(membership(isam)) + 1
173 : end do
174 :
175 : ! normalize new centers.
176 :
177 : do concurrent(icls = 1 : ncls)
178 3 : if (0_IK < size(icls)) center(icls) = center(icls) / real(size(icls), RKC)
179 : end do
180 :
181 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182 : #elif setCenter_ENABLED && D2_D2_ENABLED
183 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184 :
185 : real(RKC) :: sizeinv
186 : integer(IK) :: icls, ncls, ndim, isam, nsam
187 :
188 36 : ndim = ubound(sample, 1, IK)
189 36 : nsam = ubound(sample, 2, IK)
190 36 : ncls = ubound(center, 2, IK)
191 :
192 125080 : CHECK_ASSERTION(__LINE__, all(0._RKC <= disq), SK_"@setMember(): The condition `all(0 <= disq)` must hold. disq = "//getStr(disq))
193 252 : CHECK_ASSERTION(__LINE__, nsam == ubound(disq, 1, IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) == ubound(disq, 1)` must hold. shape(sample), ubound(disq, 1) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
194 144 : CHECK_ASSERTION(__LINE__, ndim == ubound(center, 1, IK), SK_"@setMember(): The condition `ubound(sample, 1) == ubound(center, 1)` must hold. ubound(sample, 1), ubound(center, 1) = "//getStr([ubound(sample, 1, IK), ubound(center, 1, IK)]))
195 : !CHECK_ASSERTION(__LINE__, nsam >= ubound(center, rank(center), IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) >= ubound(center, rank(center))` must hold. shape(sample), shape(center) = "//getStr([shape(sample, IK), shape(center, IK)]))
196 252 : CHECK_ASSERTION(__LINE__, nsam == ubound(membership, 1, IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) == ubound(membership, 1)` must hold. shape(sample), ubound(membership, 1) = "//getStr([shape(sample, IK), ubound(membership, 1, IK)]))
197 252 : CHECK_ASSERTION(__LINE__, ncls == ubound(potential, 1, IK), SK_"@setMember(): The condition `ubound(center, rank(center)) == ubound(potential, 1)` must hold. shape(center), ubound(potential, 1) = "//getStr([shape(center, IK), ubound(potential, 1, IK)]))
198 252 : CHECK_ASSERTION(__LINE__, ncls == ubound(size, 1, IK), SK_"@setMember(): The condition `ubound(center, rank(center)) == ubound(size, 1)` must hold. shape(center), ubound(size, 1) = "//getStr([shape(center, IK), ubound(size, 1, IK)]))
199 36 : CHECK_ASSERTION(__LINE__, 0_IK < ncls, SK_"@setMember(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
200 375240 : CHECK_ASSERTION(__LINE__, all(0 < membership .and. membership <= ncls), SK_"@setMember(): The condition `all(0 < membership .and. membership <= size(center, rank(center)))` must hold. size(center, rank(center)), membership = "//getStr([ubound(center, 2, IK), membership]))
201 :
202 : ! Initialize.
203 :
204 : do concurrent(icls = 1 : ncls)
205 489 : center(1 : ndim, icls) = 0._RKC
206 157 : potential(icls) = 0._RKC
207 193 : size(icls) = 0_IK
208 : end do
209 :
210 : ! compute new centers.
211 :
212 125080 : do isam = 1, nsam
213 375164 : center(1 : ndim, membership(isam)) = center(1 : ndim, membership(isam)) + sample(1 : ndim, isam)
214 125044 : potential(membership(isam)) = potential(membership(isam)) + disq(isam)
215 125080 : size(membership(isam)) = size(membership(isam)) + 1
216 : end do
217 :
218 : ! normalize new centers.
219 :
220 : do concurrent(icls = 1 : ncls)
221 193 : if (0_IK < size(icls)) then
222 157 : sizeinv = 1._RKC / real(size(icls), RKC)
223 489 : center(1 : ndim, icls) = center(1 : ndim, icls) * sizeinv
224 : end if
225 : end do
226 :
227 : !%%%%%%%%%%%%%%%%%%
228 : #elif setKmeansPP_ENABLED
229 : !%%%%%%%%%%%%%%%%%%
230 :
231 : real(RKC) :: temp
232 39 : real(RKC) :: cumsumdisq(0 : ubound(sample, 2, IK)) ! Sum of distance squared of points up to the given index of the vector.
233 : integer(IK) :: idim, icls, isam, ndim, nsam, cid
234 : #if Optional_ENABLED
235 : integer(IK) :: ncls
236 12 : ncls = ubound(center, 2, IK)
237 60 : CHECK_ASSERTION(__LINE__, ubound(sample, 1, IK) == ubound(center, 1, IK), SK_"@setKmeansPP(): The condition `ubound(sample, 1) == ubound(center, 1)` must hold. ubound(sample, 1), ubound(center, 1) = "//getStr([ubound(sample, 1, IK), ubound(center, 1, IK)]))
238 48 : CHECK_ASSERTION(__LINE__, ncls == ubound(potential, 1, IK), SK_"@setKmeansPP(): The condition `ncls == ubound(potential, 1)` must hold. ncls, ubound(potential, 1) = "//getStr([ncls, ubound(potential, 1, IK)]))
239 48 : CHECK_ASSERTION(__LINE__, ncls == ubound(size, 1, IK), SK_"@setKmeansPP(): The condition `ncls == ubound(size, 1)` must hold. ncls, ubound(size, 1) = "//getStr([ncls, ubound(size, 1, IK)]))
240 : #elif !Default_ENABLED
241 : #error "Unrecognized interface."
242 : #endif
243 10 : ndim = ubound(sample, 1, IK)
244 22 : nsam = ubound(sample, 2, IK)
245 22 : CHECK_ASSERTION(__LINE__, 0_IK < ncls, SK_"@setKmeansPP(): The condition `0 < ncls` must hold. ncls = "//getStr(ncls))
246 22 : CHECK_ASSERTION(__LINE__, 0_IK < ubound(sample, 2, IK), SK_"@setKmeansPP(): The condition `0 < ubound(sample, rank(sample))` must hold. shape(sample) = "//getStr(shape(sample, IK)))
247 154 : CHECK_ASSERTION(__LINE__, nsam == ubound(disq, 1, IK), SK_"@setMember(): The condition `ubound(sample, rank(sample)) == ubound(disq, 1)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
248 88 : CHECK_ASSERTION(__LINE__, nsam == ubound(membership, 1, IK), SK_"@setKmeansPP(): The condition `ubound(sample, rank(sample)) == ubound(membership, 1)` must hold. ubound(sample, rank(sample)), ubound(membership, 1) = "//getStr([nsam, ubound(membership, 1, IK)]))
249 66 : CHECK_ASSERTION(__LINE__, ncls <= nsam, SK_"@setKmeansPP(): The condition `ncls < ubound(sample, rank(sample))` must hold. ncls, ubound(sample, rank(sample)) = "//getStr([ncls, nsam]))
250 :
251 22 : cumsumdisq(0) = 0._RKC ! This element must always remain zero.
252 :
253 : ! Define the first cluster center.
254 :
255 : icls = 1
256 22 : call setUnifRand(rng, cid, 1_IK, nsam) ! This random assignment is unnecessary if the `sample` is in random order.
257 : #if Optional_ENABLED
258 52 : center(1 : ndim, icls) = sample(1 : ndim, cid)
259 : #endif
260 10095 : do isam = 1, nsam
261 : ! Find the distance-squared of each sample to the #icls cluster cid.
262 : #if INTRINSIC_ENABLED
263 : disq(isam) = sum((sample(1 : ndim, isam) - sample(1 : ndim, cid))**2)
264 : #else
265 30298 : disq(isam) = 0._RKC; do idim = 1, ndim; disq(isam) = disq(isam) + (sample(idim, isam) - sample(idim, cid))**2; end do
266 : #endif
267 10073 : cumsumdisq(isam) = cumsumdisq(isam - 1) + disq(isam)
268 10095 : membership(isam) = icls
269 : end do
270 :
271 : ! Find the second cluster center.
272 :
273 22 : call setUnifRand(rng, temp)
274 22 : temp = temp * cumsumdisq(nsam)
275 : #if INTRINSIC_ENABLED
276 : cid = minloc(cumsumdisq, 1, temp < cumsumdisq) - 1 ! `-1` compensates for the 0 starting index of `cumsumdisq`.
277 : !if (cid < 1) cid = nsam + 1
278 : #else
279 7062 : do cid = 1, nsam; if (temp < cumsumdisq(cid)) exit; end do
280 : #endif
281 : !print *, "ndim, nsam, ncls, cid, temp", ndim, nsam, ncls, cid, temp, cumsumdisq(1:nsam)
282 :
283 : ! Take care of the unusual case of only one cluster with a single sample.
284 :
285 22 : if (1_IK == ncls) then
286 : #if Optional_ENABLED
287 13 : potential(1) = sum(disq)
288 5 : size(1) = nsam
289 : #endif
290 : return
291 : end if
292 :
293 : ! Find the rest of cluster centers.
294 :
295 40 : do icls = 2, ncls - 1
296 : #if Optional_ENABLED
297 53 : center(1 : ndim, icls) = sample(1 : ndim, cid)
298 : #endif
299 : !if (nsam < cid) exit ! only if nsam < ncls.
300 30144 : do isam = 1, nsam
301 : #if INTRINSIC_ENABLED
302 : temp = sum((sample(1 : ndim, isam) - sample(1 : ndim, cid))**2)
303 : #else
304 90519 : temp = 0._RKC; do idim = 1, ndim; temp = temp + (sample(idim, isam) - sample(idim, cid))**2; end do
305 : #endif
306 30119 : if (temp < disq(isam)) then
307 10608 : disq(isam) = temp
308 10608 : membership(isam) = icls
309 : end if
310 30144 : cumsumdisq(isam) = cumsumdisq(isam - 1) + disq(isam)
311 : end do
312 25 : call setUnifRand(rng, temp)
313 25 : temp = temp * cumsumdisq(nsam)
314 : #if INTRINSIC_ENABLED
315 : cid = minloc(cumsumdisq, 1, temp < cumsumdisq) - 1
316 : !if (cid < 1) cid = nsam + 1
317 : #else
318 14958 : do cid = 1, nsam; if (temp < cumsumdisq(cid)) exit; end do
319 : #endif
320 : end do
321 :
322 : ! Initialize the output.
323 :
324 : #if Optional_ENABLED
325 : do concurrent(icls = 1 : ncls)
326 : !center(1 : ndim, icls) = 0._RKC
327 27 : potential(icls) = 0._RKC
328 34 : size(icls) = 0_IK
329 : end do
330 29 : center(1 : ndim, ncls) = sample(1 : ndim, cid)
331 : #endif
332 :
333 : ! Final re-computation of the cluster centers and sizes and sample memberships.
334 :
335 10077 : do isam = 1, nsam
336 : #if INTRINSIC_ENABLED
337 : temp = sum((sample(1 : ndim, isam) - sample(1 : ndim, cid))**2)
338 : #else
339 30253 : temp = 0._RKC; do idim = 1, ndim; temp = temp + (sample(idim, isam) - sample(idim, cid))**2; end do
340 : #endif
341 10070 : if (temp < disq(isam)) then
342 1475 : disq(isam) = temp
343 1475 : membership(isam) = ncls
344 : end if
345 : #if Optional_ENABLED
346 : !center(1 : ndim, membership(isam)) = center(1 : ndim, membership(isam)) + sample(1 : ndim, isam)
347 10025 : potential(membership(isam)) = potential(membership(isam)) + disq(isam)
348 10032 : size(membership(isam)) = size(membership(isam)) + 1
349 : #endif
350 : end do
351 :
352 : ! Normalize sample. This requires the sample size to be larger than the number of clusters to have `all(size > 0)` hold.
353 :
354 : !do concurrent(icls = 1 : ncls)
355 : ! temp = 1._RKC / real(size(icls), RKC)
356 : ! center(1 : ndim, icls) = center(1 : ndim, icls) * temp
357 : !end do
358 :
359 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
360 : #elif setKmeans_ENABLED && Init_ENABLED
361 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
362 :
363 11 : integer(IK) :: retin, sizemin, retinxam, membersnew(ubound(membership, 1, IK))
364 :
365 5051 : CHECK_ASSERTION(__LINE__, all(0._RKC <= disq), SK_"@setKmeans(): The condition `all(0 <= disq)` must hold. disq = "//getStr(disq))
366 88 : CHECK_ASSERTION(__LINE__, ubound(sample, 2, IK) == ubound(disq, 1, IK), SK_"@setKmeans(): The condition `size(sample, rank(sample)) == size(disq)` must hold. shape(sample), size(disq) = "//getStr([shape(sample, IK), ubound(disq, 1, IK)]))
367 55 : CHECK_ASSERTION(__LINE__, ubound(sample, 1, IK) == ubound(center, 1, IK), SK_"@setKmeans(): The condition `size(sample, 1) == size(center, 1)` must hold. size(sample, 1), size(center, 1) = "//getStr([ubound(sample, 1, IK), ubound(center, 1, IK)]))
368 : !CHECK_ASSERTION(__LINE__, ubound(sample, 2, IK) >= ubound(center, rank(center), IK), SK_"@setKmeans(): The condition `size(sample, rank(sample)) >= size(center, rank(center))` must hold. shape(sample), shape(center) = "//getStr([shape(sample, IK), shape(center, IK)]))
369 77 : CHECK_ASSERTION(__LINE__, ubound(sample, 2, IK) == ubound(membership, 1, IK), SK_"@setKmeans(): The condition `size(sample, rank(sample)) == size(membership)` must hold. shape(sample), size(membership) = "//getStr([shape(sample, IK), ubound(membership, 1, IK)]))
370 88 : CHECK_ASSERTION(__LINE__, ubound(center, 2, IK) == ubound(potential, 1, IK), SK_"@setKmeans(): The condition `size(center, rank(center)) == size(potential)` must hold. shape(center), size(potential) = "//getStr([shape(center, IK), ubound(potential, 1, IK)]))
371 77 : CHECK_ASSERTION(__LINE__, ubound(center, 2, IK) == ubound(size, 1, IK), SK_"@setKmeans(): The condition `size(center, rank(center)) == size(size)` must hold. shape(center), size(size) = "//getStr([shape(center, IK), ubound(size, 1, IK)]))
372 11 : CHECK_ASSERTION(__LINE__, ubound(center, 2, IK) > 0_IK, SK_"@setKmeans(): The condition `0 < size(center, rank(center))` must hold. shape(center) = "//getStr(shape(center, IK)))
373 15153 : CHECK_ASSERTION(__LINE__, all(0 < membership .and. membership <= ubound(center, 2, IK)), SK_"@setKmeans(): The condition `all(0 < membership .and. membership <= size(center, rank(center)))` must hold. size(center, rank(center)), membership = "//getStr([ubound(center, 2, IK), membership]))
374 :
375 11 : if (present(maxniter)) then
376 0 : CHECK_ASSERTION(__LINE__, 0_IK <= maxniter, SK_"@setKmeans(): The condition `0 <= maxniter` must hold. maxniter = "//getStr(maxniter))
377 : retinxam = maxniter
378 : else
379 : retinxam = 300_IK
380 : end if
381 11 : if (present(minsize)) then
382 0 : CHECK_ASSERTION(__LINE__, 0_IK <= minsize, SK_"@setKmeans(): The condition `0 <= minsize` must hold. minsize = "//getStr(minsize))
383 : sizemin = minsize
384 : else
385 : sizemin = 1_IK
386 : end if
387 :
388 : retin = 1
389 11 : do
390 :
391 : ! compute the new centers, based on the updated memberships.
392 :
393 22 : call setCenter(membership, disq, sample, center, size, potential)
394 112 : failed = any(size < sizemin)
395 22 : if (failed) exit
396 :
397 : ! compute the new memberships, based on the input initial centers.
398 :
399 22 : call setMember(membersnew, disq, sample, center)
400 5909 : if (all(membersnew == membership)) exit
401 12 : failed = retinxam < retin
402 12 : if (failed) exit
403 :
404 12 : retin = retin + 1
405 :
406 : ! compute the new centers, based on the updated memberships.
407 :
408 12 : call setCenter(membersnew, disq, sample, center, size, potential)
409 72 : failed = any(size < sizemin)
410 12 : if (failed) exit
411 :
412 : ! compute the new memberships, based on the input initial centers.
413 :
414 12 : call setMember(membership, disq, sample, center)
415 10282 : if (all(membersnew == membership)) exit
416 11 : failed = retinxam < retin
417 11 : if (failed) exit
418 :
419 22 : retin = retin + 1
420 :
421 : end do
422 11 : if (present(niter)) niter = retin
423 :
424 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
425 : #elif setKmeans_ENABLED && (RNGF_ENABLED || RNGX_ENABLED)
426 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
427 :
428 : integer(IK) :: nfail_def, itry, ncls
429 :
430 10 : if (present(nfail)) then
431 0 : CHECK_ASSERTION(__LINE__, 0_IK <= nfail, SK_"@setKmeans(): The condition `0 <= nfail` must hold. nfail = "//getStr(nfail))
432 : nfail_def = nfail
433 : else
434 : nfail_def = 1
435 : end if
436 :
437 20 : ncls = ubound(center, 2, IK)
438 10 : do itry = 1, nfail_def
439 10 : call setKmeansPP(rng, membership, disq, sample, ncls)
440 10 : call setKmeans(membership, disq, sample, center, size, potential, failed, niter, maxniter, minsize)
441 10 : if (.not. failed) exit
442 : end do
443 :
444 : #else
445 : !%%%%%%%%%%%%%%%%%%%%%%%%
446 : #error "Unrecognized interface."
447 : !%%%%%%%%%%%%%%%%%%%%%%%%
448 : #endif
449 :
450 : #undef INTRINSIC_ENABLED
|
