source/ga/cluster.inc

   1       SUBROUTINE KMNS(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,    &
   2      &ITRAN, LIVE, ITER, WSS, IFAULT)
   3 C
   4 C     ALGORITHM AS 136  APPL. STATIST. (1979) VOL.28, NO.1
   5 C
   6 C     Divide M points in N-dimensional space into K clusters so that
   7 C     the within cluster sum of squares is minimized.
   8 C
   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
  14
  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
  22
  23 C
  24 C     Define BIG to be a very large positive number
  25 C
  26       DATA BIG /1.E30/, ZERO /0.0/, ONE /1.0/
  27 C
  28       IFAULT = 3
  29       IF (K .LE. 1 .OR. K .GE. M) RETURN
  30 C
  31 C     For each point I, find its two closest centres, IC1(I) and
  32 C     IC2(I).     Assign it to IC1(I).
  33 C
  34       DO 50 I = 1, M
  35         IC1(I) = 1
  36         IC2(I) = 2
  37         DO 10 IL = 1, 2
  38           DT(IL) = ZERO
  39           DO 10 J = 1, N
  40             DA = A(I,J) - C(IL,J)
  41             DT(IL) = DT(IL) + DA*DA
  42    10   CONTINUE
  43         IF (DT(1) .GT. DT(2)) THEN
  44           IC1(I) = 2
  45           IC2(I) = 1
  46           TEMP = DT(1)
  47           DT(1) = DT(2)
  48           DT(2) = TEMP
  49         END IF
  50         DO 50 L = 3, K
  51           DB = ZERO
  52           DO 30 J = 1, N
  53             DC = A(I,J) - C(L,J)
  54             DB = DB + DC*DC
  55             IF (DB .GE. DT(2)) GO TO 50
  56    30     CONTINUE
  57           IF (DB .LT. DT(1)) GO TO 40
  58           DT(2) = DB
  59           IC2(I) = L
  60           GO TO 50
  61    40     DT(2) = DT(1)
  62           IC2(I) = IC1(I)
  63           DT(1) = DB
  64           IC1(I) = L
  65    50 CONTINUE
  66 C
  67 C     Update cluster centres to be the average of points contained
  68 C     within them.
  69 C
  70       DO 70 L = 1, K
  71         NC(L) = 0
  72         DO 60 J = 1, N
  73    60   C(L,J) = ZERO
  74    70 CONTINUE
  75       DO 90 I = 1, M
  76         L = IC1(I)
  77         NC(L) = NC(L) + 1
  78         DO 80 J = 1, N
  79    80   C(L,J) = C(L,J) + A(I,J)
  80    90 CONTINUE
  81 C
  82 C     Check to see if there is any empty cluster at this stage
  83 C
  84       DO 120 L = 1, K
  85         IF (NC(L) .EQ. 0) THEN
  86           IFAULT = 1
  87           RETURN
  88         END IF
  89         AA = NC(L)
  90         DO 110 J = 1, N
  91   110   C(L,J) = C(L,J) / AA
  92 C
  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,
  97 C              = 0 otherwise
  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.
 102 C
 103         AN2(L) = AA / (AA + ONE)
 104         AN1(L) = BIG
 105         IF (AA .GT. ONE) AN1(L) = AA / (AA - ONE)
 106         ITRAN(L) = 1
 107         NCP(L) = -1
 108   120 CONTINUE
 109       INDX = 0
 110       DO 140 IJ = 1, ITER
 111 C
 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.
 115 C
 116         CALL OPTRA(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
 117      *        ITRAN, LIVE, INDX)
 118 C
 119 C     Stop if no transfer took place in the last M optimal transfer
 120 C     steps.
 121 C
 122         IF (INDX .EQ. M) GO TO 150
 123 C
 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.
 128 C
 129         CALL QTRAN(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
 130      *       ITRAN, INDX)
 131 C
 132 C     If there are only two clusters, there is no need to re-enter the
 133 C     optimal transfer stage.
 134 C
 135         IF (K .EQ. 2) GO TO 150
 136 C
 137 C     NCP has to be set to 0 before entering OPTRA.
 138 C
 139         DO 130 L = 1, K
 140   130   NCP(L) = 0
 141   140 CONTINUE
 142 C
 143 C     Since the specified number of iterations has been exceeded, set
 144 C     IFAULT = 2.   This may indicate unforeseen looping.
 145 C
 146       IFAULT = 2
 147 C
 148 C     Compute within-cluster sum of squares for each cluster.
 149 C
 150   150 DO 160 L = 1, K
 151         WSS(L) = ZERO
 152         DO 160 J = 1, N
 153           C(L,J) = ZERO
 154   160 CONTINUE
 155       DO 170 I = 1, M
 156         II = IC1(I)
 157         DO 170 J = 1, N
 158           C(II,J) = C(II,J) + A(I,J)
 159   170 CONTINUE
 160       DO 190 J = 1, N
 161         DO 180 L = 1, K
 162   180   C(L,J) = C(L,J) / FLOAT(NC(L))
 163         DO 190 I = 1, M
 164           II = IC1(I)
 165           DA = A(I,J) - C(II,J)
 166           WSS(II) = WSS(II) + DA*DA
 167   190 CONTINUE
 168 C
 169       RETURN
 170       END
 171 C
 172 C
 173       SUBROUTINE OPTRA(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
 174      *      ITRAN, LIVE, INDX)
 175 C
 176 C     ALGORITHM AS 136.1  APPL. STATIST. (1979) VOL.28, NO.1
 177 C
 178 C     This is the optimal transfer stage.
 179 C
 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
 182 C     squares.
 183 C
 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
 186 C
 187 C     Define BIG to be a very large positive number.
 188 C
 189       DATA BIG /1.0E30/, ZERO /0.0/, ONE/1.0/
 190 C
 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.
 195 C
 196       DO 10 L = 1, K
 197         IF (ITRAN(L) .EQ. 1) LIVE(L) = M + 1
 198    10 CONTINUE
 199       DO 100 I = 1, M
 200         INDX = INDX + 1
 201         L1 = IC1(I)
 202         L2 = IC2(I)
 203         LL = L2
 204 C
 205 C     If point I is the only member of cluster L1, no transfer.
 206 C
 207         IF (NC(L1) .EQ. 1) GO TO 90
 208 C
 209 C     If L1 has not yet been updated in this stage, no need to
 210 C     re-compute D(I).
 211 C
 212         IF (NCP(L1) .EQ. 0) GO TO 30
 213         DE = ZERO
 214         DO 20 J = 1, N
 215           DF = A(I,J) - C(L1,J)
 216           DE = DE + DF*DF
 217    20   CONTINUE
 218         D(I) = DE * AN1(L1)
 219 C
 220 C     Find the cluster with minimum R2.
 221 C
 222    30   DA = ZERO
 223         DO 40 J = 1, N
 224           DB = A(I,J) - C(L2,J)
 225           DA = DA + DB*DB
 226    40   CONTINUE
 227         R2 = DA * AN2(L2)
 228         DO 60 L = 1, K
 229 C
 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.
 234 C
 235           IF (I .GE. LIVE(L1) .AND. I .GE. LIVE(L) .OR. L .EQ. L1 .OR.
 236      *        L .EQ. LL) GO TO 60
 237           RR = R2 / AN2(L)
 238           DC = ZERO
 239           DO 50 J = 1, N
 240             DD = A(I,J) - C(L,J)
 241             DC = DC + DD*DD
 242             IF (DC .GE. RR) GO TO 60
 243    50     CONTINUE
 244           R2 = DC * AN2(L)
 245           L2 = L
 246    60     CONTINUE
 247           IF (R2 .LT. D(I)) GO TO 70
 248 C
 249 C     If no transfer is necessary, L2 is the new IC2(I).
 250 C
 251           IC2(I) = L2
 252           GO TO 90
 253 C
 254 C     Update cluster centres, LIVE, NCP, AN1 & AN2 for clusters L1 and
 255 C     L2, and update IC1(I) & IC2(I).
 256 C
 257    70     INDX = 0
 258           LIVE(L1) = M + I
 259           LIVE(L2) = M + I
 260           NCP(L1) = I
 261           NCP(L2) = I
 262           AL1 = NC(L1)
 263           ALW = AL1 - ONE
 264           AL2 = NC(L2)
 265           ALT = AL2 + ONE
 266           DO 80 J = 1, N
 267             C(L1,J) = (C(L1,J) * AL1 - A(I,J)) / ALW
 268             C(L2,J) = (C(L2,J) * AL2 + A(I,J)) / ALT
 269    80     CONTINUE
 270           NC(L1) = NC(L1) - 1
 271           NC(L2) = NC(L2) + 1
 272           AN2(L1) = ALW / AL1
 273           AN1(L1) = BIG
 274           IF (ALW .GT. ONE) AN1(L1) = ALW / (ALW - ONE)
 275           AN1(L2) = ALT / AL2
 276           AN2(L2) = ALT / (ALT + ONE)
 277           IC1(I) = L2
 278           IC2(I) = L1
 279    90   CONTINUE
 280         IF (INDX .EQ. M) RETURN
 281   100 CONTINUE
 282       DO 110 L = 1, K
 283 C
 284 C     ITRAN(L) = 0 before entering QTRAN.   Also, LIVE(L) has to be
 285 C     decreased by M before re-entering OPTRA.
 286 C
 287         ITRAN(L) = 0
 288         LIVE(L) = LIVE(L) - M
 289   110 CONTINUE
 290 C
 291       RETURN
 292       END
 293 C
 294 C
 295       SUBROUTINE QTRAN(A, M, N, C, K, IC1, IC2, NC, AN1, AN2, NCP, D,
 296      *    ITRAN, INDX)
 297 C
 298 C     ALGORITHM AS 136.2  APPL. STATIST. (1979) VOL.28, NO.1
 299 C
 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
 303 C         transferred to.
 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.
 307 C
 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
 310 C
 311 C     Define BIG to be a very large positive number
 312 C
 313       DATA BIG /1.0E30/, ZERO /0.0/, ONE /1.0/
 314 C
 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.
 318 C
 319       ICOUN = 0
 320       ISTEP = 0
 321    10 DO 70 I = 1, M
 322         ICOUN = ICOUN + 1
 323         ISTEP = ISTEP + 1
 324         L1 = IC1(I)
 325         L2 = IC2(I)
 326 C
 327 C     If point I is the only member of cluster L1, no transfer.
 328 C
 329         IF (NC(L1) .EQ. 1) GO TO 60
 330 C
 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
 334 C     cluster L1.
 335 C
 336         IF (ISTEP .GT. NCP(L1)) GO TO 30
 337         DA = ZERO
 338         DO 20 J = 1, N
 339           DB = A(I,J) - C(L1,J)
 340           DA = DA + DB*DB
 341    20   CONTINUE
 342         D(I) = DA * AN1(L1)
 343 C
 344 C     If ISTEP >= both NCP(L1) & NCP(L2) there will be no transfer of
 345 C     point I at this step.
 346 C
 347    30   IF (ISTEP .GE. NCP(L1) .AND. ISTEP .GE. NCP(L2)) GO TO 60
 348         R2 = D(I) / AN2(L2)
 349         DD = ZERO
 350         DO 40 J = 1, N
 351           DE = A(I,J) - C(L2,J)
 352           DD = DD + DE*DE
 353           IF (DD .GE. R2) GO TO 60
 354    40   CONTINUE
 355 C
 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.
 359 C
 360         ICOUN = 0
 361         INDX = 0
 362         ITRAN(L1) = 1
 363         ITRAN(L2) = 1
 364         NCP(L1) = ISTEP + M
 365         NCP(L2) = ISTEP + M
 366         AL1 = NC(L1)
 367         ALW = AL1 - ONE
 368         AL2 = NC(L2)
 369         ALT = AL2 + ONE
 370         DO 50 J = 1, N
 371           C(L1,J) = (C(L1,J) * AL1 - A(I,J)) / ALW
 372           C(L2,J) = (C(L2,J) * AL2 + A(I,J)) / ALT
 373    50   CONTINUE
 374         NC(L1) = NC(L1) - 1
 375         NC(L2) = NC(L2) + 1
 376         AN2(L1) = ALW / AL1
 377         AN1(L1) = BIG
 378         IF (ALW .GT. ONE) AN1(L1) = ALW / (ALW - ONE)
 379         AN1(L2) = ALT / AL2
 380         AN2(L2) = ALT / (ALT + ONE)
 381         IC1(I) = L2
 382         IC2(I) = L1
 383 C
 384 C     If no re-allocation took place in the last M steps, return.
 385 C
 386    60   IF (ICOUN .EQ. M) RETURN
 387    70 CONTINUE
 388       GO TO 10
 389       END