1 SUBROUTINE KMNS(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D, &
2 &ITRAN, LIVE, ITER, WSS, IFAULT)
4 C ALGORITHM AS 136 APPL. STATIST. (1979) VOL.28, NO.1
6 C Divide M points in N-dimensional space into K clusters so that
7 C the within cluster sum of squares is minimized.
9 integer*4 :: I,IJ,IL,J,K,L,M,N,AA
10 integer*4 :: INDX,ITER,IFAULT,DA,DC
11 integer*4 :: IC1(M), IC2(M), NC(K), NCP(K), ITRAN(K), LIVE(K)
12 c real*8 :: A, C, D, AN1, AN2(K),WSS(K), DT(2)
13 real*8 :: BIG, ZERO, ONE, DB
15 real*8, dimension(M,N) :: A
16 real*8, dimension(K,N) :: C
17 real*8, dimension(M) :: D
18 real*8, dimension(K) :: AN1
19 real*8, dimension(K) :: AN2
20 real*8, dimension(K) :: WSS
21 real*8, dimension(2) :: DT
24 C Define BIG to be a very large positive number
26 DATA BIG /1.E30/, ZERO /0.0/, ONE /1.0/
29 IF (K .LE. 1 .OR. K .GE. M) RETURN
31 C For each point I, find its two closest centres, IC1(I) and
32 C IC2(I). Assign it to IC1(I).
41 DT(IL) = DT(IL) + DA*DA
43 IF (DT(1) .GT. DT(2)) THEN
55 IF (DB .GE. DT(2)) GO TO 50
57 IF (DB .LT. DT(1)) GO TO 40
67 C Update cluster centres to be the average of points contained
79 80 C(L,J) = C(L,J) + A(I,J)
82 C Check to see if there is any empty cluster at this stage
85 IF (NC(L) .EQ. 0) THEN
91 110 C(L,J) = C(L,J) / AA
93 C Initialize AN1, AN2, ITRAN & NCP
94 C AN1(L) = NC(L) / (NC(L) - 1)
95 C AN2(L) = NC(L) / (NC(L) + 1)
96 C ITRAN(L) = 1 if cluster L is updated in the quick-transfer stage,
98 C In the optimal-transfer stage, NCP(L) stores the step at which
99 C cluster L is last updated.
100 C In the quick-transfer stage, NCP(L) stores the step at which
101 C cluster L is last updated plus M.
103 AN2(L) = AA / (AA + ONE)
105 IF (AA .GT. ONE) AN1(L) = AA / (AA - ONE)
112 C In this stage, there is only one pass through the data. Each
113 C point is re-allocated, if necessary, to the cluster that will
114 C induce the maximum reduction in within-cluster sum of squares.
116 CALL OPTRA(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
119 C Stop if no transfer took place in the last M optimal transfer
122 IF (INDX .EQ. M) GO TO 150
124 C Each point is tested in turn to see if it should be re-allocated
125 C to the cluster to which it is most likely to be transferred,
126 C IC2(I), from its present cluster, IC1(I). Loop through the
127 C data until no further change is to take place.
129 CALL QTRAN(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
132 C If there are only two clusters, there is no need to re-enter the
133 C optimal transfer stage.
135 IF (K .EQ. 2) GO TO 150
137 C NCP has to be set to 0 before entering OPTRA.
143 C Since the specified number of iterations has been exceeded, set
144 C IFAULT = 2. This may indicate unforeseen looping.
148 C Compute within-cluster sum of squares for each cluster.
158 C(II,J) = C(II,J) + A(I,J)
162 180 C(L,J) = C(L,J) / FLOAT(NC(L))
165 DA = A(I,J) - C(II,J)
166 WSS(II) = WSS(II) + DA*DA
173 SUBROUTINE OPTRA(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
176 C ALGORITHM AS 136.1 APPL. STATIST. (1979) VOL.28, NO.1
178 C This is the optimal transfer stage.
180 C Each point is re-allocated, if necessary, to the cluster that
181 C will induce a maximum reduction in the within-cluster sum of
184 INTEGER IC1(M), IC2(M), NC(K), NCP(K), ITRAN(K), LIVE(K)
185 REAL A(M,N), D(M), C(K,N), AN1(K), AN2(K), ZERO, ONE
187 C Define BIG to be a very large positive number.
189 DATA BIG /1.0E30/, ZERO /0.0/, ONE/1.0/
191 C If cluster L is updated in the last quick-transfer stage, it
192 C belongs to the live set throughout this stage. Otherwise, at
193 C each step, it is not in the live set if it has not been updated
194 C in the last M optimal transfer steps.
197 IF (ITRAN(L) .EQ. 1) LIVE(L) = M + 1
205 C If point I is the only member of cluster L1, no transfer.
207 IF (NC(L1) .EQ. 1) GO TO 90
209 C If L1 has not yet been updated in this stage, no need to
212 IF (NCP(L1) .EQ. 0) GO TO 30
215 DF = A(I,J) - C(L1,J)
220 C Find the cluster with minimum R2.
224 DB = A(I,J) - C(L2,J)
230 C If I >= LIVE(L1), then L1 is not in the live set. If this is
231 C true, we only need to consider clusters that are in the live set
232 C for possible transfer of point I. Otherwise, we need to consider
233 C all possible clusters.
235 IF (I .GE. LIVE(L1) .AND. I .GE. LIVE(L) .OR. L .EQ. L1 .OR.
236 * L .EQ. LL) GO TO 60
242 IF (DC .GE. RR) GO TO 60
247 IF (R2 .LT. D(I)) GO TO 70
249 C If no transfer is necessary, L2 is the new IC2(I).
254 C Update cluster centres, LIVE, NCP, AN1 & AN2 for clusters L1 and
255 C L2, and update IC1(I) & IC2(I).
267 C(L1,J) = (C(L1,J) * AL1 - A(I,J)) / ALW
268 C(L2,J) = (C(L2,J) * AL2 + A(I,J)) / ALT
274 IF (ALW .GT. ONE) AN1(L1) = ALW / (ALW - ONE)
276 AN2(L2) = ALT / (ALT + ONE)
280 IF (INDX .EQ. M) RETURN
284 C ITRAN(L) = 0 before entering QTRAN. Also, LIVE(L) has to be
285 C decreased by M before re-entering OPTRA.
288 LIVE(L) = LIVE(L) - M
295 SUBROUTINE QTRAN(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
298 C ALGORITHM AS 136.2 APPL. STATIST. (1979) VOL.28, NO.1
300 C This is the quick transfer stage.
301 C IC1(I) is the cluster which point I belongs to.
302 C IC2(I) is the cluster which point I is most likely to be
304 C For each point I, IC1(I) & IC2(I) are switched, if necessary, to
305 C reduce within-cluster sum of squares. The cluster centres are
306 C updated after each step.
308 INTEGER IC1(M), IC2(M), NC(K), NCP(K), ITRAN(K)
309 REAL A(M,N), D(M), C(K,N), AN1(K), AN2(K), ZERO, ONE
311 C Define BIG to be a very large positive number
313 DATA BIG /1.0E30/, ZERO /0.0/, ONE /1.0/
315 C In the optimal transfer stage, NCP(L) indicates the step at which
316 C cluster L is last updated. In the quick transfer stage, NCP(L)
317 C is equal to the step at which cluster L is last updated plus M.
327 C If point I is the only member of cluster L1, no transfer.
329 IF (NC(L1) .EQ. 1) GO TO 60
331 C If ISTEP > NCP(L1), no need to re-compute distance from point I to
332 C cluster L1. Note that if cluster L1 is last updated exactly M
333 C steps ago, we still need to compute the distance from point I to
336 IF (ISTEP .GT. NCP(L1)) GO TO 30
339 DB = A(I,J) - C(L1,J)
344 C If ISTEP >= both NCP(L1) & NCP(L2) there will be no transfer of
345 C point I at this step.
347 30 IF (ISTEP .GE. NCP(L1) .AND. ISTEP .GE. NCP(L2)) GO TO 60
351 DE = A(I,J) - C(L2,J)
353 IF (DD .GE. R2) GO TO 60
356 C Update cluster centres, NCP, NC, ITRAN, AN1 & AN2 for clusters
357 C L1 & L2. Also update IC1(I) & IC2(I). Note that if any
358 C updating occurs in this stage, INDX is set back to 0.
371 C(L1,J) = (C(L1,J) * AL1 - A(I,J)) / ALW
372 C(L2,J) = (C(L2,J) * AL2 + A(I,J)) / ALT
378 IF (ALW .GT. ONE) AN1(L1) = ALW / (ALW - ONE)
380 AN2(L2) = ALT / (ALT + ONE)
384 C If no re-allocation took place in the last M steps, return.
386 60 IF (ICOUN .EQ. M) RETURN